File indexing completed on 2025-01-18 10:18:32
0001
0002
0003
0004 void hit_based_matching(int energy_set = 0, int Q2_set = 1){
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 = "./input/";
0097 TString run_name;
0098 TString beam_energies;
0099
0100 if(energy_set == 0 && Q2_set == 1){
0101 beam_energies = "18x275";
0102 run_name = "Q2_1/eicrecon_*.root";
0103 }
0104 else if(energy_set == 0 && Q2_set == 10){
0105 beam_energies = "18x275";
0106 run_name = "eicrecon_out_E_18_275_Q2_10.root";
0107 }
0108 else if(energy_set == 0 && Q2_set == 100){
0109 beam_energies = "18x275";
0110 run_name = "Q2_100/eicrecon_*.root";
0111 }
0112 else if(energy_set == 1 && Q2_set == 1){
0113 beam_energies = "10x100";
0114 run_name = "eicrecon_out_E_10_100_Q2_1.root";
0115 }
0116 else if(energy_set == 1 && Q2_set == 10){
0117 beam_energies = "10x100";
0118 run_name = "eicrecon_out_E_10_100_Q2_10.root";
0119 }
0120 else if(energy_set == 1 && Q2_set == 100){
0121 beam_energies = "10x100";
0122 run_name = "eicrecon_out_E_10_100_Q2_100.root";
0123 }
0124
0125
0126 TString input = path + run_name;
0127 TChain *tree = new TChain("events");
0128 tree->Add(input.Data());
0129
0130 cout<<"Analyzing "<<tree->GetEntries()<<" events!"<<endl;
0131
0132
0133 TTreeReader tr(tree);
0134
0135 TTreeReaderArray<int> gen_status(tr, "MCParticles.generatorStatus");
0136 TTreeReaderArray<int> gen_pid(tr, "MCParticles.PDG");
0137 TTreeReaderArray<float> gen_px(tr, "MCParticles.momentum.x");
0138 TTreeReaderArray<float> gen_py(tr, "MCParticles.momentum.y");
0139 TTreeReaderArray<float> gen_pz(tr, "MCParticles.momentum.z");
0140 TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass");
0141 TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
0142 TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
0143 TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
0144 TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
0145
0146 TTreeReaderArray<float> hit_assoc_weight_ts(tr,"CentralCKFTruthSeededTrackAssociations.weight");
0147 TTreeReaderArray<int> hit_assoc_mc_ts(tr, "_CentralCKFTruthSeededTrackAssociations_sim.index");
0148 TTreeReaderArray<int> hit_assoc_track_ts(tr, "_CentralCKFTruthSeededTrackAssociations_rec.index");
0149
0150 TTreeReaderArray<float> hit_assoc_weight_rs(tr,"CentralCKFTrackAssociations.weight");
0151 TTreeReaderArray<int> hit_assoc_mc_rs(tr, "_CentralCKFTrackAssociations_sim.index");
0152 TTreeReaderArray<int> hit_assoc_track_rs(tr, "_CentralCKFTrackAssociations_rec.index");
0153
0154 TTreeReaderArray<float> track_qoverp_ts(tr, "CentralCKFTruthSeededTrackParameters.qOverP");
0155 TTreeReaderArray<float> track_qoverp_rs(tr, "CentralCKFTrackParameters.qOverP");
0156
0157
0158 TLorentzVector gen_vec;
0159 TVector3 gen_vertex;
0160
0161
0162 int pure_track_counter_ts(0);
0163 int mixed_track_counter_ts(0);
0164
0165
0166 int pure_track_counter_rs(0);
0167 int mixed_track_counter_rs(0);
0168
0169 TLorentzVector rec_vec;
0170 int counter(0);
0171
0172
0173 while (tr.Next()) {
0174
0175 if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
0176 counter++;
0177
0178
0179 std::vector<std::tuple<int, int, float>> vec;
0180
0181
0182 for(int iassoc=0;iassoc<hit_assoc_weight_ts.GetSize();iassoc++){
0183
0184 auto assoc_weight = hit_assoc_weight_ts[iassoc];
0185
0186 h1a->Fill(assoc_weight);
0187
0188
0189 if(assoc_weight > 0.99)
0190 pure_track_counter_ts++;
0191 else if(iassoc==0)
0192 mixed_track_counter_ts++;
0193 else if(iassoc>0 && hit_assoc_track_ts[iassoc] != hit_assoc_track_ts[iassoc-1])
0194 mixed_track_counter_ts++;
0195
0196
0197
0198 int mc_index = hit_assoc_mc_ts[iassoc];
0199 int trk_index = hit_assoc_track_ts[iassoc];
0200
0201 gen_vec.SetXYZM(gen_px[mc_index],gen_py[mc_index],gen_pz[mc_index],gen_mass[mc_index]);
0202 auto gen_mom = gen_vec.P();
0203
0204 auto trk_mom = fabs(1./track_qoverp_ts[trk_index]);
0205
0206 h2a->Fill(gen_mom,trk_mom);
0207 if( fabs(gen_vec.Eta()) < 3 )
0208 h3a->Fill(gen_mom,trk_mom);
0209
0210
0211 int status = gen_status[mc_index];
0212
0213 if(assoc_weight > 0.99){
0214 h4a->Fill(status);
0215 }
0216 else{
0217 vec.push_back(std::make_tuple(trk_index,status,assoc_weight));
0218 }
0219
0220 }
0221
0222
0223 std::sort(vec.begin(), vec.end());
0224
0225
0226 size_t ivec = 0;
0227 while (ivec < vec.size()) {
0228
0229 int current_first = std::get<0>(vec[ivec]);
0230 float current_max_weight = std::get<2>(vec[ivec]);
0231 size_t start = ivec;
0232 size_t index_max_weight = ivec;
0233
0234
0235 while (ivec < vec.size() && std::get<0>(vec[ivec]) == current_first) {
0236 float current_weight = std::get<2>(vec[ivec]);
0237 if( current_weight > current_max_weight ){
0238 current_max_weight = current_weight;
0239 index_max_weight = ivec;
0240 }
0241 ++ivec;
0242 }
0243
0244
0245 if (ivec - start > 1) {
0246 for (size_t ii = start; ii < ivec; ++ii){
0247
0248 if(ii == index_max_weight)
0249 h5a->Fill( std::get<1>(vec[ii]) );
0250 else
0251 h6a->Fill( std::get<1>(vec[ii]) );
0252 }
0253 }
0254 }
0255
0256
0257 vec.clear();
0258 ivec = 0;
0259
0260
0261
0262 for(int iassoc=0;iassoc<hit_assoc_weight_rs.GetSize();iassoc++){
0263
0264 auto assoc_weight = hit_assoc_weight_rs[iassoc];
0265
0266 h1b->Fill(assoc_weight);
0267
0268
0269 if(assoc_weight > 0.99)
0270 pure_track_counter_rs++;
0271 else if(iassoc==0)
0272 mixed_track_counter_rs++;
0273 else if(iassoc>0 && hit_assoc_track_rs[iassoc] != hit_assoc_track_rs[iassoc-1])
0274 mixed_track_counter_rs++;
0275
0276
0277
0278 int mc_index = hit_assoc_mc_rs[iassoc];
0279 int trk_index = hit_assoc_track_rs[iassoc];
0280
0281 gen_vec.SetXYZM(gen_px[mc_index],gen_py[mc_index],gen_pz[mc_index],gen_mass[mc_index]);
0282 auto gen_mom = gen_vec.P();
0283
0284 auto trk_mom = fabs(1./track_qoverp_rs[trk_index]);
0285
0286 h2b->Fill(gen_mom,trk_mom);
0287
0288 if( fabs(gen_vec.Eta()) < 3 )
0289 h3b->Fill(gen_mom,trk_mom);
0290
0291
0292 int status = gen_status[mc_index];
0293
0294 if(assoc_weight > 0.99){
0295 h4b->Fill(status);
0296 }
0297 else{
0298 vec.push_back(std::make_tuple(trk_index,status,assoc_weight));
0299 }
0300
0301 }
0302
0303
0304 std::sort(vec.begin(), vec.end());
0305
0306
0307 while (ivec < vec.size()) {
0308
0309 int current_first = std::get<0>(vec[ivec]);
0310 float current_max_weight = std::get<2>(vec[ivec]);
0311 size_t start = ivec;
0312 size_t index_max_weight = ivec;
0313
0314
0315 while (ivec < vec.size() && std::get<0>(vec[ivec]) == current_first) {
0316 float current_weight = std::get<2>(vec[ivec]);
0317 if( current_weight > current_max_weight ){
0318 current_max_weight = current_weight;
0319 index_max_weight = ivec;
0320 }
0321 ++ivec;
0322 }
0323
0324
0325 if (ivec - start > 1) {
0326 for (size_t ii = start; ii < ivec; ++ii){
0327
0328 if(ii == index_max_weight)
0329 h5b->Fill( std::get<1>(vec[ii]) );
0330 else
0331 h6b->Fill( std::get<1>(vec[ii]) );
0332 }
0333 }
0334 }
0335
0336
0337 }
0338
0339
0340 float frac_pure_ts = ((float) pure_track_counter_ts) / ((float) (pure_track_counter_ts + mixed_track_counter_ts));
0341 float frac_pure_rs = ((float) pure_track_counter_rs) / ((float) (pure_track_counter_rs + mixed_track_counter_rs));
0342
0343
0344
0345
0346 TCanvas *c1a = new TCanvas("c1a");
0347 c1a->SetLogy();
0348 h1a->Draw();
0349
0350 TLegend *leg1a = new TLegend(0.125,0.7,0.475,0.875);
0351 leg1a->SetBorderSize(0);leg1a->SetFillStyle(0);
0352 leg1a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0353 leg1a->AddEntry(h1a, Form("Fraction of entirely pure tracks = %.3f",frac_pure_ts),"l");
0354 leg1a->Draw();
0355
0356
0357 TCanvas *c1b = new TCanvas("c1b");
0358 c1b->SetLogy();
0359 h1b->Draw();
0360
0361 TLegend *leg1b = new TLegend(0.125,0.7,0.475,0.875);
0362 leg1b->SetBorderSize(0);leg1b->SetFillStyle(0);
0363 leg1b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0364 leg1b->AddEntry(h1b, Form("Fraction of entirely pure tracks = %.3f",frac_pure_rs),"l");
0365 leg1b->Draw();
0366
0367
0368 TCanvas *c2a = new TCanvas("c2a");
0369 c2a->SetLogz();
0370 h2a->Draw("colz");
0371
0372 TLatex *tex2a = new TLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0373 tex2a->SetTextSize(0.03);
0374 tex2a->Draw();
0375
0376
0377 TCanvas *c2b = new TCanvas("c2b");
0378 c2b->SetLogz();
0379 h2b->Draw("colz");
0380
0381 TLatex *tex2b = new TLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0382 tex2b->SetTextSize(0.03);
0383 tex2b->Draw();
0384
0385
0386 TCanvas *c3a = new TCanvas("c3a");
0387 c3a->SetLogz();
0388 h3a->Draw("colz");
0389
0390 TLatex *tex3a = new TLatex();
0391 tex3a->SetTextSize(0.03);
0392 tex3a->DrawLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0393 tex3a->DrawLatex(25,5," -3 < #eta_{gen} < +3");
0394
0395
0396 TCanvas *c3b = new TCanvas("c3b");
0397 c3b->SetLogz();
0398 h3b->Draw("colz");
0399
0400 TLatex *tex3b = new TLatex();
0401 tex3b->SetTextSize(0.03);
0402 tex3b->DrawLatex(25,8,Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0403 tex3b->DrawLatex(25,5," -3 < #eta_{gen} < +3");
0404
0405
0406 TCanvas *c4a = new TCanvas("c4a");
0407 c4a->SetLogy();
0408 h4a->Draw();
0409
0410 TLegend *leg4a = new TLegend(0.125,0.7,0.475,0.875);
0411 leg4a->SetBorderSize(0);leg4a->SetFillStyle(0);
0412 leg4a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0413 leg4a->Draw();
0414
0415
0416 TCanvas *c4b = new TCanvas("c4b");
0417 c4b->SetLogy();
0418 h4b->Draw();
0419
0420 TLegend *leg4b = new TLegend(0.125,0.7,0.475,0.875);
0421 leg4b->SetBorderSize(0);leg4b->SetFillStyle(0);
0422 leg4b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0423 leg4b->Draw();
0424
0425
0426 TCanvas *c5a = new TCanvas("c5a");
0427 c5a->SetLogy();
0428 h5a->Draw();
0429
0430 TLegend *leg5a = new TLegend(0.125,0.7,0.475,0.875);
0431 leg5a->SetBorderSize(0);leg5a->SetFillStyle(0);
0432 leg5a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0433 leg5a->Draw();
0434
0435 TCanvas *c6a = new TCanvas("c6a");
0436 c6a->SetLogy();
0437 h6a->Draw();
0438
0439 TLegend *leg6a = new TLegend(0.525,0.7,0.875,0.875);
0440 leg6a->SetBorderSize(0);leg6a->SetFillStyle(0);
0441 leg6a->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0442 leg6a->Draw();
0443
0444
0445 TCanvas *c5b = new TCanvas("c5b");
0446 c5b->SetLogy();
0447 h5b->Draw();
0448
0449 TLegend *leg5b = new TLegend(0.125,0.7,0.475,0.875);
0450 leg5b->SetBorderSize(0);leg5b->SetFillStyle(0);
0451 leg5b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0452 leg5b->Draw();
0453
0454 TCanvas *c6b = new TCanvas("c6b");
0455 c6b->SetLogy();
0456 h6b->Draw();
0457
0458 TLegend *leg6b = new TLegend(0.525,0.7,0.875,0.875);
0459 leg6b->SetBorderSize(0);leg6b->SetFillStyle(0);
0460 leg6b->SetHeader(Form("Pythia8: %s GeV, Q^{2} > %d GeV^{2}",beam_energies.Data(),Q2_set));
0461 leg6b->Draw();
0462
0463
0464 c1a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf[",beam_energies.Data(),Q2_set));
0465 c1a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0466 c1b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0467 c2a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0468 c2b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0469 c3a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0470 c3b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0471 c4a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0472 c4b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0473 c5a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0474 c6a->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0475 c5b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0476 c6b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf",beam_energies.Data(),Q2_set));
0477 c6b->Print(Form("plots/hit_based_matching_%s_Q2_%d.pdf]",beam_energies.Data(),Q2_set));
0478
0479 }