zdc_neutron/analysis/fwd_neutrons_geant.C

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

0002 int get_layer_number(TVector3 pos){
0003
0004        // Get Layer position wrt proton beam direction

0005     pos.RotateY(0.025);
0006
0007        auto local_z = pos.Z() - 35.8*1000.; //in mm

0008
0009        // Convert to layer number

0010        // local_z = 22.25 + (layer_number - 1)*24.9

0011        int layer_number = (int)std::round( (local_z - 22.25)/24.9 ) + 1;
0012
0013        return layer_number;
0014 }
0015
0016 //------------------

0017 void fwd_neutrons_geant(std::string inputfile, std::string outputfile){
0018
0019     //Define Style

0020     gStyle->SetOptStat(0);
0021     gStyle->SetPadBorderMode(0);
0022     gStyle->SetFrameBorderMode(0);
0023     gStyle->SetFrameLineWidth(2);
0024     gStyle->SetLabelSize(0.035,"X");
0025     gStyle->SetLabelSize(0.035,"Y");
0026     //gStyle->SetLabelOffset(0.01,"X");

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

0028     gStyle->SetTitleXSize(0.04);
0029     gStyle->SetTitleXOffset(0.9);
0030     gStyle->SetTitleYSize(0.04);
0031     gStyle->SetTitleYOffset(0.9);
0032
0033     //Define histograms

0034     TH2 *h1_neut = new TH2D("h1_neut","Neutron true energy vs. polar angle",100,0.6,2.2,100,0,200);
0035     h1_neut->GetXaxis()->SetTitle("#theta [deg]"); h1_neut->GetXaxis()->CenterTitle();
0036     h1_neut->GetYaxis()->SetTitle("E [GeV]"); h1_neut->GetYaxis()->CenterTitle();
0037
0038     TH2 *h2_neut = new TH2D("h2_neut","Neutron true azimuthal angle vs. polar angle around p axis",100,-0.1,12,100,-200,200);
0039     h2_neut->GetXaxis()->SetTitle("#theta^{*} [mRad]"); h2_neut->GetXaxis()->CenterTitle();
0040     h2_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h2_neut->GetYaxis()->CenterTitle();
0041
0042     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);
0043     h2a_neut->GetXaxis()->SetTitle("cos(#theta^{*})"); h2a_neut->GetXaxis()->CenterTitle();
0044     h2a_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h2a_neut->GetYaxis()->CenterTitle();
0045
0046     //Require HCal hit energy sum > 1 GeV

0047     TH2 *h3_neut = new TH2D("h3_neut","Neutron true azimuthal angle vs. polar angle around p axis",100,-0.1,12,100,-200,200);
0048     h3_neut->GetXaxis()->SetTitle("#theta^{*} [mRad]"); h3_neut->GetXaxis()->CenterTitle();
0049     h3_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h3_neut->GetYaxis()->CenterTitle();
0050
0051     TH2 *h4_neut = new TH2D("h4_neut","Neutron local hit position at ZDC HCal front face",100,-300,300,100,-300,300);
0052     h4_neut->GetXaxis()->SetTitle("x [mm]"); h4_neut->GetXaxis()->CenterTitle();
0053     h4_neut->GetYaxis()->SetTitle("y [mm]"); h4_neut->GetYaxis()->CenterTitle();
0054
0055     TH1 *h1_ecal = new TH1D("h1_ecal","Total true hit energy sum",100,-0.1,2);
0056     h1_ecal->GetXaxis()->SetTitle("E [GeV]"); h1_ecal->GetXaxis()->CenterTitle();
0057     h1_ecal->SetLineColor(kRed);h1_ecal->SetLineWidth(2);
0058
0059     TH1 *h1_hcal = new TH1D("h1_hcal","Total true hit energy sum",100,-0.1,2);
0060     h1_hcal->GetXaxis()->SetTitle("E [GeV]"); h1_hcal->GetXaxis()->CenterTitle();
0061     h1_hcal->SetLineColor(kBlue);h1_hcal->SetLineWidth(2);
0062
0063     //Cut on theta*

0064     TH1 *h2_ecal = new TH1D("h2_ecal","Total true hit energy sum",100,-0.1,2);
0065     h2_ecal->GetXaxis()->SetTitle("E [GeV]"); h2_ecal->GetXaxis()->CenterTitle();
0066     h2_ecal->SetLineColor(kRed);h2_ecal->SetLineWidth(2);
0067
0068     //Cut on theta*

0069     TH1 *h2_hcal = new TH1D("h2_hcal","Total true hit energy sum",100,-0.1,2);
0070     h2_hcal->GetXaxis()->SetTitle("E [GeV]"); h2_hcal->GetXaxis()->CenterTitle();
0071     h2_hcal->SetLineColor(kBlue);h2_hcal->SetLineWidth(2);
0072
0073     //Cut on theta*

0074     TH1 *h2a_ecal = new TH1D("h2a_ecal","Total true hit energy sum",200,-0.1,30);
0075     h2a_ecal->GetXaxis()->SetTitle("E [GeV]"); h2a_ecal->GetXaxis()->CenterTitle();
0076     h2a_ecal->SetLineColor(kRed);h2a_ecal->SetLineWidth(2);
0077
0078     //Cut on theta*

0079     TH1 *h2a_hcal = new TH1D("h2a_hcal","Total true hit energy sum",200,-0.1,30);
0080     h2a_hcal->GetXaxis()->SetTitle("E [GeV]"); h2a_hcal->GetXaxis()->CenterTitle();
0081     h2a_hcal->SetLineColor(kBlue);h2a_hcal->SetLineWidth(2);
0082
0083     //Cut on theta*

0084     TH2 *h2_cal_both = new TH2D("h2_cal_both","Total true hit energy sum: ECal vs. HCal",100,-0.1,2.5,100,-0.1,30);
0085     h2_cal_both->GetXaxis()->SetTitle("E_{HCal.} [GeV]"); h2_cal_both->GetXaxis()->CenterTitle();
0086     h2_cal_both->GetYaxis()->SetTitle("E_{ECal.} [GeV]"); h2_cal_both->GetYaxis()->CenterTitle();
0087
0088     TH1 *h3a_hcal = new TH1D("h3a_hcal","Cells w/ non-zero energy deposit ",100,0,64);
0089     h3a_hcal->GetXaxis()->SetTitle("Layer Number in ZDC HCal"); h3a_hcal->GetXaxis()->CenterTitle();
0090     h3a_hcal->GetYaxis()->SetTitle("Number of cells hit per event");h3a_hcal->GetYaxis()->CenterTitle();
0091     h3a_hcal->SetLineColor(kBlue);h3a_hcal->SetLineWidth(3);
0092
0093     //Cut on theta* < 3.5 mRad

0094     TH1 *h3b_hcal = new TH1D("h3b_hcal","Cells w/ non-zero energy deposit ",100,0,64);
0095     h3b_hcal->GetXaxis()->SetTitle("Layer Number in ZDC HCal"); h3b_hcal->GetXaxis()->CenterTitle();
0096     h3b_hcal->GetYaxis()->SetTitle("Number of cells hit per event");h3b_hcal->GetYaxis()->CenterTitle();
0097     h3b_hcal->SetLineColor(kBlue);h3b_hcal->SetLineWidth(3);
0098
0099     //Cut on 3.5 mRad < theta* < 6 mRad

0100     TH1 *h3c_hcal = new TH1D("h3c_hcal","Cells w/ non-zero energy deposit ",100,0,64);
0101     h3c_hcal->GetXaxis()->SetTitle("Layer Number in ZDC HCal"); h3c_hcal->GetXaxis()->CenterTitle();
0102     h3c_hcal->GetYaxis()->SetTitle("Number of cells hit per event");h3c_hcal->GetYaxis()->CenterTitle();
0103     h3c_hcal->SetLineColor(kBlue);h3c_hcal->SetLineWidth(3);
0104
0105     //Read ROOT file

0106     TFile* file = new TFile(inputfile.c_str());
0107     TTree *tree = (TTree*) file->Get("events");
0108
0109     //Set cut

0110     float edep_cut = 1.; //In GeV

0111
0112     cout<<"Total number of events to analyze is "<<tree->GetEntries()<<endl;
0113
0114     //Create Array Reader

0115     TTreeReader tr(tree);
0116
0117     TTreeReaderArray<int>   gen_status(tr, "MCParticles.generatorStatus");
0118     TTreeReaderArray<int>   gen_pid(tr, "MCParticles.PDG");
0119     TTreeReaderArray<float> gen_px(tr, "MCParticles.momentum.x");
0120     TTreeReaderArray<float> gen_py(tr, "MCParticles.momentum.y");
0121     TTreeReaderArray<float> gen_pz(tr, "MCParticles.momentum.z");
0122     TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass");
0123     TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
0124     TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
0125     TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
0126     TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
0127
0128     //ZDC LYSO ECal

0129     TTreeReaderArray<float> ecal_hit_e(tr,"EcalFarForwardZDCHits.energy");
0130
0131     //ZDC SiPM-on-tile HCal

0132     TTreeReaderArray<float> hcal_hit_e(tr, "HcalFarForwardZDCHits.energy");
0133     TTreeReaderArray<float> hcal_hit_x(tr, "HcalFarForwardZDCHits.position.x");
0134     TTreeReaderArray<float> hcal_hit_y(tr, "HcalFarForwardZDCHits.position.y");
0135     TTreeReaderArray<float> hcal_hit_z(tr, "HcalFarForwardZDCHits.position.z");
0136
0137     //Other variables

0138     int counter(0),counter_3p5(0),counter_3p5to6(0);
0139
0140     TLorentzVector neut_true; //True neutron in lab coordinates

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

0142
0143     float ecal_e_tot(0);
0144     float hcal_e_tot(0);
0145
0146     //Loop over events

0147     while (tr.Next()) {
0148
0149     if(counter%1000==0) cout<<"Analyzing event "<<counter<<endl;
0150     counter++;
0151
0152        //Reset variables

0153        ecal_e_tot = 0;
0154        hcal_e_tot = 0;
0155
0156        //Loop over generated particles, select primary neutron

0157        for(int igen=0;igen<gen_status.GetSize();igen++){
0158             if(gen_status[igen]==1 && gen_pid[igen]==2112){
0159
0160                      neut_true.SetXYZM(gen_px[igen],gen_py[igen],gen_pz[igen],gen_mass[igen]);
0161                      h1_neut->Fill(neut_true.Theta()*TMath::RadToDeg(),neut_true.E());
0162
0163             //Wrt proton beam direction

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

0167                      h2a_neut->Fill(std::cos(neut_true_rot.Theta()),neut_true_rot.Phi()*TMath::RadToDeg());
0168
0169               }
0170        } //End loop over generated particles

0171
0172        //Loop over ECal hits

0173        for(int ihit=0;ihit<ecal_hit_e.GetSize();ihit++){
0174                ecal_e_tot +=  ecal_hit_e[ihit];
0175        }
0176
0177        //Loop over HCal hits

0178        for(int ihit=0;ihit<hcal_hit_e.GetSize();ihit++){
0179               hcal_e_tot +=  hcal_hit_e[ihit];
0180
0181               //Cell hit position in global coordinates

0182               double hit_x = (double) hcal_hit_x[ihit];
0183               double hit_y = (double) hcal_hit_y[ihit];
0184               double hit_z = (double) hcal_hit_z[ihit];
0185
0186               TVector3 hit_pos_det(hit_x,hit_y,hit_z);
0187               auto layer_number = get_layer_number(hit_pos_det);
0188
0189               h3a_hcal->Fill(layer_number);
0190               if( neut_true_rot.Theta()*1000. < 3.5 ) h3b_hcal->Fill(layer_number);
0191               if( neut_true_rot.Theta()*1000. > 3.5 &&  neut_true_rot.Theta()*1000 < 6.) h3c_hcal->Fill(layer_number);
0192
0193        } //End loop over HCal hits

0194
0195        //Fill histograms

0196        h1_ecal->Fill(ecal_e_tot);
0197        h1_hcal->Fill(hcal_e_tot);
0198
0199        if(hcal_e_tot > edep_cut){
0200               h3_neut->Fill(neut_true_rot.Theta()*1000.,neut_true_rot.Phi()*TMath::RadToDeg()); //Theta* in mRad

0201
0202               //Projection to front face of HCal

0203               auto proj_x = 35.8 * 1000. * tan(neut_true_rot.Theta()) * cos(neut_true_rot.Phi()) ;
0204               auto proj_y = 35.8 * 1000. * tan(neut_true_rot.Theta()) * sin(neut_true_rot.Phi()) ;
0205
0206               h4_neut->Fill(proj_x,proj_y);
0207        }
0208
0209        if( neut_true_rot.Theta()*1000. < 3.5 ){
0210               h2_ecal->Fill(ecal_e_tot);
0211               h2_hcal->Fill(hcal_e_tot);
0212               h2a_ecal->Fill(ecal_e_tot);
0213               h2a_hcal->Fill(hcal_e_tot);
0214
0215               h2_cal_both->Fill(hcal_e_tot,ecal_e_tot);
0216
0217               counter_3p5++;
0218        }
0219
0220        if( neut_true_rot.Theta()*1000. > 3.5 &&  neut_true_rot.Theta()*1000 < 6.){
0221               counter_3p5to6++;
0222        }
0223
0224     } //End loop over events

0225
0226     //Scale histograms

0227     h3a_hcal->Scale(1./counter);
0228     h3b_hcal->Scale(1./counter_3p5);
0229     h3c_hcal->Scale(1./counter_3p5to6);
0230
0231     //Make plots

0232     TCanvas *c1 = new TCanvas("c1");
0233     h1_neut->Draw("colz");
0234
0235     TCanvas *c2 = new TCanvas("c2");
0236     h2_neut->Draw("colz");
0237
0238     TCanvas *c2a = new TCanvas("c2a");
0239     h2a_neut->Draw("colz");
0240
0241     TCanvas *c3 = new TCanvas("c3");
0242     h1_ecal->Draw("");
0243     h1_hcal->Draw("same");
0244
0245     TLegend *leg3 = new TLegend(0.6,0.6,0.85,0.8);
0246     leg3->SetBorderSize(0);leg3->SetFillStyle(0);
0247     leg3->AddEntry(h1_ecal,"Sum of ZDC Ecal hit energies","l");
0248     leg3->AddEntry(h1_hcal,"Sum of ZDC Hcal hit energies","l");
0249     leg3->Draw();
0250
0251     TCanvas *c4 = new TCanvas("c4");
0252     h3_neut->Draw("colz");
0253
0254     TLatex *tex4 = new TLatex(7,150,Form("Hcal hit energy sum > %.1f GeV",edep_cut));
0255     tex4->SetTextSize(0.035);
0256     tex4->Draw();
0257
0258     TCanvas *c4a = new TCanvas("c4a");
0259     h4_neut->Draw("colz");
0260
0261     TLatex *tex4a = new TLatex(50,250,Form("Hcal hit energy sum > %.1f GeV",edep_cut));
0262     tex4a->SetTextSize(0.035);
0263     tex4a->Draw();
0264
0265     TCanvas *c5 = new TCanvas("c5");
0266     h2_ecal->Draw("");
0267     h2_hcal->Draw("same");
0268
0269     TLegend *leg5 = new TLegend(0.5,0.6,0.9,0.8);
0270     leg5->SetBorderSize(0);leg5->SetFillStyle(0);
0271     leg5->SetHeader("Require neutron #theta^{*} (wrt proton beam) < 3.5 mRad");
0272     leg5->AddEntry(h1_ecal,"Sum of ZDC Ecal hit energies","l");
0273     leg5->AddEntry(h1_hcal,"Sum of ZDC Hcal hit energies","l");
0274     leg5->Draw();
0275
0276     TCanvas *c5a = new TCanvas("c5a");
0277     c5a->SetLogy();
0278     h2a_ecal->Draw("");
0279     h2a_hcal->Draw("same");
0280
0281     leg5->Draw();
0282
0283     TCanvas *c5b = new TCanvas("c5b");
0284     c5b->SetLogz();
0285     h2_cal_both->Draw("colz");
0286
0287     TLatex *tex5b = new TLatex(0,28,"Require neutron #theta^{*} (wrt proton beam) < 3.5 mRad");
0288     tex5b->SetTextSize(0.03);
0289     tex5b->Draw();
0290
0291     TCanvas *c6a = new TCanvas("c6a");
0292     h3a_hcal->Draw();
0293
0294     TLegend *leg6a = new TLegend(0.5,0.6,0.9,0.8);
0295     leg6a->SetBorderSize(0);leg6a->SetFillStyle(0);
0296     leg6a->SetHeader("Average over all events");
0297     leg6a->Draw();
0298
0299     TCanvas *c6b = new TCanvas("c6b");
0300     h3b_hcal->Draw();
0301
0302     TLegend *leg6b = new TLegend(0.5,0.6,0.9,0.8);
0303     leg6b->SetBorderSize(0);leg6b->SetFillStyle(0);
0304     leg6b->SetHeader("Average over events w/ neutron #theta^{*} < 3.5 mRad");
0305     leg6b->Draw();
0306
0307     TCanvas *c6c = new TCanvas("c6c");
0308     h3c_hcal->Draw();
0309
0310     TLegend *leg6c = new TLegend(0.5,0.6,0.9,0.8);
0311     leg6c->SetBorderSize(0);leg6c->SetFillStyle(0);
0312     leg6c->SetHeader("Average over events w/ neutron 3.5 < #theta^{*} < 6 mRad");
0313     leg6c->Draw();
0314
0315     //Print plots to file

0316     c1->Print(Form("%s[",outputfile.c_str()));
0317     c1->Print(Form("%s",outputfile.c_str()));
0318     c2->Print(Form("%s",outputfile.c_str()));
0319     c2a->Print(Form("%s",outputfile.c_str()));
0320     c3->Print(Form("%s",outputfile.c_str()));
0321     c4->Print(Form("%s",outputfile.c_str()));
0322     c4a->Print(Form("%s",outputfile.c_str()));
0323     c5->Print(Form("%s",outputfile.c_str()));
0324     c5a->Print(Form("%s",outputfile.c_str()));
0325     c5b->Print(Form("%s",outputfile.c_str()));
0326     c6a->Print(Form("%s",outputfile.c_str()));
0327     c6b->Print(Form("%s",outputfile.c_str()));
0328     c6c->Print(Form("%s",outputfile.c_str()));
0329     c6c->Print(Form("%s]",outputfile.c_str()));
0330
0331 }