File indexing completed on 2026-07-16 08:11:27
0001 #include "../corrections/GeoAcceptanceCorrection.h"
0002
0003 std::string combine_name(std::string base, std::string suffix)
0004 {
0005 if(suffix!="")
0006 {
0007 return base+"_"+suffix;
0008 }
0009 else
0010 {
0011 return base;
0012 }
0013 }
0014
0015 std::string combine_title(std::string base, std::string suffix, std::string axislabel)
0016 {
0017 if(suffix!="")
0018 {
0019 return base+" "+suffix+";"+axislabel;
0020 }
0021 else
0022 {
0023 return base+";"+axislabel;
0024 }
0025 }
0026
0027 struct HistContainer
0028 {
0029 TH1F* lifetime_truth;
0030 TH1F* lifetime_reco;
0031
0032 TH1F* decaylength_truth;
0033 TH1F* decaylength_reco;
0034
0035 TH1F* lifetime_truth_signal;
0036 TH1F* lifetime_truth_background;
0037 TH1F* lifetime_reco_signal;
0038 TH1F* lifetime_reco_background;
0039
0040 HistContainer(std::string name_suffix, std::string title_suffix,
0041 int nBins_t, float min_t, float max_t,
0042 int nBins_L, float min_L, float max_L)
0043 {
0044 const std::string lifetime_label = "#tau [ps]";
0045 const std::string decaylength_label = "decay length [cm]";
0046
0047 lifetime_truth = new TH1F(combine_name("lifetime_truth",name_suffix).c_str(),
0048 combine_title("Truth lifetime",title_suffix,lifetime_label).c_str(),
0049 nBins_t,min_t,max_t);
0050 lifetime_reco = new TH1F(combine_name("lifetime_reco",name_suffix).c_str(),
0051 combine_title("Reco lifetime",title_suffix,lifetime_label).c_str(),
0052 nBins_t,min_t,max_t);
0053 decaylength_truth = new TH1F(combine_name("decaylength_truth",name_suffix).c_str(),
0054 combine_title("Truth decay length",title_suffix,decaylength_label).c_str(),
0055 nBins_L,min_L,max_L);
0056 decaylength_reco = new TH1F(combine_name("decaylength_reco",name_suffix).c_str(),
0057 combine_title("Reco decay length",title_suffix,decaylength_label).c_str(),
0058 nBins_L,min_L,max_L);
0059 lifetime_truth_signal = new TH1F(combine_name("lifetime_truth_signal",name_suffix).c_str(),
0060 combine_title("Signal truth lifetime",title_suffix,lifetime_label).c_str(),
0061 nBins_t,min_t,max_t);
0062 lifetime_truth_background = new TH1F(combine_name("lifetime_truth_background",name_suffix).c_str(),
0063 combine_title("Background truth lifetime",title_suffix,lifetime_label).c_str(),
0064 nBins_t,min_t,max_t);
0065 lifetime_reco_signal = new TH1F(combine_name("lifetime_reco_signal",name_suffix).c_str(),
0066 combine_title("Signal reco lifetime",title_suffix,lifetime_label).c_str(),
0067 nBins_t,min_t,max_t);
0068 lifetime_reco_background = new TH1F(combine_name("lifetime_reco_background",name_suffix).c_str(),
0069 combine_title("Background reco lifetime",title_suffix,lifetime_label).c_str(),
0070 nBins_t,min_t,max_t);
0071 }
0072
0073 void Write()
0074 {
0075 lifetime_truth->Write();
0076 lifetime_reco->Write();
0077 decaylength_truth->Write();
0078 decaylength_reco->Write();
0079 lifetime_truth_signal->Write();
0080 lifetime_truth_background->Write();
0081 lifetime_reco_signal->Write();
0082 lifetime_reco_background->Write();
0083 }
0084 };
0085
0086 void KFParticle_truth_lifetime()
0087 {
0088 const bool verbose = false;
0089 const float min_DIRA = 0.99;
0090 const float min_track_IPxy = 0.05;
0091 const float min_decayLength = 0.05;
0092 const float min_mother_pt = 0.8;
0093 const int min_nMVTX = 1;
0094 const int min_nINTT = 1;
0095 const int min_nTPC = 1;
0096 const float max_pv_quality = 20.;
0097 const float max_daughter_DCA = 0.5;
0098
0099 const std::string particle_name = "Lambda0";
0100 const int mother_pdgid = 3122;
0101 const float mother_mass = 1.115;
0102
0103 const float speed = 2.99792458e-2;
0104
0105 GeoAcceptanceCorrection geoacc_corr("/sphenix/u/cdean/analysis/LightFlavorRatios/geometric_acceptance/analysis/plots/Lambda0_geometric_acceptance_ratio_pT.root","Lambda0_inGeo_pT");
0106
0107
0108 TFile* f = TFile::Open("/sphenix/tg/tg01/hf/mjpeters/LightFlavorResults/outputKFParticle_Lambda_reco_swum_weighted.root");
0109 TTree* t = (TTree*)f->Get("DecayTree");
0110
0111 TFile* fout = new TFile("KFParticle_truth_h.root","RECREATE");
0112
0113 const int nBins_lifetime = 100;
0114 const float min_lifetime = 0.;
0115 const float max_lifetime = 500.;
0116
0117 const int nBins_decaylength = 100;
0118 const float min_decaylength = 0.;
0119 const float max_decaylength = 10.;
0120
0121
0122 HistContainer h("","",
0123 nBins_lifetime,min_lifetime,max_lifetime,nBins_decaylength,min_decaylength,max_decaylength);
0124
0125 HistContainer h_motherpt("motherptcut","with mother pT>800MeV",
0126 nBins_lifetime,min_lifetime,max_lifetime,nBins_decaylength,min_decaylength,max_decaylength);
0127
0128 HistContainer h_geoacc("geoacc","with mother pT>800MeV, geo. acc. corrected",
0129 nBins_lifetime,min_lifetime,max_lifetime,nBins_decaylength,min_decaylength,max_decaylength);
0130
0131 HistContainer h_kfcut("kfcut","with mother pT>800MeV, geo. acc. corrected, with KFParticle selection cuts",
0132 nBins_lifetime,min_lifetime,max_lifetime,nBins_decaylength,min_decaylength,max_decaylength);
0133
0134 HistContainer h_swim("swim","with mother pT>800 MeV, geo. acc. corrected, with KFParticle selection cuts and swimming correction",
0135 nBins_lifetime,min_lifetime,max_lifetime,nBins_decaylength,min_decaylength,max_decaylength);
0136
0137 TH1F* lifetime_diff = new TH1F("lifetime_diff","lifetime_diff",100,-20.,20.);
0138 TH1F* decaylength_diff = new TH1F("decaylength_diff","decaylength_diff",100,-2.,2.);
0139
0140 float decayTime;
0141 float decayLength;
0142 float mother_p;
0143 float mother_pt;
0144 float mother_DIRA;
0145 float mother_PV_x;
0146 float mother_PV_y;
0147 float mother_PV_z;
0148
0149 std::vector<int>* track1_pdgid_history = nullptr;
0150 float track1_EV_x;
0151 float track1_EV_y;
0152 float track1_EV_z;
0153 float track1_PV_x;
0154 float track1_PV_y;
0155 float track1_PV_z;
0156 float track1_gpx;
0157 float track1_gpy;
0158 float track1_gpz;
0159 float track1_IP_xy;
0160 unsigned int track1_nMVTX;
0161 unsigned int track1_nINTT;
0162 unsigned int track1_nTPC;
0163
0164 std::vector<int>* track2_pdgid_history = nullptr;
0165 float track2_EV_x;
0166 float track2_EV_y;
0167 float track2_EV_z;
0168 float track2_PV_x;
0169 float track2_PV_y;
0170 float track2_PV_z;
0171 float track2_gpx;
0172 float track2_gpy;
0173 float track2_gpz;
0174 float track2_IP_xy;
0175 unsigned int track2_nMVTX;
0176 unsigned int track2_nINTT;
0177 unsigned int track2_nTPC;
0178
0179 float pv_chi2;
0180 unsigned int pv_ndof;
0181 float daughter_DCA;
0182
0183 std::vector<float>* turningPoints_TAU = nullptr;
0184 float swim_weight;
0185
0186 t->SetBranchAddress((particle_name+"_decayTime").c_str(),&decayTime);
0187 t->SetBranchAddress((particle_name+"_decayLength").c_str(),&decayLength);
0188 t->SetBranchAddress((particle_name+"_p").c_str(),&mother_p);
0189 t->SetBranchAddress((particle_name+"_pT").c_str(),&mother_pt);
0190 t->SetBranchAddress((particle_name+"_DIRA").c_str(),&mother_DIRA);
0191 t->SetBranchAddress("primary_vertex_x",&mother_PV_x);
0192 t->SetBranchAddress("primary_vertex_y",&mother_PV_y);
0193 t->SetBranchAddress("primary_vertex_z",&mother_PV_z);
0194
0195 t->SetBranchAddress("track_1_true_track_history_PDG_ID",&track1_pdgid_history);
0196 t->SetBranchAddress("track_1_true_EV_x",&track1_EV_x);
0197 t->SetBranchAddress("track_1_true_EV_y",&track1_EV_y);
0198 t->SetBranchAddress("track_1_true_EV_z",&track1_EV_z);
0199 t->SetBranchAddress("track_1_true_PV_x",&track1_PV_x);
0200 t->SetBranchAddress("track_1_true_PV_y",&track1_PV_y);
0201 t->SetBranchAddress("track_1_true_PV_z",&track1_PV_z);
0202 t->SetBranchAddress("track_1_true_px",&track1_gpx);
0203 t->SetBranchAddress("track_1_true_py",&track1_gpy);
0204 t->SetBranchAddress("track_1_true_pz",&track1_gpz);
0205 t->SetBranchAddress("track_1_IP_xy",&track1_IP_xy);
0206 t->SetBranchAddress("track_1_MVTX_nHits",&track1_nMVTX);
0207 t->SetBranchAddress("track_1_INTT_nHits",&track1_nINTT);
0208 t->SetBranchAddress("track_1_TPC_nHits",&track1_nTPC);
0209
0210 t->SetBranchAddress("track_2_true_track_history_PDG_ID",&track2_pdgid_history);
0211 t->SetBranchAddress("track_2_true_EV_x",&track2_EV_x);
0212 t->SetBranchAddress("track_2_true_EV_y",&track2_EV_y);
0213 t->SetBranchAddress("track_2_true_EV_z",&track2_EV_z);
0214 t->SetBranchAddress("track_2_true_PV_x",&track2_PV_x);
0215 t->SetBranchAddress("track_2_true_PV_y",&track2_PV_y);
0216 t->SetBranchAddress("track_2_true_PV_z",&track2_PV_z);
0217 t->SetBranchAddress("track_2_true_px",&track2_gpx);
0218 t->SetBranchAddress("track_2_true_py",&track2_gpy);
0219 t->SetBranchAddress("track_2_true_pz",&track2_gpz);
0220 t->SetBranchAddress("track_2_IP_xy",&track2_IP_xy);
0221 t->SetBranchAddress("track_2_MVTX_nHits",&track2_nMVTX);
0222 t->SetBranchAddress("track_2_INTT_nHits",&track2_nINTT);
0223 t->SetBranchAddress("track_2_TPC_nHits",&track2_nTPC);
0224
0225 t->SetBranchAddress("primary_vertex_chi2",&pv_chi2);
0226 t->SetBranchAddress("primary_vertex_ndof",&pv_ndof);
0227 t->SetBranchAddress("track_1_track_2_DCA",&daughter_DCA);
0228
0229 t->SetBranchAddress("swim_weight",&swim_weight);
0230 t->SetBranchAddress("turningPoints_TAU",&turningPoints_TAU);
0231
0232 for(int i=0;i<t->GetEntries();i++)
0233 {
0234 if(i % 1000 == 0) std::cout << i << std::endl;
0235 t->GetEntry(i);
0236
0237 if(track1_nMVTX<min_nMVTX || track2_nMVTX<min_nMVTX ||
0238 track1_nINTT<min_nINTT || track2_nINTT<min_nINTT ||
0239 track1_nTPC<min_nTPC || track2_nTPC<min_nTPC ||
0240 pv_chi2/pv_ndof>max_pv_quality ||
0241 daughter_DCA>max_daughter_DCA)
0242 {
0243 continue;
0244 }
0245
0246 if(verbose)
0247 {
0248 std::cout << "--------------------------------------------" << std::endl;
0249 std::cout << "entry " << i << std::endl;
0250 std::cout << "reco decayLength: " << decayLength << " decayTime: " << decayTime << std::endl;
0251 std::cout << "p: " << mother_p << std::endl;
0252 std::cout << "reco PV: (" << mother_PV_x << ", " << mother_PV_y << ", " << mother_PV_z << ")" << std::endl;
0253 }
0254
0255 const bool is_signal = (track1_pdgid_history && track2_pdgid_history &&
0256 track1_pdgid_history->size()>0 && track2_pdgid_history->size()>0 &&
0257 ((track1_pdgid_history->at(0)==mother_pdgid && track2_pdgid_history->at(0)==mother_pdgid) ||
0258 (track1_pdgid_history->at(0)==-mother_pdgid && track2_pdgid_history->at(0)==-mother_pdgid)));
0259
0260 if(verbose)
0261 {
0262 if(is_signal)
0263 {
0264 std::cout << "is signal" << std::endl;
0265 }
0266 else
0267 {
0268 std::cout << "is not signal" << std::endl;
0269 if(track1_pdgid_history == nullptr)
0270 {
0271 std::cout << "history unavailable for track 1" << std::endl;
0272 }
0273 if(track2_pdgid_history == nullptr)
0274 {
0275 std::cout << "history unavailable for track 2" << std::endl;
0276 }
0277 if(track1_pdgid_history && track2_pdgid_history && track1_pdgid_history->size()==0)
0278 {
0279 std::cout << "history size is zero for track 1" << std::endl;
0280 }
0281 if(track1_pdgid_history && track2_pdgid_history && track2_pdgid_history->size()==0)
0282 {
0283 std::cout << "history size is zero for track 2" << std::endl;
0284 }
0285 if(track1_pdgid_history && track2_pdgid_history && track1_pdgid_history->size()>0 && track2_pdgid_history->size()>0)
0286 {
0287 std::cout << "direct descendant PIDs are " << track1_pdgid_history->at(0) << " and " << track2_pdgid_history->at(0) << std::endl;
0288 }
0289 }
0290 }
0291
0292 std::array<float,3> true_SV;
0293
0294 if(fabs(track1_PV_x+99.)>0.001)
0295 {
0296 true_SV = {track1_PV_x, track1_PV_y, track1_PV_z};
0297 }
0298 else if(fabs(track2_PV_x+99.)>0.001)
0299 {
0300 true_SV = {track2_PV_x, track2_PV_y, track2_PV_z};
0301 }
0302 else
0303 {
0304 true_SV = {track1_EV_x,track1_EV_y,track1_EV_z};
0305 }
0306
0307 if(verbose)
0308 {
0309 std::cout << "SV: (" << true_SV[0] << ", " << true_SV[1] << ", " << true_SV[2] << ")" << std::endl;
0310 std::cout << "PV: (" << mother_PV_x << ", " << mother_PV_y << ", " << mother_PV_z << ")" << std::endl;
0311 }
0312
0313 const float true_decayLength = sqrt(pow(true_SV[0]-mother_PV_x,2)+pow(true_SV[1]-mother_PV_y,2)+pow(true_SV[2]-mother_PV_z,2));
0314 const float true_p = sqrt(pow(track1_gpx+track2_gpx,2)+pow(track1_gpy+track2_gpy,2)+pow(track1_gpz+track2_gpz,2));
0315 const float true_decayTime = 1./speed * mother_mass*true_decayLength/true_p;
0316
0317 if(verbose)
0318 {
0319 std::cout << "true decay length: " << true_decayLength << std::endl;
0320 std::cout << "true p: " << true_p << std::endl;
0321 std::cout << "true decay time: " << true_decayTime << std::endl;
0322 }
0323
0324 h.lifetime_reco->Fill(decayTime);
0325 h.decaylength_reco->Fill(decayLength);
0326 h.decaylength_truth->Fill(true_decayLength);
0327 decaylength_diff->Fill(decayLength-true_decayLength);
0328 h.lifetime_truth->Fill(true_decayTime);
0329 lifetime_diff->Fill(decayTime-true_decayTime);
0330
0331 if(is_signal)
0332 {
0333 h.lifetime_truth_signal->Fill(true_decayTime);
0334 h.lifetime_reco_signal->Fill(decayTime);
0335 }
0336 else
0337 {
0338 h.lifetime_truth_background->Fill(true_decayTime);
0339 h.lifetime_reco_background->Fill(decayTime);
0340 }
0341
0342 if(mother_pt>min_mother_pt)
0343 {
0344 h_motherpt.lifetime_truth->Fill(true_decayTime);
0345 h_motherpt.decaylength_truth->Fill(true_decayLength);
0346 h_motherpt.lifetime_reco->Fill(decayTime);
0347 h_motherpt.decaylength_reco->Fill(decayLength);
0348
0349 const float geoacc_weight = 1./geoacc_corr.h_corr->Interpolate(mother_pt);
0350 h_geoacc.lifetime_truth->Fill(true_decayTime,geoacc_weight);
0351 h_geoacc.decaylength_truth->Fill(true_decayLength,geoacc_weight);
0352 h_geoacc.lifetime_reco->Fill(decayTime,geoacc_weight);
0353 h_geoacc.decaylength_reco->Fill(decayLength,geoacc_weight);
0354
0355 if(is_signal)
0356 {
0357 h_motherpt.lifetime_truth_signal->Fill(true_decayTime);
0358 h_motherpt.lifetime_reco_signal->Fill(decayTime);
0359
0360 h_geoacc.lifetime_truth_signal->Fill(true_decayTime,geoacc_weight);
0361 h_geoacc.lifetime_reco_signal->Fill(decayTime,geoacc_weight);
0362 }
0363 else
0364 {
0365 h_motherpt.lifetime_truth_background->Fill(true_decayTime);
0366 h_motherpt.lifetime_reco_background->Fill(decayTime);
0367
0368 h_geoacc.lifetime_truth_background->Fill(true_decayTime,geoacc_weight);
0369 h_geoacc.lifetime_reco_background->Fill(decayTime,geoacc_weight);
0370 }
0371
0372 if(mother_DIRA>min_DIRA && decayLength>min_decayLength && fabs(track1_IP_xy)>min_track_IPxy && fabs(track2_IP_xy)>min_track_IPxy)
0373 {
0374 h_kfcut.lifetime_truth->Fill(true_decayTime,geoacc_weight);
0375 h_kfcut.decaylength_truth->Fill(true_decayLength,geoacc_weight);
0376 h_kfcut.lifetime_reco->Fill(decayTime,geoacc_weight);
0377 h_kfcut.decaylength_reco->Fill(decayLength,geoacc_weight);
0378
0379 h_swim.lifetime_truth->Fill(true_decayTime,geoacc_weight*swim_weight);
0380 h_swim.decaylength_truth->Fill(true_decayLength,geoacc_weight*swim_weight);
0381 h_swim.lifetime_reco->Fill(decayTime,geoacc_weight*swim_weight);
0382 h_swim.decaylength_reco->Fill(decayLength,geoacc_weight*swim_weight);
0383
0384 if(is_signal)
0385 {
0386 h_kfcut.lifetime_truth_signal->Fill(true_decayTime,geoacc_weight);
0387 h_kfcut.lifetime_reco_signal->Fill(decayTime,geoacc_weight);
0388
0389 h_swim.lifetime_truth_signal->Fill(true_decayTime,geoacc_weight*swim_weight);
0390 h_swim.lifetime_reco_signal->Fill(decayTime,geoacc_weight*swim_weight);
0391 }
0392 else
0393 {
0394 h_kfcut.lifetime_truth_background->Fill(true_decayTime,geoacc_weight);
0395 h_kfcut.lifetime_reco_background->Fill(decayTime,geoacc_weight);
0396
0397 h_swim.lifetime_truth_background->Fill(true_decayTime,geoacc_weight*swim_weight);
0398 h_swim.lifetime_reco_background->Fill(decayTime,geoacc_weight*swim_weight);
0399 }
0400 }
0401 }
0402 }
0403
0404 h.Write();
0405 h_motherpt.Write();
0406 h_geoacc.Write();
0407 h_kfcut.Write();
0408 h_swim.Write();
0409
0410 lifetime_diff->Write();
0411 decaylength_diff->Write();
0412 }
