analysis/FrawleyTrackEvalMacros/draw_acts_eval.C

0001 #include "TFile.h"
0002 #include "TStyle.h"
0003 #include "TTree.h"
0004 #include "TTreeReader.h"
0005 #include "TTreeReaderValue.h"
0006 #include "TTreeReaderArray.h"
0007 #include "TH1D.h"
0008 #include "TH2D.h"
0009 #include "TCanvas.h"
0010 #include "TChain.h"
0011 #include "TLegend.h"
0012
0013
0014 void draw_acts_eval()
0015 {
0016   gStyle->SetOptStat(0);
0017
0018
0019   unsigned int minhits = 20;
0020   bool verbose1 = false;
0021   bool verbose2 = false;
0022   double ptmax = 40.0;
0023   double maxtracks = 200;
0024
0025   TChain *theTree = new TChain("tracktree");
0026   //theTree->Add("/sphenix/user/frawley/acts_qa/macros/detectors/sPHENIX/G4sPHENIX_g4svtx_eval.root_acts.root");
0027   theTree->Add("/sphenix/user/frawley/new_oct23/macros/detectors/sPHENIX/G4sPHENIX_g4svtx_eval.root_acts.root");
0028
0029   /*
0030   for(int i=0;i<836;++i)
0031     {
0032       char name[500];
0033       sprintf(name,"/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root_acts.root",i);
0034       theTree->Add(name);
0035       cout << "Added file " << name << endl;
0036     }
0037   */
0038   TH1D *hitr[2];
0039   hitr[0] = new TH1D("hitr0","track number",100,1,maxtracks+1);
0040   hitr[1] = new TH1D("hitr1","track number",100,1,maxtracks+1);
0041
0042   TH1D *hpt = new TH1D("hpt","truth_pt",100,0,ptmax);
0043   TH1D *hpt_fit_seed = new TH1D("hpt_fit_seed","fit p_{T} - seed p_{T}",100,-10.0, 10.0);
0044   TH1D *hpt_flt_over_truth = new TH1D("hpt_flt_over_truth","flt p_{T} / truth pT",200,0, 2.0);
0045   TH1D *hpt_actsfit_over_truth = new TH1D("hpt_actsfit_over_truth","actsfit p_{T} / truth pT",200,0, 2.0);
0046   TH2D *hpt_seed_truth = new TH2D("hpt_seed_truth","seed/fit p_{T} vs truth p_{T}",100, 0, ptmax, 100, 0, ptmax);
0047   TH2D *hpt_flt_truth = new TH2D("hpt_flt_truth","flt p_{T}/truth p_{T} vs truth p_{T}",100,0, ptmax, 500, 0.5, 1.5);
0048   TH2D *hpt_actsfit_truth = new TH2D("hpt_actsfit_truth","acts fit p_{T}/truth p_{T} vs truth p_{T}",100,0, ptmax, 500, 0.5, 1.5);
0049
0050   TH2D *hxdiff_KF_truth_xlocal = new TH2D("hxdiff_KF_truth_xlocal","local KF x projection vs truth x (last layer, mm)",100, 0, ptmax, 100,-2.0,2.0);
0051   TH2D *hxdiff_KF_truth_ylocal = new TH2D("hxdiff_KF_truth_ylocal","local KF y projection vs truth y (last layer, mm)",100, 0, ptmax, 100,-2.0,2.0);
0052   TH2D *hxdiff_KF_truth_xglobal = new TH2D("hxdiff_KF_truth_xglobal","global KF x projection vs truth x (last layer, mm)",100, 0, ptmax, 100,-2.0,2.0);
0053   TH2D *hxdiff_KF_truth_yglobal = new TH2D("hxdiff_KF_truth_yglobal","global KF y projection vs truth y (last layer, mm)",100, 0, ptmax, 100,-2.0,2.0);
0054   TH2D *hxdiff_KF_truth_zglobal = new TH2D("hxdiff_KF_truth_zglobal","global KF z projection vs truth z (last_layer, mm)",100, 0, ptmax, 100,-2.0,2.0);
0055
0056   TH1D *hxerr_KF = new TH1D("hxerr_KF"," err_eLOC0_flt", 100, -3, 3);
0057   TH1D *hxpull_KF = new TH1D("hxpull_KF"," pull_eLOC0_flt outermost", 100, -5, 5);
0058   TH1D *hypull_KF = new TH1D("hypull_KF"," pull_eLOC1_flt outermost", 100, -5, 5);
0059   TH2D *hxpull_KF_radius = new TH2D("hxpull_KF_radius"," pull_eLOC0_flt vs radius", 860, 0, 860, 100, -5, 5);
0060   TH2D *hzpull_KF_radius = new TH2D("hzpull_KF_radius"," pull_eLOC1_flt vs radius", 860, 0, 860, 100, -5, 5);
0061
0062   TH1D *hpt_smt = new TH1D("hpt_smt","state/truth p_{T}",100,0, ptmax);
0063   TH2D *htxy = new TH2D("htxy","y vs x truth",4000,-860,860,4000,-860,860);
0064   TH2D *hxy_flt = new TH2D("hxy_flt","hxy_flt",4000,-860,860,4000,-860,860);
0065   TH2D *hxy_smt = new TH2D("hxy_smt","hxy_smt",4000,-860,860,4000,-860,860);
0066   TH2D *hxy_prt = new TH2D("hxy_prt","hxy_prt",4000,-860,860,4000,-860,860);
0067
0068   TH1D *hdcaxy[2];
0069   TH1D *hdcaz[2];
0070   hdcaxy[0] = new TH1D("hdcaxy0","DCA xy", 500, -2.0, 2);
0071   hdcaxy[0]->GetXaxis()->SetTitle("DCA (mm)");
0072   hdcaxy[1] = new TH1D("hdcaxy1","DCA xy", 500, -2.0, 2);
0073   hdcaz[0] = new TH1D("hdcaz0","DCA z", 500, -2.0, 2.0);
0074   hdcaz[1] = new TH1D("hdcaz1","DCA z", 500, -2.0, 2.0);
0075
0076
0077   TH1D *hprotox[2];
0078   TH1D *hprotoy[2];
0079   TH1D *hprotoz[2];
0080   double range = 3.0;
0081   hprotox[0] = new TH1D("hprotox0","proto x - vertex x", 100, -range, range);
0082   hprotox[1] = new TH1D("hprotox1","proto x - vertex x", 100, -range, range);
0083   hprotoy[0] = new TH1D("hprotoy0","proto y - vertex y", 100, -range, range);
0084   hprotoy[1] = new TH1D("hprotoy1","proto y - vertex y", 100, -range, range);
0085   hprotoz[0] = new TH1D("hprotoz0","proto z - vertex z", 100, -range, range);
0086   hprotoz[1] = new TH1D("hprotoz1","proto z - vertex z", 100, -range, range);
0087
0088   TH1D *hhit[5];
0089   hhit[0] = new TH1D("hhitx_global","global x hit - truth hit (last layer, mm)", 200, -0.5,0.5);
0090   hhit[1] = new TH1D("hhity_global","global y hit - truth (last layer, mm)", 200, -0.5,0.5);
0091   hhit[2] = new TH1D("hhitz_global","global z hit - truth (last layer, mm)", 200, -3,3);
0092   hhit[3] = new TH1D("hhitx_local","local x hit - truth (last layer, mm)", 200, -0.6,0.6);
0093   hhit[4] = new TH1D("hhity_local","local y hit - truth (last_layer, mm)", 200, -3,3);
0094
0095
0096   //TTree* theTree = nullptr;
0097   //fin->GetObject("tracktree",theTree);
0098   //theTree->Print();
0099   TTreeReader theReader(theTree);
0100
0101   TTreeReaderValue<int> event(theReader, "event_nr");
0102
0103   // Truth
0104   TTreeReaderArray<float> t_x(theReader, "t_x");
0105   TTreeReaderArray<float> t_y(theReader, "t_y");
0106   TTreeReaderArray<float> t_z(theReader, "t_z");
0107   TTreeReaderValue<float> t_pT(theReader, "t_pT");
0108   TTreeReaderValue<float> t_vx(theReader, "t_vx");
0109   TTreeReaderValue<float> t_vy(theReader, "t_vy");
0110   TTreeReaderValue<float> t_vz(theReader, "t_vz");
0111   TTreeReaderValue<unsigned long long> t_barcode(theReader, "t_barcode");
0112   TTreeReaderArray<float> t_eLOC0(theReader, "t_eLOC0");
0113   TTreeReaderArray<float> t_eLOC1(theReader, "t_eLOC1");
0114   //TTreeReaderArray<float> layer(theReader, "layer_id");
0115
0116   // proto track
0117   TTreeReaderValue<float> x_proto(theReader, "g_protoTrackX");
0118   TTreeReaderValue<float> y_proto(theReader, "g_protoTrackY");
0119   TTreeReaderValue<float> z_proto(theReader, "g_protoTrackZ");
0120   TTreeReaderValue<float> px_proto(theReader, "g_protoTrackPx");
0121   TTreeReaderValue<float> py_proto(theReader, "g_protoTrackPy");
0122   TTreeReaderValue<float> pz_proto(theReader, "g_protoTrackPz");
0123
0124   // Predicted states
0125   TTreeReaderArray<float> pT_prt(theReader, "pT_prt");
0126   TTreeReaderArray<float> g_x_prt(theReader, "g_x_prt");
0127   TTreeReaderArray<float> g_y_prt(theReader, "g_y_prt");
0128   TTreeReaderArray<float> g_z_prt(theReader, "g_z_prt");
0129   TTreeReaderArray<float> eLOC0_prt(theReader, "eLOC0_prt");
0130   TTreeReaderArray<float> eLOC1_prt(theReader, "eLOC1_prt");
0131   TTreeReaderArray<float> pull_eLOC0_prt(theReader, "pull_eLOC0_prt");
0132   TTreeReaderArray<float> pull_eLOC1_prt(theReader, "pull_eLOC1_prt");
0133   TTreeReaderArray<float> err_eLOC0_prt(theReader, "err_eLOC0_prt");
0134   TTreeReaderArray<float> err_eLOC1_prt(theReader, "err_eLOC1_prt");
0135   TTreeReaderArray<float> res_eLOC0_prt(theReader, "res_eLOC0_prt");
0136   TTreeReaderArray<float> res_eLOC1_prt(theReader, "res_eLOC1_prt");
0137   TTreeReaderArray<float> px_prt(theReader, "px_prt");
0138   TTreeReaderArray<float> py_prt(theReader, "py_prt");
0139   TTreeReaderArray<float> pz_prt(theReader, "pz_prt");
0140
0141   // Filtered states
0142   TTreeReaderArray<float> g_x_flt(theReader, "g_x_flt");
0143   TTreeReaderArray<float> g_y_flt(theReader, "g_y_flt");
0144   TTreeReaderArray<float> g_z_flt(theReader, "g_z_flt");
0145   TTreeReaderArray<float> eLOC0_flt(theReader, "eLOC0_flt");
0146   TTreeReaderArray<float> eLOC1_flt(theReader, "eLOC1_flt");
0147   TTreeReaderArray<float> pull_eLOC0_flt(theReader, "pull_eLOC0_flt");
0148   TTreeReaderArray<float> pull_eLOC1_flt(theReader, "pull_eLOC1_flt");
0149   TTreeReaderArray<float> err_eLOC0_flt(theReader, "err_eLOC0_flt");
0150   TTreeReaderArray<float> err_eLOC1_flt(theReader, "err_eLOC1_flt");
0151   TTreeReaderArray<float> res_eLOC0_flt(theReader, "res_eLOC0_flt");
0152   TTreeReaderArray<float> res_eLOC1_flt(theReader, "res_eLOC1_flt");
0153   TTreeReaderArray<float> pT_flt(theReader, "pT_flt");
0154   TTreeReaderArray<float> px_flt(theReader, "px_flt");
0155   TTreeReaderArray<float> py_flt(theReader, "py_flt");
0156   TTreeReaderArray<float> pz_flt(theReader, "pz_flt");
0157
0158   // Smoothed states
0159   TTreeReaderArray<float> g_x_smt(theReader, "g_x_smt");
0160   TTreeReaderArray<float> g_y_smt(theReader, "g_y_smt");
0161   TTreeReaderArray<float> g_z_smt(theReader, "g_z_smt");
0162   TTreeReaderArray<float> eLOC0_smt(theReader, "eLOC0_smt");
0163   TTreeReaderArray<float> eLOC1_smt(theReader, "eLOC1_smt");
0164   TTreeReaderArray<float> pull_eLOC0_smt(theReader, "pull_eLOC0_smt");
0165   TTreeReaderArray<float> pull_eLOC1_smt(theReader, "pull_eLOC1_smt");
0166   TTreeReaderArray<float> err_eLOC0_smt(theReader, "err_eLOC0_smt");
0167   TTreeReaderArray<float> err_eLOC1_smt(theReader, "err_eLOC1_smt");
0168   TTreeReaderArray<float> res_eLOC0_smt(theReader, "res_eLOC0_smt");
0169   TTreeReaderArray<float> res_eLOC1_smt(theReader, "res_eLOC1_smt");
0170   TTreeReaderArray<float> pT_smt(theReader, "pT_smt");
0171   TTreeReaderArray<float> px_smt(theReader, "px_smt");
0172   TTreeReaderArray<float> py_smt(theReader, "py_smt");
0173   TTreeReaderArray<float> pz_smt(theReader, "pz_smt");
0174
0175   // Fitted
0176   TTreeReaderValue<float> dcaxy_fit(theReader, "g_dca3Dxy_fit");
0177   TTreeReaderValue<float> dcaz_fit(theReader, "g_dca3Dz_fit");
0178   TTreeReaderValue<float> x_fit(theReader, "g_x_fit");
0179   TTreeReaderValue<float> y_fit(theReader, "g_y_fit");
0180   TTreeReaderValue<float> z_fit(theReader, "g_z_fit");
0181   TTreeReaderValue<float> px_fit(theReader, "g_px_fit");
0182   TTreeReaderValue<float> py_fit(theReader, "g_py_fit");
0183   TTreeReaderValue<float> pz_fit(theReader, "g_pz_fit");
0184
0185
0186   // Actual hits
0187   TTreeReaderArray<float> g_x_hit(theReader, "g_x_hit");
0188   TTreeReaderArray<float> g_y_hit(theReader, "g_y_hit");
0189   TTreeReaderArray<float> g_z_hit(theReader, "g_z_hit");
0190   TTreeReaderArray<float> l_x_hit(theReader, "l_x_hit");
0191   TTreeReaderArray<float> l_y_hit(theReader, "l_y_hit");
0192   TTreeReaderArray<float> err_x_hit(theReader, "err_x_hit");
0193   TTreeReaderArray<float> err_y_hit(theReader, "err_y_hit");
0194   TTreeReaderArray<float> res_x_hit(theReader, "res_x_hit");
0195   TTreeReaderArray<float> res_y_hit(theReader, "res_y_hit");
0196   TTreeReaderArray<float> pull_x_hit(theReader, "pull_x_hit");
0197   TTreeReaderArray<float> pull_y_hit(theReader, "pull_y_hit");
0198
0199   int itr = 0;
0200    while(theReader.Next()){
0201
0202      //if( !( (*t_barcode == 1) || (*t_barcode == 23) ) ) continue;
0203
0204      if(pT_flt.GetSize() < 1)
0205        {
0206      cout << " pT_smt not found " << endl;
0207      continue;
0208        }
0209
0210      if(pT_flt.GetSize() < minhits) continue;
0211
0212      double inner_radius = sqrt(pow(t_x[pT_smt.GetSize()-1], 2) + pow(t_y[pT_smt.GetSize()-1], 2));
0213      //if(inner_radius > 80)   continue;
0214
0215      if(pT_smt.GetSize() != pT_flt.GetSize())
0216        cout << " ***********   smt size " << pT_smt.GetSize() << " flt size " << pT_flt.GetSize() << endl;
0217
0218      // cout << " t_barcode " << *t_barcode << endl;
0219      //if(*t_barcode > maxtracks)
0220      //{
0221      //cout << "skip track with barcode " << *t_barcode << endl;
0222      //continue;
0223      //}
0224
0225      float pT_fit =sqrt(pow(*px_fit,2) + pow(*py_fit,2));
0226      float pT_proto =sqrt(pow(*px_proto,2) + pow(*py_proto,2));
0227
0228
0229      double xdiff_global_last = g_x_flt[0] - t_x[0];
0230      double ydiff_global_last = g_y_flt[0] - t_y[0];
0231      double zdiff_global_last = g_z_flt[0] - t_z[0];
0232      double xdiff_local_last = eLOC0_flt[0] - t_eLOC0[0];
0233      double ydiff_local_last = eLOC1_flt[0] - t_eLOC1[0];
0234
0235
0236      if(verbose1)
0237        {
0238      cout << " new track: " << itr << " truth Z vertex " <<  *t_vz  << " inner radius " << inner_radius <<  " nhits " << pT_flt.GetSize() << endl;
0239      //cout << "    Truth pT " << *t_pT << " seed pT " << pT_prt[pT_flt.GetSize() - 1] <<  " KF pT " << pT_flt[0] << " SM pT " << pT_smt[pT_smt.GetSize()-1]
0240      cout << "       track " << itr <<  " truth barcode " << *t_barcode<< " Truth pT " << *t_pT << " seed pT " << pT_prt[pT_flt.GetSize() - 1] <<  " fit pT " << pT_fit << " inner radius " << inner_radius << " nhits " << pT_smt.GetSize();
0241      if(pT_fit / *t_pT < 0.8 || pT_fit/ *t_pT > 1.2)
0242        cout << "    ---- bad pT " << endl;
0243      else
0244        cout << "    ---- good pT " << endl;
0245
0246      cout << "           px_proto " << *px_proto << " py_proto " << *py_proto << " pz_proto " << *pz_proto << " pT_proto " << pT_proto << endl;
0247      //cout << "           px_prt " << px_prt[pT_prt.GetSize() - 1] << "  py_prt " << py_prt[pT_prt.GetSize() - 1]  << " pz_prt " << pz_prt[pT_prt.GetSize() - 1]
0248      //      << " pT_prt " << pT_prt[pT_prt.GetSize() - 1] << endl;
0249      //cout << "           px_flt " << px_flt[0] << "  py_flt " << py_flt[0]  << " pz_flt " << pz_flt[0] << " pT_flt " << pT_flt[0] << endl;
0250      //cout << "           px_smt " << px_smt[0] << "  py_smt " << py_smt[0]  << " pz_smt " << pz_smt[0] << " pT_smt " << pT_smt[0] << endl;
0251      //cout << "           px_smt " << px_smt[pT_smt.GetSize() - 1] << "  py_smt " << py_smt[pT_smt.GetSize() - 1]  << " pz_smt " << pz_smt[pT_smt.GetSize() - 1] << " pT_smt " << pT_smt[pT_smt.GetSize() - 1] << endl;
0252      cout << "           px_fit " << *px_fit << "  py_fit " << *py_fit  << " pz_fit " << *pz_fit << " pT_fit " << pT_fit << endl;
0253      cout << "           x_proto " << *x_proto << " y_proto " << *y_proto << " z_proto " << *z_proto << endl;
0254      cout << "           x_fit   " << *x_fit << " y_fit   " << *y_fit << " z_fit   " << *z_fit << endl;
0255        }
0256
0257      double pT_truth = *t_pT;
0258      double pT_fitted = pT_flt[0];
0259      double pT_actsfit = pT_fit;  // has problems
0260      double pT_seed = pT_prt[pT_prt.GetSize() - 1];
0261
0262      hpt_seed_truth->Fill(pT_truth, pT_fitted);
0263      hpt_flt_truth->Fill(pT_truth, pT_fitted / pT_truth);
0264      hpt_actsfit_truth->Fill(pT_truth, pT_actsfit / pT_truth);
0265      hpt_fit_seed->Fill(pT_fitted -  pT_seed);
0266      hpt_flt_over_truth->Fill(pT_fitted /  pT_truth);
0267      hpt_actsfit_over_truth->Fill(pT_actsfit /  pT_truth);
0268      hxdiff_KF_truth_xlocal->Fill(pT_truth, xdiff_local_last);
0269      hxdiff_KF_truth_ylocal->Fill(pT_truth, ydiff_local_last);
0270      hxdiff_KF_truth_xglobal->Fill(pT_truth, xdiff_global_last);
0271      hxdiff_KF_truth_yglobal->Fill(pT_truth, ydiff_global_last);
0272      hxdiff_KF_truth_zglobal->Fill(pT_truth, zdiff_global_last);
0273      hhit[0]->Fill(g_x_hit[0] - t_x[0]);
0274      hhit[1]->Fill(g_y_hit[0] - t_y[0]);
0275      hhit[2]->Fill(g_z_hit[0] - t_z[0]);
0276      hhit[3]->Fill(l_x_hit[0] - t_eLOC0[0]);
0277      hhit[4]->Fill(l_y_hit[0] - t_eLOC1[0]);
0278
0279      hxerr_KF->Fill( err_eLOC0_flt[0]);
0280      hxpull_KF->Fill( pull_eLOC0_flt[0]);
0281      hypull_KF->Fill( pull_eLOC1_flt[0]);
0282     if(inner_radius < 80)
0283       {
0284     hitr[0]->Fill(*t_barcode);
0285     hdcaxy[0]->Fill(*dcaxy_fit);
0286     //hdcaz[0]->Fill(*dcaz_fit-*t_vz);
0287     hdcaz[0]->Fill(*dcaz_fit);
0288
0289     hprotox[0]->Fill(*x_proto - *t_vx);
0290     hprotoy[0]->Fill(*y_proto - *t_vy);
0291     hprotoz[0]->Fill(*z_proto - *t_vz);
0292       }
0293     else
0294       {
0295     hitr[1]->Fill(*t_barcode);
0296     hdcaxy[1]->Fill(*dcaxy_fit);
0297     //hdcaz[1]->Fill(*dcaz_fit-*t_vz);
0298     hdcaz[1]->Fill(*dcaz_fit);
0299
0300     hprotox[1]->Fill(*x_proto - *t_vx);
0301     hprotoy[1]->Fill(*y_proto - *t_vy);
0302     hprotoz[1]->Fill(*z_proto - *t_vz);
0303
0304       }
0305
0306      for(unsigned int i = 0; i < pT_smt.GetSize(); ++i)
0307        {
0308      htxy->Fill(t_x[i], t_y[i]);
0309      hxy_flt->Fill(g_x_flt[i], g_y_flt[i]);
0310      hxy_smt->Fill(g_x_smt[i], g_y_smt[i]);
0311      hxy_prt->Fill(g_x_prt[i], g_y_prt[i]);
0312
0313      hpt_smt->Fill(pT_smt[i]);
0314      hpt->Fill(*t_pT);
0315
0316      double radius = sqrt(pow(t_x[i], 2) + pow(t_y[i], 2));
0317      hxpull_KF_radius->Fill(radius, pull_eLOC0_flt[i]);
0318      hzpull_KF_radius->Fill( radius, pull_eLOC1_flt[i]);
0319      //if(i == pT_smt.GetSize()-1)
0320        if(verbose2)
0321        {
0322          // Identify the hit
0323          cout << " hit " << i << " rad " << radius
0324           << " flt rphi pull " << pull_eLOC0_flt[i]
0325           << " smt rphi pull " << pull_eLOC0_smt[i]
0326           << endl;
0327          cout << "    pT_truth " << pT_truth << " pT_prt " << pT_prt[i] << " pT_flt " << pT_flt[i] << " pT_smt " << pT_smt[i] << endl;
0328
0329          // global hit positions
0330          cout << "    truth global hit : "  << "t_x     " << t_x[i] << " t_y     " << t_y[i] << " t_z      " << t_z[i]  << endl;
0331          cout << "    Global hit:        g_x_hit " << g_x_hit[i] << " g_y_hit " << g_y_hit[i]  << " g_z_hit " << g_z_hit[i]
0332           << "    dx " << g_x_hit[i] - t_x[i]
0333           << " dy " << g_y_hit[i] - t_y[i]
0334           << " dz " << g_z_hit[i] - t_z[i]
0335           << endl;
0336          /*
0337          cout << "    Global predicted:  g_x_prt " << g_x_prt[i] << " g_y_prt " << g_y_prt[i] << " g_z_prt " << g_z_prt[i]
0338           << "     dx " << g_x_prt[i] - t_x[i]
0339           << " dy " << g_y_prt[i] - t_y[i]
0340           << " dz " << g_z_prt[i] - t_z[i]
0341           << endl;
0342          cout << "    Global filter:     g_x_flt " << g_x_flt[i] << " g_y_flt " << g_y_flt[i]  << " g_z_flt " << g_z_flt[i]
0343           << "     dx " << g_x_flt[i] - t_x[i]
0344           << " dy " << g_y_flt[i] - t_y[i]
0345           << " dy " << g_z_flt[i] - t_z[i]
0346           << endl;
0347          */
0348          cout << "    Global smoothed:   g_x_smt " << g_x_smt[i] << " g_y_smt " << g_y_smt[i]  << " g_z_smt " << g_z_smt[i]
0349           << "     dx " << g_x_smt[i] - t_x[i]
0350           << " dy " << g_y_smt[i] - t_y[i]
0351           << " dy " << g_z_smt[i] - t_z[i]
0352           << endl;
0353
0354          /*
0355          // local hit positions
0356          cout << "    X:      local truth :       t_eLOC0 " << t_eLOC0[i] << endl;
0357
0358          cout << "            Local hit:          l_x_hit " << l_x_hit[i] << "       err_x_hit " << err_x_hit[i] << "    res_x_hit " << res_x_hit[i] << " pull_x_hit " << pull_x_hit[i] << endl;
0359          cout << "            Local predicted:  elOC0_prt " << eLOC0_prt[i] << " err_eLOC0_prt " << err_eLOC0_prt[i] << " res_eLOC0_prt " << res_eLOC0_prt[i] << " pull_eLOC0_prt " << pull_eLOC0_prt[i] << endl;
0360          cout << "            Local filter:     elOC0_flt " << eLOC0_flt[i] << " err_eLOC0_flt " << err_eLOC0_flt[i] << " res_eLOC0_flt " << res_eLOC0_flt[i] << " pull_eLOC0_flt " << pull_eLOC0_flt[i] << endl;
0361          cout << "            Local smoothed:   eLOC0_smt " << eLOC0_smt[i] << " err_eLOC0_smt " << err_eLOC0_smt[i] << " res_eLOC0_smt " << res_eLOC0_smt[i] << " pull_eLOC0_smt " << pull_eLOC0_smt[i] << endl;
0362
0363          cout << "    Y:     local truth :        t_eLOC1 " << t_eLOC1[i] << endl;
0364          cout << "            Local hit:          l_y_hit " << l_y_hit[i] <<  " err_y_hit " << err_y_hit[i] << " res_y_hit " << res_y_hit[i] << " pull_y_hit " << pull_y_hit[i] << endl;
0365          cout << "            Local predicted:  eLOC1_prt " << eLOC1_prt[i] << " err_eLOC1_prt " << err_eLOC1_prt[i] << " res_eLOC1_prt " << res_eLOC1_prt[i] << " pull_eLOC1_prt " << pull_eLOC1_prt[i] << endl;
0366          cout << "            Local filter:     eLOC1_flt " << eLOC1_flt[i] << " err_eLOC1_flt " << err_eLOC1_flt[i] << " res_eLOC1_flt " << res_eLOC1_flt[i] << " pull_eLOC1_flt " << pull_eLOC1_flt[i] << endl;
0367          cout << "            Local smoothed:   eLOC1_smt " << eLOC1_smt[i] << " err_eLOC1_smt " << err_eLOC1_smt[i] << " res_eLOC1_smt " << res_eLOC1_smt[i] << " pull_eLOC1_smt " << pull_eLOC1_smt[i] << endl;
0368          */
0369          /*
0370          // Smoothed track projuections
0371          cout << "    Global smoothed: g_x_smt " << g_x_smt[i] << " g_y_smt " << g_y_smt[i] << " g_z_smt " << g_z_smt[i]
0372          << " dx " << g_x_smt[i] - t_x[i]
0373          << " dy " << g_y_smt[i] - t_y[i]
0374          << " dz " << g_z_smt[i] - t_z[i]
0375          << endl;
0376          cout << "    Local smoothed: elOC0_smt " << eLOC0_smt[i] << " eLOC1_smt " << eLOC1_smt[i] << endl;
0377          cout << "                  Local smoothed: err_eLOC0_smt " << err_eLOC0_smt[i] << " res_eLOC0_smt " << res_eLOC0_smt[i] << " pull_eLOC0_smt " << pull_eLOC0_smt[i] << endl;
0378          cout << "                  Local smoothed: err_eLOC1_smt " << err_eLOC1_smt[i] << " res_eLOC1_smt " << res_eLOC1_smt[i] << " pull_eLOC1_smt " << pull_eLOC1_smt[i] << endl;
0379          */
0380        }
0381        }
0382
0383      // track info
0384      if(verbose2)
0385        {
0386      cout << " track " << itr
0387           << " truth pT = " << *t_pT
0388           << endl;
0389
0390      cout << " dcaxy_fit " << *dcaxy_fit << " dcaz_fit " << *dcaz_fit << endl;
0391        }
0392
0393
0394      itr++;
0395    }
0396
0397    // track vertex position
0398
0399    TCanvas *c = new TCanvas("c","c",5,5,800,800);
0400
0401   TH2D *hvxy = new TH2D("hvxy","y vs x track vertex",4000,-860,860,4000,-860,860);
0402   hvxy->Fill(*t_vx, *t_vy);
0403
0404   /*
0405    // now get the truth clusters and plot those too
0406   TChain* ntp_g4cluster = new TChain("ntp_g4cluster","truth clusters");
0407   ntp_g4cluster->Add("/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output/g4svtx_eval_0.root_g4svtx_eval.root");
0408
0409   TH2D *hclusxy = new TH2D("hclusxy","y vs x (magenta=all-hits, black=acts-hits,  red=prt, green=flt)",4000,-800,800,4000,-800,800);
0410   ntp_g4cluster->Draw("(10.0*gy):(10.0)*gx>>hclusxy");
0411
0412    hpt_smt->Draw();
0413    hpt->SetLineColor(kRed);
0414    hpt->Draw("same");
0415
0416    TCanvas *ch= new TCanvas("ch","ch",5,5,800,800);
0417    ch->SetLeftMargin(0.15);
0418
0419    hclusxy->SetMarkerStyle(20);
0420    hclusxy->SetMarkerSize(0.5);
0421    hclusxy->SetMarkerColor(kMagenta);
0422    hclusxy->GetYaxis()->SetTitle("y (mm)");
0423    hclusxy->GetXaxis()->SetTitle("x (mm)");
0424    hclusxy->Draw();
0425                     */
0426
0427    hvxy->SetMarkerStyle(20);
0428    hvxy->SetMarkerSize(1.0);
0429    hvxy->SetMarkerColor(kCyan);
0430    hvxy->Draw("same");
0431
0432    htxy->SetMarkerStyle(20);
0433    htxy->SetMarkerSize(0.5);
0434    htxy->Draw("same");
0435    hxy_prt->SetMarkerColor(kRed);
0436    hxy_prt->SetMarkerStyle(20);
0437    hxy_prt->SetMarkerSize(0.5);
0438    hxy_prt->Draw("same");
0439    hxy_flt->SetMarkerColor(kGreen+2);
0440    hxy_flt->SetMarkerStyle(20);
0441    hxy_flt->SetMarkerSize(0.5);
0442    hxy_flt->Draw("same");
0443    hxy_smt->SetMarkerColor(kBlue);
0444    hxy_smt->SetMarkerStyle(20);
0445    hxy_smt->SetMarkerSize(0.5);
0446    //hxy_smt->Draw("same");
0447
0448    //
0449    TCanvas *ctruth = new TCanvas("ctruth","ctruth",5,5,1500,800);
0450    ctruth->Divide(3,1);
0451
0452    ctruth->cd(1);
0453    gPad->SetLeftMargin(0.12);
0454    hpt_seed_truth->GetXaxis()->SetTitle("Truth p_{T}");
0455    hpt_seed_truth->GetYaxis()->SetTitle("Seed or fit p_{T}");
0456    hpt_seed_truth->SetMarkerStyle(20);
0457    hpt_seed_truth->SetMarkerSize(0.4);
0458    hpt_seed_truth->GetYaxis()->SetTitleSize(0.05);
0459    hpt_seed_truth->GetXaxis()->SetTitleSize(0.05);
0460    hpt_seed_truth->GetYaxis()->SetLabelSize(0.05);
0461    hpt_seed_truth->GetXaxis()->SetLabelSize(0.05);
0462    hpt_seed_truth->GetXaxis()->SetNdivisions(505);
0463    //hpt_seed_truth->Draw();
0464
0465    TLegend *trleg = new TLegend(0.2, 0.75, 0.55, 0.85, "", "NDC");
0466    trleg->AddEntry(hpt_actsfit_truth,"seed", "p");
0467    trleg->AddEntry(hpt_flt_truth, "fit ","p");
0468    trleg->SetTextSize(0.05);
0469    trleg->Draw();
0470
0471    cout << "get truth plots" << endl;
0472
0473    hpt_flt_truth->SetMarkerStyle(20);
0474    hpt_flt_truth->SetMarkerSize(0.4);
0475    hpt_flt_truth->SetMarkerColor(kRed);
0476    hpt_flt_truth->Draw("colz");
0477
0478    ctruth->cd(2);
0479    gPad->SetLeftMargin(0.12);
0480    hpt_actsfit_over_truth->GetXaxis()->SetTitle("actsfit p_{T}  / truth p_{T}");
0481    hpt_actsfit_over_truth->SetTitleSize(0.05);
0482    hpt_actsfit_over_truth->GetYaxis()->SetLabelSize(0.05);
0483    hpt_actsfit_over_truth->GetXaxis()->SetLabelSize(0.05);
0484    hpt_actsfit_over_truth->GetXaxis()->SetNdivisions(505);
0485    hpt_actsfit_over_truth->Draw();
0486
0487    TLegend *tkleg = new TLegend(0.15, 0.79, 0.48, 0.85, "", "NDC");
0488    tkleg->AddEntry(hpt_actsfit_over_truth, "actsfit/truth","l");
0489    tkleg->SetTextSize(0.05);
0490    tkleg->Draw();
0491
0492    ctruth->cd(3);
0493    gPad->SetLeftMargin(0.12);
0494    hpt_flt_over_truth->SetMarkerStyle(20);
0495    hpt_flt_over_truth->SetMarkerSize(0.2);
0496    hpt_flt_over_truth->SetLineColor(kRed);
0497    hpt_flt_over_truth->GetXaxis()->SetTitle("flt p_{T} / truth p_{T}");
0498    hpt_flt_over_truth->GetYaxis()->SetTitleSize(0.05);
0499    hpt_flt_over_truth->GetXaxis()->SetTitleSize(0.05);
0500    hpt_flt_over_truth->GetYaxis()->SetLabelSize(0.05);
0501    hpt_flt_over_truth->GetXaxis()->SetLabelSize(0.05);
0502    hpt_flt_over_truth->GetXaxis()->SetNdivisions(505);
0503    hpt_flt_over_truth->Draw();
0504
0505    int binlo =  hpt_flt_over_truth->GetXaxis()->FindBin(0.8);
0506    int binhi =  hpt_flt_over_truth->GetXaxis()->FindBin(1.2);
0507    cout << " binlo " << binlo << " binhi " << binhi << endl;
0508    cout << " hpt_flt_over_truth integral 0.8 to 1.2 " << hpt_flt_over_truth->Integral(binlo, binhi) << endl;;
0509    cout << " hpt_flt_over_truth integral 0 to 0.8  " << hpt_flt_over_truth->Integral(1, binlo) << endl;;
0510    cout << " hpt_actsfit_over_truth integral 0.8 to 1.2 " << hpt_actsfit_over_truth->Integral(binlo, binhi) << endl;;
0511
0512
0513    TLegend *tfleg = new TLegend(0.15, 0.79, 0.48, 0.85, "", "NDC");
0514    tfleg->AddEntry(hpt_flt_over_truth, "flt/truth p_{T}","l");
0515    tfleg->SetTextSize(0.05);
0516    tfleg->Draw();
0517
0518    cout << "get pT res plots" << endl;
0519
0520    hpt_flt_truth->FitSlicesY();
0521    TH1D*hptres = (TH1D*)gDirectory->Get("hpt_flt_truth_2");
0522    hptres->GetYaxis()->SetRangeUser(0.0, 0.15);
0523    TH1D*hptcent = (TH1D*)gDirectory->Get("hpt_flt_truth_1");
0524    hptcent->GetYaxis()->SetRangeUser(0.8, 1.2);
0525
0526    cout << "plot pT res plots" << endl;
0527
0528    TCanvas *cptres = new TCanvas("cptres","cptres", 8,8,1200,800);
0529    cptres->Divide(2,1);
0530    cptres->cd(1);
0531    cout << "plot res plot" << endl;
0532    hptres->Draw("p");
0533    cptres->cd(2);
0534    cout << "plot cent plot" << endl;
0535    hptcent->Draw("p");
0536
0537    TCanvas *cdiff = new TCanvas("cdiff","cdiff",5,5,1600,800);
0538    cdiff->Divide(3,2);
0539    TH1D *hxdiff[5];
0540
0541
0542    hxdiff[0] = hxdiff_KF_truth_xglobal->ProjectionY();
0543    hxdiff[0]->SetTitle("global KF filt x - truth x");
0544    hxdiff[0]->GetXaxis()->SetTitle("global KF filt x - truth x (mm)");
0545
0546    hxdiff[1] = hxdiff_KF_truth_yglobal->ProjectionY();
0547    hxdiff[1]->SetTitle("global KF filt y - truth y");
0548    hxdiff[1]->GetXaxis()->SetTitle("global KF filt y - truth y (mm)");
0549
0550    hxdiff[2] = hxdiff_KF_truth_zglobal->ProjectionY();
0551    hxdiff[2]->SetTitle("global KF filt z - truth z");
0552    hxdiff[2]->GetXaxis()->SetTitle("global KF filt z - truth z (mm)");
0553
0554    hxdiff[3] = hxdiff_KF_truth_xlocal->ProjectionY();
0555    hxdiff[3]->SetTitle("local KF filt x - truth x");
0556    hxdiff[3]->GetXaxis()->SetTitle("local KF filt x - truth x (mm)");
0557
0558    hxdiff[4] = hxdiff_KF_truth_ylocal->ProjectionY();
0559    hxdiff[4]->SetTitle("local KF filt y - truth y");
0560    hxdiff[4]->GetXaxis()->SetTitle("local KF filt y - truth y (mm)");
0561
0562
0563    for(int i=0;i<5;++i)
0564      {
0565        cdiff->cd(i+1);
0566
0567        hxdiff[i]->SetTitleSize(0.05);
0568        hxdiff[i]->GetYaxis()->SetTitleSize(0.05);
0569        hxdiff[i]->GetXaxis()->SetTitleSize(0.05);
0570        hxdiff[i]->GetYaxis()->SetLabelSize(0.05);
0571        hxdiff[i]->GetXaxis()->SetLabelSize(0.05);
0572        hxdiff[i]->GetXaxis()->SetNdivisions(505);
0573        hxdiff[i]->DrawCopy();
0574        cout << " hxdiff " << i << " RMS =  " << hxdiff[i]->GetRMS() << endl;
0575      }
0576
0577
0578    TLegend *txleg = new TLegend(0.15, 0.75, 0.85, 0.85, "", "NDC");
0579    txleg->AddEntry(hxdiff[0],"KF  - truth  outer", "l");
0580    //txleg->AddEntry(hxerr_KF, "KF  error outer","l");
0581    txleg->SetTextSize(0.05);
0582    //   txleg->Draw();
0583
0584    TCanvas *chit = new TCanvas("chit","chit",5,5,1600,800);
0585    chit->Divide(3,2);
0586
0587    for(int i=0;i<5;++i)
0588      {
0589        chit->cd(i+1);
0590        hhit[i]->DrawCopy();
0591
0592      }
0593
0594    TCanvas *cpull = new TCanvas("cpull","local X/Y Pull",5,5,1600,800);
0595    cpull->Divide(2,1);
0596    cpull->cd(1);
0597    hxpull_KF->Draw();
0598    cpull->cd(2);
0599    hypull_KF->Draw();
0600
0601
0602    TCanvas *cdca = new TCanvas("cdca","DCA",5,5,1200,800);
0603    cdca->Divide(2,1);
0604
0605    TLegend *dcaleg = new TLegend(0.15, 0.75, 0.45, 0.85, "", "NDC");
0606    dcaleg->AddEntry(hdcaxy[0],"MVTX match", "l");
0607    dcaleg->AddEntry(hdcaxy[1],"no MVTX match", "l");
0608    dcaleg->SetTextSize(0.03);
0609
0610    cdca->cd(1);
0611    hdcaxy[0]->Draw();
0612    hdcaxy[1]->SetLineColor(kRed);
0613    hdcaxy[1]->Draw("same");
0614    dcaleg->Draw();
0615    cdca->cd(2);
0616    hdcaz[0]->Draw();
0617    hdcaz[1]->SetLineColor(kRed);
0618    hdcaz[1]->Draw("same");
0619    dcaleg->Draw();
0620
0621    TCanvas *cproto = new TCanvas("cproto","proto pos - truth pos",5,5,1200,800);
0622    cproto->Divide(3,1);
0623    cproto->cd(1);
0624    hprotox[0]->GetXaxis()->SetTitle("mm");
0625    hprotox[0]->Draw();
0626    hprotox[1]->SetLineColor(kRed);
0627    hprotox[1]->Draw("same");
0628    dcaleg->Draw();
0629
0630    cproto->cd(2);
0631    hprotoy[0]->GetXaxis()->SetTitle("mm");
0632    hprotoy[0]->Draw();
0633    hprotoy[1]->SetLineColor(kRed);
0634    hprotoy[1]->Draw("same");
0635    dcaleg->Draw();
0636
0637    cproto->cd(3);
0638    hprotoz[0]->GetXaxis()->SetTitle("mm");
0639    hprotoz[0]->Draw();
0640    hprotoz[1]->SetLineColor(kRed);
0641    hprotoz[1]->Draw("same");
0642    dcaleg->Draw();
0643
0644    TCanvas *citr = new TCanvas("citr","track number",5,5,1200,800);
0645    hitr[1]->SetMinimum(0);
0646    hitr[1]->Draw();
0647    hitr[0]->SetLineColor(kRed);
0648    hitr[0]->Draw("same");
0649
0650 }