File indexing completed on 2025-07-01 08:59:06
0001
0002
0003
0004 void hit_based_matching(int energy_set = 0, int Q2_set = 1, bool use_campaign = 0){
0005
0006
0007 gStyle->SetOptStat(0);
0008 gStyle->SetPadBorderMode(0);
0009 gStyle->SetFrameBorderMode(0);
0010 gStyle->SetFrameLineWidth(2);
0011 gStyle->SetLabelSize(0.035,"X");
0012 gStyle->SetLabelSize(0.035,"Y");
0013
0014
0015 gStyle->SetTitleXSize(0.04);
0016 gStyle->SetTitleXOffset(0.9);
0017 gStyle->SetTitleYSize(0.04);
0018 gStyle->SetTitleYOffset(0.9);
0019
0020
0021
0022
0023 TH1 *h1a = new TH1D("h1a","Truth-seeded tracks",100,0,1.1);
0024 h1a->GetXaxis()->SetTitle("Fraction of track measurements from a given MC Particle");h1a->GetXaxis()->CenterTitle();
0025 h1a->GetYaxis()->SetTitle("Number of Tracks"); h1a->GetYaxis()->CenterTitle();
0026 h1a->SetLineColor(kBlack);h1a->SetLineWidth(2);
0027
0028 TH1 *h1b = new TH1D("h1b","Real-seeded tracks",100,0,1.1);
0029 h1b->GetXaxis()->SetTitle("Fraction of track measurements from a given MC Particle");h1b->GetXaxis()->CenterTitle();
0030 h1b->GetYaxis()->SetTitle("Number of Tracks"); h1b->GetYaxis()->CenterTitle();
0031 h1b->SetLineColor(kBlack);h1b->SetLineWidth(2);
0032
0033
0034 TH2 *h2a = new TH2D("h2a","Truth-seeded tracks",100,0,50,100,0,50);
0035 h2a->GetXaxis()->SetTitle("Associated generated particle momentum [GeV/c]"); h2a->GetXaxis()->CenterTitle();
0036 h2a->GetYaxis()->SetTitle("Reconstructed track momentum [GeV/c]"); h2a->GetYaxis()->CenterTitle();
0037
0038 TH2 *h2b = new TH2D("h2b","Real-seeded tracks",100,0,50,100,0,50);
0039 h2b->GetXaxis()->SetTitle("Associated generated particle momentum [GeV/c]"); h2b->GetXaxis()->CenterTitle();
0040 h2b->GetYaxis()->SetTitle("Reconstructed track momentum [GeV/c]"); h2b->GetYaxis()->CenterTitle();
0041
0042
0043 TH2 *h3a = new TH2D("h3a","Truth-seeded tracks",100,0,50,100,0,50);
0044 h3a->GetXaxis()->SetTitle("Associated generated particle momentum [GeV/c]"); h3a->GetXaxis()->CenterTitle();
0045 h3a->GetYaxis()->SetTitle("Reconstructed track momentum [GeV/c]"); h3a->GetYaxis()->CenterTitle();
0046
0047 TH2 *h3b = new TH2D("h3b","Real-seeded tracks",100,0,50,100,0,50);
0048 h3b->GetXaxis()->SetTitle("Associated generated particle momentum [GeV/c]"); h3b->GetXaxis()->CenterTitle();
0049 h3b->GetYaxis()->SetTitle("Reconstructed track momentum [GeV/c]"); h3b->GetYaxis()->CenterTitle();
0050
0051
0052 TH1 *h4a = new TH1D("h4a","Truth-seeded tracks: Entirely pure tracks",2,0,2);
0053 h4a->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0054 h4a->GetXaxis()->SetBinLabel(2,"Primary Particle");
0055 h4a->GetXaxis()->SetTitle("Status of associated MC particle");h4a->GetXaxis()->CenterTitle();
0056 h4a->GetYaxis()->SetTitle("Number of Tracks"); h4a->GetYaxis()->CenterTitle();
0057 h4a->SetLineColor(kBlack);h4a->SetLineWidth(2);
0058
0059 TH1 *h4b = new TH1D("h4b","Real-seeded tracks: Entirely pure tracks",2,0,2);
0060 h4b->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0061 h4b->GetXaxis()->SetBinLabel(2,"Primary Particle");
0062 h4b->GetXaxis()->SetTitle("Status of associated MC particle");h4b->GetXaxis()->CenterTitle();
0063 h4b->GetYaxis()->SetTitle("Number of Tracks"); h4b->GetYaxis()->CenterTitle();
0064 h4b->SetLineColor(kBlack);h4b->SetLineWidth(2);
0065
0066
0067 TH1 *h5a = new TH1D("h5a","Truth-seeded tracks: Status of most-associated particle for mixed tracks",2,0,2);
0068 h5a->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0069 h5a->GetXaxis()->SetBinLabel(2,"Primary Particle");
0070 h5a->GetXaxis()->SetTitle("Status of associated MC particle");h5a->GetXaxis()->CenterTitle();
0071 h5a->GetYaxis()->SetTitle("Number of Tracks"); h5a->GetYaxis()->CenterTitle();
0072 h5a->SetLineColor(kBlack);h5a->SetLineWidth(2);
0073
0074 TH1 *h6a = new TH1D("h6a","Truth-seeded tracks: Status of less-associated particles for mixed tracks",2,0,2);
0075 h6a->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0076 h6a->GetXaxis()->SetBinLabel(2,"Primary Particle");
0077 h6a->GetXaxis()->SetTitle("Status of associated MC particle");h6a->GetXaxis()->CenterTitle();
0078 h6a->GetYaxis()->SetTitle("Number of Tracks"); h6a->GetYaxis()->CenterTitle();
0079 h6a->SetLineColor(kBlack);h6a->SetLineWidth(2);
0080
0081 TH1 *h5b = new TH1D("h5b","Real-seeded tracks: Status of most-associated particle for mixed tracks",2,0,2);
0082 h5b->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0083 h5b->GetXaxis()->SetBinLabel(2,"Primary Particle");
0084 h5b->GetXaxis()->SetTitle("Status of associated MC particle");h5b->GetXaxis()->CenterTitle();
0085 h5b->GetYaxis()->SetTitle("Number of Tracks"); h5b->GetYaxis()->CenterTitle();
0086 h5b->SetLineColor(kBlack);h5b->SetLineWidth(2);
0087
0088 TH1 *h6b = new TH1D("h6b","Real-seeded tracks: Status of less-associated particles for mixed tracks",2,0,2);
0089 h6b->GetXaxis()->SetBinLabel(1,"Secondary Particle");
0090 h6b->GetXaxis()->SetBinLabel(2,"Primary Particle");
0091 h6b->GetXaxis()->SetTitle("Status of associated MC particle");h6b->GetXaxis()->CenterTitle();
0092 h6b->GetYaxis()->SetTitle("Number of Tracks"); h6b->GetYaxis()->CenterTitle();
0093 h6b->SetLineColor(kBlack);h6b->SetLineWidth(2);
0094
0095
0096 TString path;
0097 if(use_campaign) path = "./campaign_input/";
0098 else path = "./input/";
0099
0100 TString beam_energies;
0101 if(energy_set == 0) beam_energies = "18x275";
0102 else if(energy_set == 1) beam_energies = "10x100";
0103
0104 TString run_name;
0105 if(use_campaign) run_name.Form("list_%s_Q2_%d.txt",beam_energies.Data(),Q2_set);
0106 else run_name.Form("%s/Q2_%d/eicrecon_*root",beam_energies.Data(),Q2_set);
0107
0108
0109 TString input = path + run_name;
0110 TChain *tree = new TChain("events");
0111
0112 if(use_campaign){
0113 std::ifstream in(input);
0114 std::string file("");
0115 while (in >> file)
0116 tree->Add(file.data());
0117 in.close();
0118 }
0119 else{
0120 tree->Add(input.Data());
0121 }
0122
0123 cout<<"Running file "<<input<<"!"<<endl;
0124 cout<<"Analyzing "<<tree->GetEntries()<<" events!"<<endl;
0125
0126
0127 TTreeReader tr(tree);
0128
0129 TTreeReaderArray<int> gen_status(tr, "MCParticles.generatorStatus");
0130 TTreeReaderArray<int> gen_pid(tr, "MCParticles.PDG");
0131 TTreeReaderArray<double> gen_px(tr, "MCParticles.momentum.x");
0132 TTreeReaderArray<double> gen_py(tr, "MCParticles.momentum.y");
0133 TTreeReaderArray<double> gen_pz(tr, "MCParticles.momentum.z");
0134 TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass");
0135 TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
0136 TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
0137 TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
0138 TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
0139
0140 TTreeReaderArray<float> hit_assoc_weight_ts(tr,"CentralCKFTruthSeededTrackAssociations.weight");
0141 TTreeReaderArray<int> hit_assoc_mc_ts(tr, "_CentralCKFTruthSeededTrackAssociations_sim.index");
0142 TTreeReaderArray<int> hit_assoc_track_ts(tr, "_CentralCKFTruthSeededTrackAssociations_rec.index");
0143
0144 TTreeReaderArray<float> hit_assoc_weight_rs(tr,"CentralCKFTrackAssociations.weight");
0145 TTreeReaderArray<int> hit_assoc_mc_rs(tr, "_CentralCKFTrackAssociations_sim.index");
0146 TTreeReaderArray<int> hit_assoc_track_rs(tr, "_CentralCKFTrackAssociations_rec.index");
0147
0148 TTreeReaderArray<float> track_qoverp_ts(tr, "CentralCKFTruthSeededTrackParameters.qOverP");
0149 TTreeReaderArray<float> track_qoverp_rs(tr, "CentralCKFTrackParameters.qOverP");
0150
0151
0152 TLorentzVector gen_vec;
0153 TVector3 gen_vertex;
0154
0155
0156 int pure_track_counter_ts(0);
0157 int mixed_track_counter_ts(0);
0158
0159
0160 int pure_track_counter_rs(0);
0161 int mixed_track_counter_rs(0);
0162
0163 TLorentzVector rec_vec;
0164 int counter(0);
0165
0166
0167 while (tr.Next()) {
0168
0169 if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
0170 counter++;
0171
0172
0173 std::vector<std::tuple<int, int, float>> vec;
0174
0175
0176 for(int iassoc=0;iassoc<hit_assoc_weight_ts.GetSize();iassoc++){
0177
0178 auto assoc_weight = hit_assoc_weight_ts[iassoc];
0179
0180 h1a->Fill(assoc_weight);
0181
0182
0183 if(assoc_weight > 0.99)
0184 pure_track_counter_ts++;
0185 else if(iassoc==0)
0186 mixed_track_counter_ts++;
0187 else if(iassoc>0 && hit_assoc_track_ts[iassoc] != hit_assoc_track_ts[iassoc-1])
0188 mixed_track_counter_ts++;
0189
0190
0191
0192 int mc_index = hit_assoc_mc_ts[iassoc];
0193 int trk_index = hit_assoc_track_ts[iassoc];
0194
0195 gen_vec.SetXYZM(gen_px[mc_index],gen_py[mc_index],gen_pz[mc_index],gen_mass[mc_index]);
0196 auto gen_mom = gen_vec.P();
0197
0198 auto trk_mom = fabs(1./track_qoverp_ts[trk_index]);
0199
0200 h2a->Fill(gen_mom,trk_mom);
0201 if( fabs(gen_vec.Eta()) < 3 )
0202 h3a->Fill(gen_mom,trk_mom);
0203
0204
0205 int status = gen_status[mc_index];
0206
0207 if(assoc_weight > 0.99){
0208 h4a->Fill(status);
0209 }
0210 else{
0211 vec.push_back(std::make_tuple(trk_index,status,assoc_weight));
0212 }
0213
0214 }
0215
0216
0217 std::sort(vec.begin(), vec.end());
0218
0219
0220 size_t ivec = 0;
0221 while (ivec < vec.size()) {
0222
0223 int current_first = std::get<0>(vec[ivec]);
0224 float current_max_weight = std::get<2>(vec[ivec]);
0225 size_t start = ivec;
0226 size_t index_max_weight = ivec;
0227
0228
0229 while (ivec < vec.size() && std::get<0>(vec[ivec]) == current_first) {
0230 float current_weight = std::get<2>(vec[ivec]);
0231 if( current_weight > current_max_weight ){
0232 current_max_weight = current_weight;
0233 index_max_weight = ivec;
0234 }
0235 ++ivec;
0236 }
0237
0238
0239 if (ivec - start > 1) {
0240 for (size_t ii = start; ii < ivec; ++ii){
0241
0242 if(ii == index_max_weight)
0243 h5a->Fill( std::get<1>(vec[ii]) );
0244 else
0245 h6a->Fill( std::get<1>(vec[ii]) );
0246 }
0247 }
0248 }
0249
0250
0251 vec.clear();
0252 ivec = 0;
0253
0254
0255
0256 for(int iassoc=0;iassoc<hit_assoc_weight_rs.GetSize();iassoc++){
0257
0258 auto assoc_weight = hit_assoc_weight_rs[iassoc];
0259
0260 h1b->Fill(assoc_weight);
0261
0262
0263 if(assoc_weight > 0.99)
0264 pure_track_counter_rs++;
0265 else if(iassoc==0)
0266 mixed_track_counter_rs++;
0267 else if(iassoc>0 && hit_assoc_track_rs[iassoc] != hit_assoc_track_rs[iassoc-1])
0268 mixed_track_counter_rs++;
0269
0270
0271
0272 int mc_index = hit_assoc_mc_rs[iassoc];
0273 int trk_index = hit_assoc_track_rs[iassoc];
0274
0275 gen_vec.SetXYZM(gen_px[mc_index],gen_py[mc_index],gen_pz[mc_index],gen_mass[mc_index]);
0276 auto gen_mom = gen_vec.P();
0277
0278 auto trk_mom = fabs(1./track_qoverp_rs[trk_index]);
0279
0280 h2b->Fill(gen_mom,trk_mom);
0281
0282 if( fabs(gen_vec.Eta()) < 3 )
0283 h3b->Fill(gen_mom,trk_mom);
0284
0285
0286 int status = gen_status[mc_index];
0287
0288 if(assoc_weight > 0.99){
0289 h4b->Fill(status);
0290 }
0291 else{
0292 vec.push_back(std::make_tuple(trk_index,status,assoc_weight));
0293 }
0294
0295 }
0296
0297
0298 std::sort(vec.begin(), vec.end());
0299
0300
0301 while (ivec < vec.size()) {
0302
0303 int current_first = std::get<0>(vec[ivec]);
0304 float current_max_weight = std::get<2>(vec[ivec]);
0305 size_t start = ivec;
0306 size_t index_max_weight = ivec;
0307
0308
0309 while (ivec < vec.size() && std::get<0>(vec[ivec]) == current_first) {
0310 float current_weight = std::get<2>(vec[ivec]);
0311 if( current_weight > current_max_weight ){
0312 current_max_weight = current_weight;
0313 index_max_weight = ivec;
0314 }
0315 ++ivec;
0316 }
0317
0318
0319 if (ivec - start > 1) {
0320 for (size_t ii = start; ii < ivec; ++ii){
0321
0322 if(ii == index_max_weight)
0323 h5b->Fill( std::get<1>(vec[ii]) );
0324 else
0325 h6b->Fill( std::get<1>(vec[ii]) );
0326 }
0327 }
0328 }
0329
0330
0331 }
0332
0333
0334 float frac_pure_ts = ((float) pure_track_counter_ts) / ((float) (pure_track_counter_ts + mixed_track_counter_ts));
0335 float frac_pure_rs = ((float) pure_track_counter_rs) / ((float) (pure_track_counter_rs + mixed_track_counter_rs));
0336
0337
0338
0339
0340 TCanvas *c1a = new TCanvas("c1a");
0341 c1a->SetLogy();
0342 h1a->Draw();
0343
0344 TLegend *leg1a = new TLegend(0.125,0.7,0.475,0.875);
0345 leg1a->SetBorderSize(0);leg1a->SetFillStyle(0);
0346 leg1a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0347 leg1a->AddEntry(h1a, Form("Fraction of entirely pure tracks = %.3f",frac_pure_ts),"l");
0348 leg1a->Draw();
0349
0350
0351 TCanvas *c1b = new TCanvas("c1b");
0352 c1b->SetLogy();
0353 h1b->Draw();
0354
0355 TLegend *leg1b = new TLegend(0.125,0.7,0.475,0.875);
0356 leg1b->SetBorderSize(0);leg1b->SetFillStyle(0);
0357 leg1b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0358 leg1b->AddEntry(h1b, Form("Fraction of entirely pure tracks = %.3f",frac_pure_rs),"l");
0359 leg1b->Draw();
0360
0361
0362 TCanvas *c2a = new TCanvas("c2a");
0363 c2a->SetLogz();
0364 h2a->Draw("colz");
0365
0366 TLatex *tex2a = new TLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0367 tex2a->SetTextSize(0.03);
0368 tex2a->Draw();
0369
0370
0371 TCanvas *c2b = new TCanvas("c2b");
0372 c2b->SetLogz();
0373 h2b->Draw("colz");
0374
0375 TLatex *tex2b = new TLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0376 tex2b->SetTextSize(0.03);
0377 tex2b->Draw();
0378
0379
0380 TCanvas *c3a = new TCanvas("c3a");
0381 c3a->SetLogz();
0382 h3a->Draw("colz");
0383
0384 TLatex *tex3a = new TLatex();
0385 tex3a->SetTextSize(0.03);
0386 tex3a->DrawLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0387 tex3a->DrawLatex(25,5," -3 < #eta_{gen} < +3");
0388
0389
0390 TCanvas *c3b = new TCanvas("c3b");
0391 c3b->SetLogz();
0392 h3b->Draw("colz");
0393
0394 TLatex *tex3b = new TLatex();
0395 tex3b->SetTextSize(0.03);
0396 tex3b->DrawLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0397 tex3b->DrawLatex(25,5," -3 < #eta_{gen} < +3");
0398
0399
0400 TCanvas *c4a = new TCanvas("c4a");
0401 c4a->SetLogy();
0402 h4a->Draw();
0403
0404 TLegend *leg4a = new TLegend(0.125,0.7,0.475,0.875);
0405 leg4a->SetBorderSize(0);leg4a->SetFillStyle(0);
0406 leg4a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0407 leg4a->Draw();
0408
0409
0410 TCanvas *c4b = new TCanvas("c4b");
0411 c4b->SetLogy();
0412 h4b->Draw();
0413
0414 TLegend *leg4b = new TLegend(0.125,0.7,0.475,0.875);
0415 leg4b->SetBorderSize(0);leg4b->SetFillStyle(0);
0416 leg4b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0417 leg4b->Draw();
0418
0419
0420 TCanvas *c5a = new TCanvas("c5a");
0421 c5a->SetLogy();
0422 h5a->Draw();
0423
0424 TLegend *leg5a = new TLegend(0.125,0.7,0.475,0.875);
0425 leg5a->SetBorderSize(0);leg5a->SetFillStyle(0);
0426 leg5a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0427 leg5a->Draw();
0428
0429 TCanvas *c6a = new TCanvas("c6a");
0430 c6a->SetLogy();
0431 h6a->Draw();
0432
0433 TLegend *leg6a = new TLegend(0.525,0.7,0.875,0.875);
0434 leg6a->SetBorderSize(0);leg6a->SetFillStyle(0);
0435 leg6a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0436 leg6a->Draw();
0437
0438
0439 TCanvas *c5b = new TCanvas("c5b");
0440 c5b->SetLogy();
0441 h5b->Draw();
0442
0443 TLegend *leg5b = new TLegend(0.125,0.7,0.475,0.875);
0444 leg5b->SetBorderSize(0);leg5b->SetFillStyle(0);
0445 leg5b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0446 leg5b->Draw();
0447
0448 TCanvas *c6b = new TCanvas("c6b");
0449 c6b->SetLogy();
0450 h6b->Draw();
0451
0452 TLegend *leg6b = new TLegend(0.525,0.7,0.875,0.875);
0453 leg6b->SetBorderSize(0);leg6b->SetFillStyle(0);
0454 leg6b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0455 leg6b->Draw();
0456
0457
0458 c1a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf[",beam_energies.Data(),Q2_set));
0459 c1a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0460 c1b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0461 c2a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0462 c2b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0463 c3a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0464 c3b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0465 c4a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0466 c4b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0467 c5a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0468 c6a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0469 c5b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0470 c6b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0471 c6b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf]",beam_energies.Data(),Q2_set));
0472
0473 }