analysis/FrawleyTrackEvalMacros/track_calorimeter_matching.C

0001 #include <TChain.h>
0002 #include <TCanvas.h>
0003 #include <TH1.h>
0004 #include <TH1F.h>
0005 #include <TH2D.h>
0006 #include <TF1.h>
0007 #include <TFile.h>
0008 #include <TGraph.h>
0009 #include <TStyle.h>
0010 #include <TROOT.h>
0011 #include <TLegend.h>
0012 #include <TLine.h>
0013 #include <TLatex.h>
0014 #include <TRandom1.h>
0015 #include <TPolyLine.h>
0016 #include <iostream>
0017 #include <fstream>
0018 #include <TMath.h>
0019 #include <TLorentzVector.h>
0020 #include <TEfficiency.h>
0021 #include <TGraphAsymmErrors.h>
0022
0023 void track_calorimeter_matching()
0024 {
0025   int verbose = 2;
0026
0027   int n_maps_layer = 3;
0028   int n_intt_layer = 4;
0029   int n_tpc_layer = 40;
0030
0031   int nlayers = n_maps_layer+n_intt_layer+n_tpc_layer;  // maximum number of tracking layers for tpc+intt+maps
0032   static const int nmissed = 30;  // maximum number of missed layers to allow for a track
0033
0034   int embed_flag = 1;
0035
0036  double quality_cut = 3.0;
0037   double dca_cut = 0.1; //  cm
0038
0039   TH2D * hEoverP_cemc[2];
0040   hEoverP_cemc[0] = new TH2D("hEoverP0_cemc","",200, 0, 20, 200, 0, 2.0);
0041   hEoverP_cemc[0]->SetMarkerStyle(20);
0042   hEoverP_cemc[0]->SetMarkerSize(0.5);
0043   hEoverP_cemc[1] = new TH2D("hEoverP1_cemc","",200, 0, 20, 200, 0, 2.0);
0044   hEoverP_cemc[1]->SetMarkerStyle(20);
0045   hEoverP_cemc[1]->SetMarkerSize(0.5);
0046
0047
0048   int n_embed_rtrack = 0;
0049   int naccept = 0;
0050
0051  // The condor job output files
0052   for(int i=0;i<1000;i++)
0053     //for(int i=0;i<2;i++)
0054     {
0055       // Open the evaluator output file
0056
0057       TChain* ntp_track = new TChain("ntp_track","reco tracks");
0058       TChain* ntp_gtrack = new TChain("ntp_gtrack","g4 tracks");
0059       TChain* ntp_vertex = new TChain("ntp_vertex","events");
0060       TChain *ntp_cluster = new TChain("ntp_cluster","clusters");
0061
0062       // This include file contains the definitions of the ntuple variables, and the chain definitions
0063 #include "ntuple_variables.C"
0064
0065       char name[500];
0066       sprintf(name,"/sphenix/user/frawley/latest/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",i);
0067
0068       cout << "Adding file number " << i << " with name " << name << endl;
0069
0070       ntp_vertex->Add(name);
0071       ntp_track->Add(name);
0072       ntp_gtrack->Add(name);
0073
0074       // skip this file if there are no tracks
0075       if(!ntp_gtrack->GetEntries())
0076     continue;
0077
0078       if(verbose> 0)
0079     cout <<  " ntp_vertex entries: " << ntp_vertex->GetEntries()
0080          << " ntp_gtrack entries: " << ntp_gtrack->GetEntries()
0081          << " ntp_track entries: " << ntp_track->GetEntries()
0082          << endl;
0083
0084       //==================================
0085       // Begin event loop
0086       //==================================
0087
0088       // These keep track of the starting position in the track ntuples for each event
0089       int nr = 0;
0090       //int ng = 0;
0091
0092       for(int iev=0;iev<ntp_vertex->GetEntries();iev++)
0093     {
0094       if(verbose) cout << " iev = " << iev << " nr " << nr << endl;
0095
0096       int recoget = ntp_vertex->GetEntry(iev);
0097       if(!recoget)
0098         {
0099           cout << "Failed to get ntp_vertex entry " << iev << endl;
0100           exit(1);
0101         }
0102       if(iev%1 == 0)
0103         if(verbose)
0104           cout << "Get event " << iev << " with vertex x " << evx
0105            << " vertex y " << evy << " vertex z " << evz
0106            << " ngtracks " << ngtracks << " ntracks " << ntracks
0107            << endl;
0108
0109       // loop over tracks
0110
0111       for(int ir=nr;ir<nr+ntracks;ir++)
0112         {
0113           int recoget = ntp_track->GetEntry(ir);
0114           if(!recoget)
0115         {
0116           if(verbose) cout << "Failed to get ntp_track entry " << ir << endl;
0117           break;
0118         }
0119           if(verbose > 0) cout << " ir = " << ir << " rtrackid " << rtrackid << " revent " << revent
0120                    << " rquality " << rquality << "  rgtrackid = " << rgtrackid << " rgnhits " << rgnhits
0121                    << " rnhits " << rnhits << " rnfromtruth " << rnfromtruth << endl;
0122
0123           // Make some overall reconstructed track cuts
0124           if(rgnhits < nlayers-nmissed)
0125         {
0126           //skip tracks that do not pass through enough layers in truth, same cut made on truth tracks earlier
0127           cout << "   --- skip because rgnhits too small " << endl;
0128           continue;
0129         }
0130
0131           if(rquality > quality_cut)
0132         {
0133           if(verbose > 0) cout << "   --- failed quality cut - rejected " << endl;
0134           continue;
0135         }
0136
0137           if(fabs(rvz - evz) > 0.1)
0138         {
0139           if(verbose) cout << "skip because failed z vertex cut with rvz = " << rvz << " evz = " << evz << endl;
0140           continue;
0141         }
0142
0143           if( isnan(rdca2d) )
0144         {
0145           if(verbose) cout << "skip because dca2d is nan" << endl;
0146           continue;
0147         }
0148
0149           if(rdca2d > 0.1)
0150         {
0151           if(verbose) cout << "skip because failed dca2d cut " << endl;
0152           continue;
0153         }
0154
0155           //double rgpT = sqrt(rgpx*rgpx+rgpy*rgpy);
0156           //double rpT = sqrt(rpx*rpx+rpy*rpy);
0157
0158           if(rgembed == embed_flag)
0159         {
0160           n_embed_rtrack++;
0161           naccept++;
0162
0163           double reta = asinh(rpz/sqrt(rpx*rpx+rpy*rpy));
0164           double rmom = sqrt(rpx*rpx+rpy*rpy+rpz*rpz);
0165           double EoverP_cemc = cemc_e / rmom;
0166           cout << "        reta " << reta << " rmom " << rmom << " EoverP " << EoverP_cemc << endl;
0167
0168           if(rgflavor ==11)
0169             hEoverP_cemc[0]->Fill(rmom, EoverP_cemc);
0170           else
0171             hEoverP_cemc[1]->Fill(rmom, EoverP_cemc);
0172
0173         }
0174
0175         }  // end loop over reco'd tracks
0176       if(verbose > 0)  cout << " Done with tracks loop: n_embed_rtrack = " << n_embed_rtrack << " naccept = " << naccept << endl;
0177       //====================================================
0178
0179       // set the starting positions in the track ntuples for the next event
0180       nr += ntracks;
0181       //ng += ngtracks;
0182     }  // end loop over events
0183
0184
0185       delete ntp_gtrack;
0186       delete ntp_track;
0187       delete ntp_cluster;
0188       delete ntp_vertex;
0189
0190     } // end loop over files
0191
0192   TCanvas *cep_cemc = new TCanvas("cep_cemc","",5,5,800,500);
0193   cep_cemc->Divide(2,1);
0194
0195   cep_cemc->cd(1);
0196   hEoverP_cemc[0]->Draw();
0197
0198   cep_cemc->cd(2);
0199   hEoverP_cemc[1]->Draw();
0200
0201   TCanvas *cproj_cemc = new TCanvas("cproj_cemc","",20,20,1400,800);
0202   cproj_cemc->Divide(3,2);
0203
0204   double dp = 2.0;
0205   int np = 12.0 / dp;
0206   for(int ip=0;ip<np;ip++)
0207     {
0208       cproj_cemc->cd(ip+1);
0209
0210       double plo = ip * dp;
0211       double phi = ip * dp + dp;
0212       int binlo = hEoverP_cemc[0]->GetXaxis()->FindBin(plo);
0213       int binhi = hEoverP_cemc[0]->GetXaxis()->FindBin(phi);
0214
0215       char label[500];
0216       sprintf(label,"p = %.1f  - %.1f GeV/c",plo, phi);
0217       TH1D *h = new TH1D("h",label, 200, 0 , 2.0);
0218       hEoverP_cemc[0]->ProjectionY("h",binlo,binhi);
0219       h->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0220       h->GetXaxis()->SetTitle("E/p for CEMC");
0221       h->GetXaxis()->SetTitleOffset(1.0);
0222
0223       h->DrawCopy();
0224
0225       delete h;
0226     }
0227
0228
0229 }