analysis/FrawleyTrackEvalMacros/purity.C

0001 #include <TChain.h>
0002 #include <TNtuple.h>
0003 #include <TCanvas.h>
0004 #include <TH1.h>
0005 #include <TH1F.h>
0006 #include <TH2D.h>
0007 #include <TF1.h>
0008 #include <TFile.h>
0009 #include <TGraph.h>
0010 #include <TStyle.h>
0011 #include <TROOT.h>
0012 #include <TLegend.h>
0013 #include <TLine.h>
0014 #include <TLatex.h>
0015 #include <TRandom1.h>
0016 #include <TPolyLine.h>
0017 #include <iostream>
0018 #include <fstream>
0019 #include <TMath.h>
0020 #include <TLorentzVector.h>
0021 #include <TEfficiency.h>
0022 #include <TGraphAsymmErrors.h>
0023
0024 using namespace std;
0025 #define use_events
0026
0027 void purity()
0028 {
0029   gROOT->SetStyle("Plain");
0030   gStyle->SetOptStat(0);
0031   gStyle->SetOptFit(0);
0032   gStyle->SetOptTitle(0);
0033
0034   int verbose = -1;
0035
0036   // Setup parameters
0037   //=================
0038   int embed_flag = 2;  // set to match value used when events were generated
0039
0040   // these should match the values used when the files were created
0041   //int n_maps_layer = 3;
0042   int n_maps_layer = 0;
0043   //int n_intt_layer = 4;
0044   int n_intt_layer = 0;
0045   int n_tpc_layer = 48; // 16 inner
0046   int nlayers = n_maps_layer+n_intt_layer+n_tpc_layer;  // maximum number of tracking layers for tpc+intt+maps
0047
0048   // track cuts
0049   //static const int nmissed = 20;  // maximum number of missed layers to allow for a track
0050   double quality_cut = 3.0;
0051    double dca2d_cut = 0.1;  // cm
0052   double dcaz_cut = 0.1;  // cm
0053   double gnhits_cut = 20;
0054   //double gnhits_cut = 10;
0055   double truth_hits_cut = 0.8;
0056
0057   // other parameters
0058   double ptmax = 12.2;  // used for Hijing tracks only
0059   //double hptmax = 10.0;
0060   double hptmax = 50.0;
0061
0062   //===============================
0063
0064
0065   //=======================================
0066   // define histograms
0067
0068
0069   TH1D *hnhits = new TH1D("hnhits","hnhits",100,0,99);
0070   TH2D *hpt_nfake = new TH2D("hpt_nfake","hpt_nfake",200,0,hptmax, 73, -5, 68.0);
0071   TH2D *hpt_nmissed_maps_layers = new TH2D("hpt_nmissed_maps_layers","hpt_nmissed_maps_layers",200,0,hptmax, 53, -5, 48.0);
0072   TH2D *hpt_nmissed_intt_layers = new TH2D("hpt_nmissed_intt_layers","hpt_nmissed_intt_layers",200,0,hptmax, 53, -5, 48.0);
0073   TH2D *hpt_nmissed_tpc_layers = new TH2D("hpt_nmissed_tpc_layers","hpt_nmissed_tpc_layers",200,0,hptmax, 53, -5, 48.0);
0074   TH2D *hcorr_nfake_nmaps = new TH2D("hcorr_nfake_nmaps","hcorr_nfake_nmaps",40, -1, 5, 40, -1, 4);
0075   TH1D *hzevt = new TH1D("hzevt","hzevt",200,-35.0, 35.0);
0076   TH1D *hzvtx_res = new TH1D("hzvtx_res","hzvtx_res", 1000, -0.1,  0.1);
0077   TH1D *hxvtx_res = new TH1D("hxvtx_res","hxvtx_res", 1000, -0.04,  0.04);
0078   TH1D *hyvtx_res = new TH1D("hyvtx_res","hyvtx_res", 1000, -0.04,  0.04);
0079   TH1D *hdcavtx_res = new TH1D("hdcavtx_res","hdcavtx_res", 1000, -0.04,  0.04);
0080
0081   static const int NPTDCA = 3;
0082
0083   TH1D *hdca2d[NPTDCA];
0084   for(int ipt=0;ipt<NPTDCA;ipt++)
0085     {
0086       char hname[500];
0087       sprintf(hname,"hdca2d%i",ipt);
0088
0089       hdca2d[ipt] = new TH1D(hname,hname,2000,-0.1,0.1);
0090       hdca2d[ipt]->GetXaxis()->SetTitle("DCA (cm)");
0091       hdca2d[ipt]->GetXaxis()->SetTitleSize(0.055);
0092       hdca2d[ipt]->GetXaxis()->SetLabelSize(0.055);
0093
0094       hdca2d[ipt]->GetYaxis()->SetTitleSize(0.06);
0095       hdca2d[ipt]->GetYaxis()->SetLabelSize(0.055);
0096     }
0097
0098   TH1D *hdca2dsigma = new TH1D("hdca2dsigma","hdca2dsigma",2000,-5.0,5.0);
0099   hdca2dsigma->GetXaxis()->SetTitle("dca2d / dca2dsigma");
0100
0101   TH1D *hZdca = new TH1D("hZdca","hZdca",2000,-0.5,0.5);
0102   hZdca->GetXaxis()->SetTitle("Z DCA (cm)");
0103
0104   TH1D *hquality = new TH1D("hquality","hquality",2000,0.0,20.0);
0105   hquality->GetXaxis()->SetTitle("quality");
0106
0107   TH1D *hgeta = new TH1D("hgeta","hgeta",100,-1.2,1.2);
0108   TH1D *hreta = new TH1D("hreta","hreta",100,-1.2,1.2);
0109
0110   static const int NVARBINS = 36;
0111   double xbins[NVARBINS+1] = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,
0112                 1.1, 1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,
0113                 2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,
0114                   4.5,5.0,5.5,6.0,7.0,8.0,10.0};
0115
0116   TH1D *hpt_truth = new TH1D("hpt_truth","hpt_truth",500, 0.0, hptmax);
0117   TH2D *hpt_compare = new TH2D("hpt_compare","hpt_compare",500, 0.0, hptmax, 2000, 0.0, 2.0);
0118   hpt_compare->GetXaxis()->SetTitle("p_{T}");
0119   hpt_compare->GetYaxis()->SetTitle("#Delta p_{T}/p_{T}");
0120
0121   TH2D *hpt_dca2d = new TH2D("hpt_dca2d","hpt_dca2d",500, 0.0, hptmax, 2000, -0.1, 0.1);
0122   hpt_dca2d->GetXaxis()->SetTitle("p_{T}");
0123   hpt_dca2d->GetYaxis()->SetTitle("dca2d");
0124
0125   TH2D *hpt_dcaZ = new TH2D("hpt_dcaZ","hpt_dcaZ",500, 0.0, hptmax, 2000, -0.1, 0.1);
0126   hpt_dcaZ->GetXaxis()->SetTitle("p_{T}");
0127   hpt_dcaZ->GetYaxis()->SetTitle("dca_Z");
0128
0129   TH1D *hpt_hijing_truth = new TH1D("hpt_hijing_truth","hpt_hijing_truth",500, 0.0, 10.0);
0130   TH2D *hpt_hijing_compare = new TH2D("hpt_hijing_compare","hpt_hijing_compare",500, 0.0, 10.0, 2000, 0.0, 2.0);
0131   hpt_hijing_compare->GetXaxis()->SetTitle("p_{T}");
0132   hpt_hijing_compare->GetYaxis()->SetTitle("#Delta p_{T}/p_{T}");
0133   TH2D *hpt_hijing_dca2d = new TH2D("hpt_hijing_dca2d","hpt_hijing_dca2d",500, 0.0, 10.0, 2000, -0.1, 0.1);
0134   TH2D *hpt_hijing_dcaZ = new TH2D("hpt_hijing_dcaZ","hpt_hijing_dcaZ",500, 0.0, 10.0, 2000, -0.1, 0.1);
0135
0136   TH1D *hptreco[NVARBINS];
0137   for(int ipt=0;ipt<NVARBINS;ipt++)
0138     {
0139       char hn[1000];
0140       sprintf(hn,"hptreco%i",ipt);
0141       hptreco[ipt] = new TH1D(hn, hn, NVARBINS, xbins);
0142     }
0143
0144   TH2D *hg4ntrack = new TH2D("hg4ntrack","hg4ntrack",200,0,3000.0, 200, 0., 2000);
0145
0146   // end of histogram definitions
0147   //===================================================
0148
0149   //============================================================
0150   // Loop over events
0151   //   Loop over reco'd tracks
0152   //     Make quality cut
0153   //       Fill dca2d histos
0154   //       Make loose dca2d cut (usually < 1 mm)
0155   //         Fill Z dca histo
0156   //         Drop tracks outside 1 mm in Z dca from evt vertex
0157   //           Add this reco'd track 2D histo of (rpT-pT) vs pT
0158   //============================================================
0159
0160   // The condor job output files
0161   for(int i=0;i<1200;i++)
0162     {
0163       // Open the evaluator output file
0164
0165       TChain* ntp_track = new TChain("ntp_track","reco tracks");
0166       TChain* ntp_gtrack = new TChain("ntp_gtrack","g4 tracks");
0167       TChain* ntp_vertex = new TChain("ntp_vertex","events");
0168       TChain *ntp_cluster = new TChain("ntp_cluster","clusters");
0169       TChain *ntp_g4hit = new TChain("ntp_g4hit","g4hits");
0170
0171       // This include file contains the definitions of the ntuple variables, and the chain definitions
0172 #include "ntuple_variables.C"
0173
0174       char name[1000];
0175       // latest files
0176       sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
0177
0178       // for CD1 plots with no MVTX
0179       //sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_noINTT_80ns_100pions_clusterizer_fixed/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
0180       //sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_noINTT_80ns_100pions_clusterizer_fixed/g4svtx_eval_%i.root_g4svtx_eval.root",i);
0181
0182       //sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_ups1s_100pions_80ns/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
0183
0184
0185       bool checkfiles = false;
0186       int mintracks = 10;
0187       if(checkfiles)
0188     {
0189       // some QA on the input files
0190       TFile*f=TFile::Open(name);
0191       if(!f)
0192         continue;
0193       TNtuple*test_ntp_track=(TNtuple*)f->Get("ntp_track");
0194       if(!test_ntp_track)
0195         {
0196           f->Close();
0197           continue;
0198         }
0199       if(test_ntp_track->GetEntries() < mintracks)
0200         continue;
0201
0202       f->Close();
0203     }
0204
0205       cout << "Adding file number " << i << " with name " << name << endl;
0206
0207       ntp_vertex->Add(name);
0208       ntp_track->Add(name);
0209       ntp_gtrack->Add(name);
0210
0211
0212
0213       // Use the number of g4 tracks (with some cuts) as a measure of centrality
0214       // Examination of the output file shows that the 10% most central events have 1700 or more g4 tracks with these cuts
0215       double ng4tr = 0;
0216       for(int ig = 0;ig < ntp_gtrack->GetEntries();ig++)
0217     {
0218       int recoget = ntp_gtrack->GetEntry(ig);
0219       if(tembed == 0 && tprimary == 1 && fabs(teta) < 1.0 && sqrt(tpx*tpx+tpy*tpy) > 0.2)
0220         ng4tr++;
0221     }
0222
0223       double nrtr = 0;
0224       for(int ir = 0;ir < ntp_gtrack->GetEntries();ir++)
0225     {
0226       int recoget = ntp_track->GetEntry(ir);
0227       double reta = asinh(rpz/rpt);
0228       if(rgembed == 0 && rprimary == 1 && fabs(reta) < 1.0 && rpt > 0.2)
0229         nrtr++;
0230     }
0231       hg4ntrack->Fill(ng4tr, nrtr);
0232
0233       //if(ng4tr < 1250)  // keep 0-20%
0234       //if(ng4tr > 700)  // keep 60-100%
0235       if(ng4tr > 10000) // keep 0-100%
0236     {
0237       delete ntp_gtrack;
0238       delete ntp_track;
0239       delete ntp_cluster;
0240       delete ntp_vertex;
0241       continue;
0242     }
0243       if(verbose> -1)
0244     cout
0245       <<  " ntp_vertex entries: " << ntp_vertex->GetEntries()
0246          << " ntp_gtrack entries: " << ntp_gtrack->GetEntries()
0247          << " ntp_track entries: " << ntp_track->GetEntries()
0248          << endl;
0249
0250       //==================================
0251       // Begin event loop
0252       //==================================
0253
0254       // These keep track of the starting position in the track ntuples for each event
0255       int nr = 0;
0256       int ng = 0;
0257       int nev = ntp_vertex->GetEntries();
0258       //int nev = 1;
0259
0260       for(int iev=0;iev<nev;iev++)
0261     {
0262
0263       if(verbose > 0) cout << " iev = " << iev << " ng " << ng << " nr " << nr << endl;
0264
0265       int recoget = ntp_vertex->GetEntry(iev);
0266       if(!recoget)
0267         {
0268           cout << "Failed to get ntp_vertex entry " << iev << endl;
0269           continue;
0270         }
0271
0272       /*
0273       if(egvt != 0)
0274         {
0275           // this is a pileup event, we want only the triggered event, so skip it
0276           // we need to advance the pointers to the track ntuples
0277           if(isnan(ntracks))
0278         ntracks  = 0;
0279           nr += ntracks;
0280           ng += ngtracks;
0281
0282           if(verbose>0)
0283         cout << "   skip pileup event with ngtracks = " << ngtracks << " and ntracks = " << ntracks
0284              << " advance ng to " << ng << " and nr to " << nr
0285              << endl;
0286
0287           continue;
0288         }
0289       */
0290
0291       if(iev%1 == 0)
0292         if(verbose>0)
0293           cout << "Get event " << iev
0294            << " with vertex x " << evx
0295            << " vertex y " << evy << " vertex z " << evz
0296            << " egvz " << egvz
0297            << " ngtracks " << ngtracks << " ntracks " << ntracks
0298            << " gtarcks start at " << ng << " ntracks start at " << nr
0299            << endl;
0300
0301       if(fabs(evz - egvz) > 0.1)
0302         {
0303           cout << "     ALERT! event vertex is not correct: egvz = " << egvz << " evz = " << evz << " evx = " << evx << " evy = " << evy << endl;
0304           continue;
0305         }
0306
0307
0308       // do not consider events that are not within the full acceptance
0309       if(fabs(egvz) > 10.0)
0310         {
0311           cout << " Event vertex is outside 10 cm, reject this event - egvz = " << egvz << endl;
0312           if(isnan(ntracks))
0313         ntracks  = 0;
0314           nr += ntracks;
0315           ng += ngtracks;
0316           continue;
0317         }
0318
0319       hzevt->Fill(evz);
0320       hzvtx_res->Fill(evz-egvz);
0321       hxvtx_res->Fill(evx-egvx);
0322       hyvtx_res->Fill(evy-egvy);
0323       hdcavtx_res->Fill(sqrt(evx*evx+evy*evy) - sqrt(egvx*egvx+egvy*egvy));
0324
0325       //====================================================
0326       // ntp_gtracks
0327       //====================================================
0328       // ngtracks is defined in ntuple_variables.C and is the number of g4 tracks
0329       // ntracks is defined in ntuple_variables.C and is the number of reco'd tracks
0330
0331       if(verbose > 0) cout << "Process truth tracks:" << endl;
0332
0333       int n_embed_gtrack = 0;
0334       //for(int ig=ng;ig<ng+ngtracks;ig++)
0335       for(int ig=ng;ig<ntp_gtrack->GetEntries();ig++)
0336         {
0337           int recoget = ntp_gtrack->GetEntry(ig);
0338           if(recoget == 0)
0339         {
0340           //if(verbose) cout << "Failed to get ntp_gtrack entry " << ig << " in ntp_gtrack" << endl;
0341           cout << "Failed to get ntp_gtrack entry " << ig << " in ntp_gtrack" << endl;
0342           break;
0343         }
0344
0345           //if(tembed != embed_flag)
0346           //    continue;
0347
0348           /*
0349           // if the scan_for_embedded flag is set to true, the number of reco tracks recorded in ntp_vertex will be all reco'd tracks, but only embedded tracks will be present
0350           // check for change of event number to detect this
0351           if(tevent != iev)
0352         {
0353           if(verbose > 0) cout << "Change of event number from " << iev << " to " << tevent <<  " go  to next event" << endl;
0354           ngtracks = ig - ng;
0355           break;
0356         }
0357
0358           if(verbose > 0) cout << " ig = " << ig << " tevent " << tevent << " tgtrackid " << tgtrackid
0359                    << " ttrackid " << ttrackid << " tgnhits " << tgnhits << " tnhits " << tnhits << " tnfromtruth " << tnfromtruth << " tembed " << tembed  << endl;
0360           */
0361
0362           // ntp_gtrack track cuts
0363           //===================
0364           // skip tracks that do not pass through enough layers (judged using truth track value of nhits)
0365           if(tgnhits <= gnhits_cut)
0366         {
0367           if(verbose>0) cout << "  -------- failed gnhits cut " << endl;
0368           continue;
0369         }
0370
0371           if(tnfromtruth / tnhits <= truth_hits_cut)
0372         {
0373           if(verbose>0) cout << " -------  failed nfromtruth/nhits cut " << endl;
0374           continue;
0375         }
0376           if(tgtrackid < 0)
0377         {
0378           if(verbose>0) cout << "  ------- failed tgtrackid cut " << endl;
0379           continue;
0380         }
0381
0382           if(tquality > quality_cut)
0383         {
0384           if(verbose > 0) cout << "   --------  failed quality cut - rejected " << endl;
0385           continue;
0386         }
0387
0388           if(fabs(trdca3dz) > dcaz_cut)
0389         {
0390           if(verbose>0) cout << "   --------  skip because track failed z vertex cut with dca3dz  = " << trdca3dz << endl;
0391           continue;
0392         }
0393
0394           if(trdca2d > dca2d_cut)
0395         {
0396           if(verbose>0) cout << "  -------- skip because failed dca2d cut, rdca2d =  " << rdca2d << endl;
0397           continue;
0398         }
0399           // get the truth pT
0400           double tgpT = sqrt(tpx*tpx+tpy*tpy);
0401
0402           if(tembed == embed_flag)  // take only embedded pions
0403         {
0404           double geta = asinh(tpz/sqrt(tpx*tpx+tpy*tpy));
0405           hgeta->Fill(geta); // optional cut
0406           if(fabs(geta) < 1.0)
0407             {
0408               n_embed_gtrack++;
0409               //cout << "fill htp_truth with tgpT = " << tgpT << endl;
0410               hpt_truth->Fill(tgpT);
0411
0412               // Capture hpt_compare from ntp_gtrack
0413               hpt_compare->Fill(tpt,trpt/tpt);
0414               hpt_dca2d->Fill(tpt, trdca2d);
0415               hpt_dcaZ->Fill(tpt, trdca3dz);
0416             }
0417         }
0418
0419           // record embedded pions and Hijing tracks separately
0420           if(tembed == 0)
0421         {
0422           hpt_hijing_truth->Fill(tgpT);
0423         }
0424
0425         } // end loop over truth tracks
0426       if(verbose>0) cout << "n_embed_gtrack = " << n_embed_gtrack << endl;
0427       //===================================================
0428
0429
0430       //====================================================
0431       // ntp_tracks
0432       //====================================================
0433
0434       if(verbose > 0) cout << "Process reco tracks:" << endl;
0435
0436       int n_embed_rtrack = 0;
0437       int naccept = 0;
0438       //for(int ir=nr;ir<nr+ntracks;ir++)
0439       for(int ir=nr;ir<ntp_track->GetEntries();ir++)
0440         {
0441           int recoget = ntp_track->GetEntry(ir);
0442           if(!recoget)
0443         {
0444           if(verbose>0) cout << "Failed to get ntp_track entry " << ir << endl;
0445           break;
0446         }
0447
0448           if(rgembed != embed_flag)
0449         continue;
0450
0451           /*
0452           // if the scan_for_embedded flag is set to true, the number of reco tracks recorded in ntp_vertex will be all reco'd tracks, but only embedded tracks will be present
0453           // check for change of event number to detect this
0454           if(revent != iev)
0455         {
0456           if(verbose > 0) cout << "Change of event number from " << iev << " to " << revent <<  " go  to next event" << endl;
0457           ntracks = ir - nr;
0458           break;
0459         }
0460
0461           if(verbose > 0) cout << " ir = " << ir << " rtrackid " << rtrackid << " revent " << revent
0462                    << " rquality " << rquality << "  rgtrackid = " << rgtrackid << " rgnhits " << rgnhits
0463                    << " rnhits " << rnhits << " rnfromtruth " << rnfromtruth << " rgembed " << rgembed << endl;
0464           */
0465
0466           //  ntp_track track cuts
0467           //================================
0468           if(rgnhits <= gnhits_cut)
0469         {
0470           //skip tracks that do not pass through enough layers in truth, same cut made on truth tracks earlier
0471           cout << "   -------- skip because rgnhits too small " << endl;
0472           continue;
0473         }
0474
0475           if(rnfromtruth / rnhits <= truth_hits_cut)
0476         {
0477           if(verbose>0) cout << " -------  failed rnfromtruth/rnhits cut " << endl;
0478           continue;
0479         }
0480
0481           if(rgtrackid < 0)
0482         {
0483           if(verbose>0) cout << " -------  failed rgtrackid < 0 cut " << endl;
0484           continue;
0485         }
0486           if(rquality > quality_cut)
0487         {
0488           if(verbose > 0) cout << "   --------  failed quality cut - rejected " << endl;
0489           continue;
0490         }
0491
0492           if(fabs(rpcaz - rvz) > dcaz_cut)
0493         {
0494           if(verbose>0) cout << "   --------  skip because track " << ir <<  " failed z vertex cut with rvz = " << rpcaz << " gvz = " << rvz << endl;
0495           continue;
0496         }
0497
0498           if(rdca2d > dca2d_cut)
0499         {
0500           if(verbose>0) cout << "  -------- skip because failed dca2d cut, rdca2d =  " << rdca2d << endl;
0501           continue;
0502         }
0503
0504           /*
0505           if(rnmaps != n_maps_layer)
0506         {
0507           if(verbose>0) cout << "skip track " << ir << " because nmaps = " << rnmaps << endl;
0508           continue;
0509         }
0510           */
0511
0512           //=============================
0513
0514           double rdcaZ = rpcaz - rvz;
0515
0516           if(rgembed == embed_flag)  // take only embedded pions
0517         {
0518           // get the quality and Zdca histos before vertex cuts
0519           hquality->Fill(rquality);
0520           hZdca->Fill(rpcaz - rvz);
0521
0522           if(verbose > 0) cout << "   accepted: rgembed = " << rgembed << " rgpt = " << rgpt << endl;
0523
0524           naccept++;
0525
0526           if(verbose > 0) cout << " rdca2d = " << rdca2d << " rdcaZ = " << rdcaZ << endl;
0527
0528           double reta = asinh(rpz/rpt);
0529           hreta->Fill(reta);
0530           if(fabs(reta) < 1.0)  // optional cut
0531             {
0532               n_embed_rtrack++;
0533               //hpt_compare->Fill(rgpt,rpt/rgpt);
0534               if(rnmaps == n_maps_layer)  // require all maps layers
0535             {
0536               //hpt_dca2d->Fill(rgpt, rdca2d);
0537               //hpt_dcaZ->Fill(rgpt,  rdcaZ);
0538             }
0539             }
0540           hnhits->Fill(rnhits);
0541
0542           double nfake = rnhits - rnfromtruth;
0543           hpt_nfake->Fill(rgpt, nfake);
0544
0545           double nmissed_maps_layers = rgnlmaps - rnlmaps;
0546           hpt_nmissed_maps_layers->Fill(rgpt,nmissed_maps_layers);
0547
0548           double nmissed_intt_layers = rgnlintt - rnlintt;
0549           hpt_nmissed_intt_layers->Fill(rgpt,nmissed_intt_layers);
0550
0551
0552           double nmissed_tpc_layers = rgnltpc - rnltpc;
0553           hpt_nmissed_tpc_layers->Fill(rgpt,nmissed_tpc_layers);
0554
0555           hcorr_nfake_nmaps->Fill(nfake,nmissed_maps_layers);
0556         }
0557
0558           if(rgembed == 1)
0559         {
0560           // for the non-embedded tracks from Hijing, we want to be able to do pt_sigma cuts later
0561           hpt_hijing_compare->Fill(rgpt,rpt/rgpt);
0562           if(rnmaps == n_maps_layer) // require all maps layers
0563             {
0564               hpt_hijing_dca2d->Fill(rgpt,rdca2d);
0565               hpt_hijing_dcaZ->Fill(rgpt,rdcaZ);
0566             }
0567
0568           for(int ipt=0;ipt<NVARBINS-1;ipt++)
0569             {
0570               if(rpt > xbins[ipt] && rpt < xbins[ipt+1])
0571             hptreco[ipt]->Fill(rpt);
0572             }
0573
0574           // Add to the 3 panel  dca2d histos
0575           if(rgpt > 0.5 && rgpt <= 1.0)  hdca2d[0]->Fill(rdca2d);
0576           if(rgpt > 1.0 && rgpt <= 2.0)   hdca2d[1]->Fill(rdca2d);
0577           if(rgpt > 2.0)  hdca2d[2]->Fill(rdca2d);
0578         }
0579
0580         }  // end loop over reco'd tracks
0581       if(verbose > 0)  cout << " Done with loop: n_embed_gtrack = " << n_embed_gtrack << " n_embed_rtrack = " << n_embed_rtrack << " naccept = " << naccept << endl;
0582       //====================================================
0583
0584       // set the starting positions in the track ntuples for the next event
0585       nr += ntracks;
0586       ng += ngtracks;
0587       if(verbose > 0) cout << "  embedded tracks this event = " <<  n_embed_gtrack << " accepted tracks this event = " << naccept << endl;
0588     }  // end loop over events
0589
0590
0591       delete ntp_gtrack;
0592       delete ntp_track;
0593       delete ntp_cluster;
0594       delete ntp_vertex;
0595
0596     } // end loop over files
0597
0598   //==============
0599   // Output the results
0600
0601   TFile *fout = new TFile("root_files/purity_out.root","recreate");
0602
0603   hnhits->Write();
0604
0605   // Three dca2d histos made from Hijing tracks
0606   hdca2d[0]->Write();
0607   hdca2d[1]->Write();
0608   hdca2d[2]->Write();
0609
0610   // truth pT distributions for embedded and non-embedded particles
0611   hpt_truth->Write();
0612   hpt_hijing_truth->Write();
0613
0614   // 2D histos made with embedded pions
0615   hpt_compare->Write();
0616   hpt_dca2d->Write();
0617   hpt_dcaZ->Write();
0618   hpt_nfake->Write();
0619   hpt_nmissed_maps_layers->Write();
0620   hpt_nmissed_intt_layers->Write();
0621   hpt_nmissed_tpc_layers->Write();
0622   hcorr_nfake_nmaps->Write();
0623
0624   // 2d histo made with hijing tracks only
0625   hpt_hijing_compare->Write();
0626   hpt_hijing_dca2d->Write();
0627   hpt_hijing_dcaZ->Write();
0628
0629   // efficiency plot made with Hijing tracks
0630   for(int ipt=0;ipt<NVARBINS;ipt++)
0631     {
0632       hptreco[ipt]->Write();
0633     }
0634
0635   hquality->Write();
0636   hZdca->Write();
0637
0638   hzevt->Write();
0639   hzvtx_res->Write();
0640   hxvtx_res->Write();
0641   hyvtx_res->Write();
0642   hdcavtx_res->Write();
0643
0644   hreta->Write();
0645   hgeta->Write();
0646   hg4ntrack->Write();
0647
0648   fout->Close();
0649
0650
0651 }
0652