zdc_neutron/analysis/fwd_neutrons_recon.C

0001 //------------------

0002 void fwd_neutrons_recon(std::string inputfile, std::string outputfile){
0003
0004     //Define Style

0005     gStyle->SetOptStat(0);
0006     gStyle->SetPadBorderMode(0);
0007     gStyle->SetFrameBorderMode(0);
0008     gStyle->SetFrameLineWidth(2);
0009     gStyle->SetLabelSize(0.035,"X");
0010     gStyle->SetLabelSize(0.035,"Y");
0011     //gStyle->SetLabelOffset(0.01,"X");

0012     //gStyle->SetLabelOffset(0.01,"Y");

0013     gStyle->SetTitleXSize(0.04);
0014     gStyle->SetTitleXOffset(0.9);
0015     gStyle->SetTitleYSize(0.04);
0016     gStyle->SetTitleYOffset(0.9);
0017
0018     //Define histograms

0019     TH2 *h1_neut = new TH2D("h1_neut","Neutron true energy vs. polar angle",100,0.6,2.2,100,0,200);
0020     h1_neut->GetXaxis()->SetTitle("#theta [deg]"); h1_neut->GetXaxis()->CenterTitle();
0021     h1_neut->GetYaxis()->SetTitle("E [GeV]"); h1_neut->GetYaxis()->CenterTitle();
0022
0023     TH2 *h2_neut = new TH2D("h2_neut","Neutron true azimuthal angle vs. polar angle around p axis",100,-0.1,12,100,-200,200);
0024     h2_neut->GetXaxis()->SetTitle("#theta^{*} [mRad]"); h2_neut->GetXaxis()->CenterTitle();
0025     h2_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h2_neut->GetYaxis()->CenterTitle();
0026
0027     TH2 *h2a_neut = new TH2D("h2a_neut","Neutron true azimuthal angle vs. cosine of polar angle around p axis",100,0.99996,1,100,-200,200);
0028     h2a_neut->GetXaxis()->SetTitle("cos(#theta^{*})"); h2a_neut->GetXaxis()->CenterTitle();
0029     h2a_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h2a_neut->GetYaxis()->CenterTitle();
0030
0031     TH1 *h1_ecal_adc = new TH1D("h1_ecal_adc","ECal ADC amplitude spectrum",1000,-0.1,35000);
0032     h1_ecal_adc->GetXaxis()->SetTitle("ADC Channel"); h1_ecal_adc->GetXaxis()->CenterTitle();
0033     h1_ecal_adc->SetLineColor(kRed);h1_ecal_adc->SetLineWidth(2);
0034
0035     TH1 *h1_hcal_adc = new TH1D("h1_hcal_adc","HCal ADC amplitude spectrum",1000,-0.1,35000);
0036     h1_hcal_adc->GetXaxis()->SetTitle("ADC Channel"); h1_hcal_adc->GetXaxis()->CenterTitle();
0037     h1_hcal_adc ->SetLineColor(kBlue);h1_hcal_adc->SetLineWidth(2);
0038
0039     TH1 *h1_ecal = new TH1D("h1_ecal","Total reconstructed hit energy sum",100,-0.1,2);
0040     h1_ecal->GetXaxis()->SetTitle("E [GeV]"); h1_ecal->GetXaxis()->CenterTitle();
0041     h1_ecal->SetLineColor(kRed);h1_ecal->SetLineWidth(2);
0042
0043     TH1 *h1_hcal = new TH1D("h1_hcal","Total reconstructed hit energy sum",100,-0.1,2);
0044     h1_hcal->GetXaxis()->SetTitle("E [GeV]"); h1_hcal->GetXaxis()->CenterTitle();
0045     h1_hcal->SetLineColor(kBlue);h1_hcal->SetLineWidth(2);
0046
0047     //Cut on theta*

0048     TH1 *h2_ecal = new TH1D("h2_ecal","Total reconstructed hit energy sum",100,-0.1,2);
0049     h2_ecal->GetXaxis()->SetTitle("E [GeV]"); h2_ecal->GetXaxis()->CenterTitle();
0050     h2_ecal->SetLineColor(kRed);h2_ecal->SetLineWidth(2);
0051
0052     //Cut on theta*

0053     TH1 *h2_hcal = new TH1D("h2_hcal","Total reconstructed hit energy sum",100,-0.1,2);
0054     h2_hcal->GetXaxis()->SetTitle("E [GeV]"); h2_hcal->GetXaxis()->CenterTitle();
0055     h2_hcal->SetLineColor(kBlue);h2_hcal->SetLineWidth(2);
0056
0057     //Cut on theta*

0058     TH1 *h3_ecal = new TH1D("h3_ecal","Number of reconstructed clusters in ECal",5,0,5);
0059     h3_ecal->GetXaxis()->SetTitle("Number of clusters"); h3_ecal->GetXaxis()->CenterTitle();
0060     h3_ecal->GetXaxis()->SetNdivisions(405);h3_ecal->GetXaxis()->CenterLabels();
0061     h3_ecal->SetLineColor(kRed);h3_ecal->SetLineWidth(2);
0062
0063     //Cut on theta*

0064     TH1 *h3_hcal = new TH1D("h3_hcal","Number of reconstructed clusters in HCal",5,0,5);
0065     h3_hcal->GetXaxis()->SetTitle("Number of clusters"); h3_hcal->GetXaxis()->CenterTitle();
0066     h3_hcal->GetXaxis()->SetNdivisions(405);h3_hcal->GetXaxis()->CenterLabels();
0067     h3_hcal->SetLineColor(kBlue);h3_hcal->SetLineWidth(2);
0068
0069     //Cut on theta*

0070     TH1 *h4_hcal = new TH1D("h4_hcal","HCal cluster energy: Events w/ 1 HCal Cluster",100,-0.1,120);
0071     h4_hcal->GetXaxis()->SetTitle("E [GeV]"); h4_hcal->GetXaxis()->CenterTitle();
0072     h4_hcal->SetLineColor(kBlue);h4_hcal->SetLineWidth(2);
0073
0074     //Cut on generated neutron theta* + 1 HCal Cluster

0075     TH1 *h1_neut_rec = new TH1D("h1_neut_rec","Reconstructed Neutron Energy",100,-0.1,120);
0076     h1_neut_rec->GetXaxis()->SetTitle("E [GeV]"); h1_neut_rec->GetXaxis()->CenterTitle();
0077     h1_neut_rec->SetLineColor(kBlue);h1_neut_rec->SetLineWidth(2);
0078
0079     //Cut on generated neutron theta* + 1 HCal Cluster

0080     TH1 *h2_neut_rec = new TH2D("h2_neut_rec","Reconstructed Neutron local hit position at ZDC HCal front face",100,-200,200,100,-200,200);
0081     h2_neut_rec->GetXaxis()->SetTitle("x [mm]"); h2_neut_rec->GetXaxis()->CenterTitle();
0082     h2_neut_rec->GetYaxis()->SetTitle("y [mm]"); h2_neut_rec->GetYaxis()->CenterTitle();
0083
0084     //Read ROOT file

0085     TFile* file = new TFile(inputfile.c_str());
0086     TTree *tree = (TTree*) file->Get("events");
0087
0088     cout<<"Total number of events to analyze is "<<tree->GetEntries()<<endl;
0089
0090     //Create Array Reader

0091     TTreeReader tr(tree);
0092
0093     TTreeReaderArray<int>   gen_status(tr, "MCParticles.generatorStatus");
0094     TTreeReaderArray<int>   gen_pid(tr, "MCParticles.PDG");
0095     TTreeReaderArray<float> gen_px(tr, "MCParticles.momentum.x");
0096     TTreeReaderArray<float> gen_py(tr, "MCParticles.momentum.y");
0097     TTreeReaderArray<float> gen_pz(tr, "MCParticles.momentum.z");
0098     TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass");
0099     TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
0100     TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
0101     TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
0102     TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
0103
0104     //ZDC LYSO ECal

0105     TTreeReaderArray<int> ecal_adc(tr,"EcalFarForwardZDCRawHits.amplitude");
0106     TTreeReaderArray<float> ecal_hit_e(tr,"EcalFarForwardZDCRecHits.energy");
0107     TTreeReaderArray<float> ecal_cluster_energy(tr, "EcalFarForwardZDCClusters.energy");
0108
0109     //ZDC SiPM-on-tile HCal

0110     TTreeReaderArray<int> hcal_adc(tr,"HcalFarForwardZDCRawHits.amplitude");
0111     TTreeReaderArray<float> hcal_hit_e(tr, "HcalFarForwardZDCRecHits.energy");
0112     TTreeReaderArray<float> hcal_cluster_energy(tr, "HcalFarForwardZDCClusters.energy");
0113
0114     //Reconstructed neutron quantity

0115     TTreeReaderArray<float> rec_neutron_energy(tr,"ReconstructedFarForwardZDCNeutrons.energy");
0116     TTreeReaderArray<float> rec_neutron_px(tr,"ReconstructedFarForwardZDCNeutrons.momentum.x");
0117     TTreeReaderArray<float> rec_neutron_py(tr,"ReconstructedFarForwardZDCNeutrons.momentum.y");
0118     TTreeReaderArray<float> rec_neutron_pz(tr,"ReconstructedFarForwardZDCNeutrons.momentum.z");
0119
0120     //Other variables

0121     int counter(0);
0122
0123     TLorentzVector neut_true; //True neutron in lab coordinates

0124     TLorentzVector neut_true_rot; //True neutron wrt proton beam direction

0125
0126     float ecal_e_tot(0);
0127     float hcal_e_tot(0);
0128
0129     //Loop over events

0130     while (tr.Next()) {
0131
0132     if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
0133     counter++;
0134
0135        //Reset variables

0136        ecal_e_tot = 0;
0137        hcal_e_tot = 0;
0138
0139        //Loop over generated particles, select primary neutron

0140        for(int igen=0;igen<gen_status.GetSize();igen++){
0141             if(gen_status[igen]==1 && gen_pid[igen]==2112){
0142
0143                      neut_true.SetXYZM(gen_px[igen],gen_py[igen],gen_pz[igen],gen_mass[igen]);
0144                      h1_neut->Fill(neut_true.Theta()*TMath::RadToDeg(),neut_true.E());
0145
0146             //Wrt proton beam direction

0147             neut_true_rot = neut_true;
0148             neut_true_rot.RotateY(0.025);
0149             h2_neut->Fill(neut_true_rot.Theta()*1000.,neut_true_rot.Phi()*TMath::RadToDeg()); //Theta* in mRad

0150                         h2a_neut->Fill(std::cos(neut_true_rot.Theta()),neut_true_rot.Phi()*TMath::RadToDeg());
0151
0152               }
0153        } //End loop over generated particles

0154
0155        //Loop over ECal ADC hits (Raw hits may have different length than Rec hits due to zero supression)

0156        for(int iadc=0;iadc<ecal_adc.GetSize();iadc++){
0157                h1_ecal_adc->Fill(ecal_adc[iadc]);
0158        }
0159
0160        //Loop over ECal hits

0161        for(int ihit=0;ihit<ecal_hit_e.GetSize();ihit++){
0162                ecal_e_tot +=  ecal_hit_e[ihit];
0163        }
0164
0165        //Loop over ECal ADC hits (Raw hits may have different length than Rec hits due to zero supression)

0166        for(int iadc=0;iadc<hcal_adc.GetSize();iadc++){
0167                h1_hcal_adc->Fill(hcal_adc[iadc]);
0168        }
0169
0170        //Loop over HCal hits

0171        for(int ihit=0;ihit<hcal_hit_e.GetSize();ihit++){
0172                hcal_e_tot +=  hcal_hit_e[ihit];
0173        }
0174
0175        //ECal cluster size

0176        int ecal_clus_size = ecal_cluster_energy.GetSize();
0177
0178        //HCal cluster size

0179        int hcal_clus_size = hcal_cluster_energy.GetSize();
0180
0181        //Fill histograms for total energy and clusters

0182        h1_ecal->Fill(ecal_e_tot);
0183        h1_hcal->Fill(hcal_e_tot);
0184
0185        if( neut_true_rot.Theta()*1000. < 3.5 ){
0186                 h2_ecal->Fill(ecal_e_tot);
0187                 h2_hcal->Fill(hcal_e_tot);
0188
0189                 h3_ecal->Fill(ecal_clus_size);
0190                 h3_hcal->Fill(hcal_clus_size);
0191
0192                 //HCal cluster energy -- 1 cluster events

0193                 if(hcal_clus_size==1) h4_hcal->Fill(hcal_cluster_energy[0]);
0194        }
0195
0196        //Reconstructed neutron(s)

0197        for(int irec=0;irec<rec_neutron_energy.GetSize();irec++){
0198                 if(neut_true_rot.Theta()*1000. < 3.5 && hcal_clus_size==1){
0199                         h1_neut_rec->Fill(rec_neutron_energy[irec]);
0200
0201                         TVector3 neut_rec(rec_neutron_px[irec],rec_neutron_py[irec],rec_neutron_pz[irec]);
0202
0203                         //Wrt proton beam direction

0204             TVector3 neut_rec_rot = neut_rec;
0205             neut_rec_rot.RotateY(0.025);
0206
0207                         //Projection of reconstructed neutron to front face of HCal

0208                         auto proj_x = 35.8 * 1000. * tan(neut_rec_rot.Theta()) * cos(neut_rec_rot.Phi()) ;
0209                         auto proj_y = 35.8 * 1000. * tan(neut_rec_rot.Theta()) * sin(neut_rec_rot.Phi()) ;
0210                         h2_neut_rec->Fill(proj_x,proj_y);
0211                 }
0212        } //End loop over reconstructed neutrons

0213
0214     } //End loop over events

0215
0216     //Make plots

0217     TCanvas *c1 = new TCanvas("c1");
0218     h1_neut->Draw("colz");
0219
0220     TCanvas *c2 = new TCanvas("c2");
0221     h2_neut->Draw("colz");
0222
0223     TCanvas *c2a = new TCanvas("c2a");
0224     h2a_neut->Draw("colz");
0225
0226     TCanvas *c3a = new TCanvas("c3a");
0227     c3a->SetLogy();
0228     h1_ecal_adc->Draw();
0229
0230     TCanvas *c3b = new TCanvas("c3b");
0231     c3b->SetLogy();
0232     h1_hcal_adc->Draw();
0233
0234     TCanvas *c4 = new TCanvas("c4");
0235     h1_ecal->Draw("");
0236     h1_hcal->Draw("same");
0237
0238     TLegend *leg4 = new TLegend(0.6,0.6,0.85,0.8);
0239     leg4->SetBorderSize(0);leg4->SetFillStyle(0);
0240     leg4->AddEntry(h1_ecal,"Sum of digitized ZDC Ecal hit energies","l");
0241     leg4->AddEntry(h1_hcal,"Sum of digitized ZDC Hcal hit energies","l");
0242     leg4->Draw();
0243
0244     TCanvas *c5 = new TCanvas("c5");
0245     h2_ecal->Draw("");
0246     h2_hcal->Draw("same");
0247
0248     TLegend *leg5 = new TLegend(0.5,0.6,0.9,0.8);
0249     leg5->SetBorderSize(0);leg5->SetFillStyle(0);
0250     leg5->SetHeader("Require neutron #theta^{*} (wrt proton beam) < 3.5 mRad");
0251     leg5->AddEntry(h1_ecal,"Sum of digitized ZDC Ecal hit energies","l");
0252     leg5->AddEntry(h1_hcal,"Sum of digitized ZDC Hcal hit energies","l");
0253     leg5->Draw();
0254
0255     TCanvas *c6a = new TCanvas("c6a");
0256     h3_ecal->Draw();
0257
0258     TLegend *leg6 = new TLegend(0.5,0.7,0.9,0.9);
0259     leg6->SetBorderSize(0);leg6->SetFillStyle(0);
0260     leg6->SetHeader("Require neutron #theta^{*} (wrt proton beam) < 3.5 mRad");
0261     leg6->Draw();
0262
0263     TCanvas *c6b = new TCanvas("c6b");
0264     h3_hcal->Draw();
0265     leg6->Draw();
0266
0267     TCanvas *c7 = new TCanvas("c7");
0268     h4_hcal->Draw();
0269
0270     TLegend *leg7 = new TLegend(0.15,0.7,0.5,0.9);
0271     leg7->SetBorderSize(0);leg7->SetFillStyle(0);
0272     leg7->SetHeader("Require neutron #theta^{*} (wrt proton beam) < 3.5 mRad");
0273     leg7->Draw();
0274
0275     TCanvas *c8 = new TCanvas("c8");
0276     h1_neut_rec->Draw();
0277
0278     TLegend *leg8 = new TLegend(0.15,0.7,0.7,0.9);
0279     leg8->SetBorderSize(0);leg8->SetFillStyle(0);
0280     leg8->SetHeader("Generated neutron #theta^{*} < 3.5 mRad + 1 HCal cluster");
0281     leg8->Draw();
0282
0283     TCanvas *c9 = new TCanvas("c9");
0284     h2_neut_rec->Draw("colz");
0285
0286     TLatex *tex9 = new TLatex(-150,150,"Generated neutron #theta^{*} < 3.5 mRad + 1 HCal cluster");
0287     tex9->SetTextSize(0.03);
0288     tex9->Draw();
0289
0290     //Print plots to file

0291     c1->Print(Form("%s[",outputfile.c_str()));
0292     c1->Print(Form("%s",outputfile.c_str()));
0293     c2->Print(Form("%s",outputfile.c_str()));
0294     c2a->Print(Form("%s",outputfile.c_str()));
0295     c3a->Print(Form("%s",outputfile.c_str()));
0296     c3b->Print(Form("%s",outputfile.c_str()));
0297     c4->Print(Form("%s",outputfile.c_str()));
0298     c5->Print(Form("%s",outputfile.c_str()));
0299     c6a->Print(Form("%s",outputfile.c_str()));
0300     c6b->Print(Form("%s",outputfile.c_str()));
0301     c7->Print(Form("%s",outputfile.c_str()));
0302     c8->Print(Form("%s",outputfile.c_str()));
0303     c9->Print(Form("%s",outputfile.c_str()));
0304     c9->Print(Form("%s]",outputfile.c_str()));
0305
0306 }