From 0c5110a177d36a48140be4f9adb6bddddb6d339f Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Tue, 12 Sep 2023 10:22:21 -0500 Subject: [PATCH 1/3] Updating demmcsim so that trkprimary is always first --- src/InfoMCStructHelper.cc | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/InfoMCStructHelper.cc b/src/InfoMCStructHelper.cc index de7d32e..118eadb 100644 --- a/src/InfoMCStructHelper.cc +++ b/src/InfoMCStructHelper.cc @@ -129,17 +129,6 @@ namespace mu2e { auto trkprimaryptr = kseedmc.simParticle().simParticle(_spcH); auto trkprimary = trkprimaryptr->originParticle(); - // Add all the primary particles first - SimInfo sim_info; - for(auto const& spp : primary.primarySimParticles()){ - // check whether we already put this primary in - fillSimInfo(spp, sim_info); - sim_info.trkrel = MCRelationship(spp, trkprimaryptr); - sim_info.prirel = MCRelationship(spp, spp); - siminfos.push_back(sim_info); - } - - auto current_sim_particle_ptr = trkprimaryptr; auto current_sim_particle = trkprimary; if (n_generations == -1) { // means do all generations @@ -162,10 +151,7 @@ namespace mu2e { } sim_info.prirel = bestrel; - // We already added all the primaries - if (sim_info.prirel != MCRelationship::same) { - siminfos.push_back(sim_info); - } + siminfos.push_back(sim_info); if (current_sim_particle.parent().isNonnull()) { current_sim_particle_ptr = current_sim_particle.parent(); current_sim_particle = current_sim_particle_ptr->originParticle(); @@ -174,6 +160,19 @@ namespace mu2e { break; // this particle doesn't have a parent } } + + // Now add all the primary particles + SimInfo sim_info; + for(auto const& spp : primary.primarySimParticles()){ + // check whether we already put this primary in + fillSimInfo(spp, sim_info); + sim_info.trkrel = MCRelationship(spp, trkprimaryptr); + sim_info.prirel = MCRelationship(spp, spp); + if (sim_info.trkrel != MCRelationship::same) { + siminfos.push_back(sim_info); + } + } + } From 85d437271f1008f2a1ee38a09bca83626d25f031 Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Tue, 12 Sep 2023 10:36:16 -0500 Subject: [PATCH 2/3] Updating example genealogy scripts --- example-analysis-scripts/Genealogy.C | 40 +++++++++++++-------------- example-analysis-scripts/Genealogy.py | 36 ++++++++++++------------ 2 files changed, 38 insertions(+), 38 deletions(-) diff --git a/example-analysis-scripts/Genealogy.C b/example-analysis-scripts/Genealogy.C index 5e7071c..0944515 100644 --- a/example-analysis-scripts/Genealogy.C +++ b/example-analysis-scripts/Genealogy.C @@ -11,8 +11,8 @@ void Genealogy() { c1->Divide(2,3); c1->cd(1); - trkana->Draw("demmcsim.mom.R()>>hMomGen0(100,100,110)", "demmcsim.generation==0", "PLC"); // "PLC" means new colors will be picked automatically - trkana->Draw("demmcsim.mom.R()>>hMomGen0_Ce(100,100,110)", "demmcsim.generation==0 && demmcsim.gen==167", "PFC SAME"); + trkana->Draw("demmcsim.mom.R()>>hMomGen0(100,100,110)", "demmcsim.trkrel._rem==0", "PLC"); // "PLC" means new colors will be picked automatically + trkana->Draw("demmcsim.mom.R()>>hMomGen0_Ce(100,100,110)", "demmcsim.trkrel._rem==0 && demmcsim.gen==167", "PFC SAME"); TH1D* hMomGen0 = ((TH1D*)c1->GetPad(1)->GetPrimitive("hMomGen0")); hMomGen0->SetTitle(""); TH1D* hMomGen0_Ce = ((TH1D*)c1->GetPad(1)->GetPrimitive("hMomGen0_Ce")); @@ -20,13 +20,13 @@ void Genealogy() { leg_1->SetBorderSize(0); leg_1->SetTextSize(0.035); leg_1->SetFillColor(kWhite); - leg_1->AddEntry(hMomGen0, "generation = -1", "l"); - leg_1->AddEntry(hMomGen0_Ce, "generation = -1 and proc=CeEndpoint", "f"); + leg_1->AddEntry(hMomGen0, "generation = 0", "l"); + leg_1->AddEntry(hMomGen0_Ce, "generation = 0 and proc=CeEndpoint", "f"); leg_1->Draw(); c1->cd(2); - trkana->Draw("demmcsim.mom.R()>>hMomGen_min1(100,0,100)", "demmcsim.generation==-1", "PLC"); // "PLC" means new colors will be picked automatically - trkana->Draw("demmcsim.mom.R()>>hMomGen_min1_mu(100,0,100)", "demmcsim.generation==-1 && demmcsim.pdg==13", "PFC SAME"); + trkana->Draw("demmcsim.mom.R()>>hMomGen_min1(100,0,100)", "demmcsim.trkrel._rem==1", "PLC"); // "PLC" means new colors will be picked automatically + trkana->Draw("demmcsim.mom.R()>>hMomGen_min1_mu(100,0,100)", "demmcsim.trkrel._rem==1 && demmcsim.pdg==13", "PFC SAME"); TH1D* hMomGen_min1 = ((TH1D*)c1->GetPad(2)->GetPrimitive("hMomGen_min1")); hMomGen_min1->SetTitle(""); TH1D* hMomGen_min1_mu = ((TH1D*)c1->GetPad(2)->GetPrimitive("hMomGen_min1_mu")); @@ -34,13 +34,13 @@ void Genealogy() { leg_2->SetBorderSize(0); leg_2->SetTextSize(0.035); leg_2->SetFillColor(kWhite); - leg_2->AddEntry(hMomGen_min1, "generation = -1", "l"); - leg_2->AddEntry(hMomGen_min1_mu, "generation = -1 and pdg=mu-", "f"); + leg_2->AddEntry(hMomGen_min1, "generation = 1", "l"); + leg_2->AddEntry(hMomGen_min1_mu, "generation = 1 and pdg=mu-", "f"); leg_2->Draw(); c1->cd(3); - trkana->Draw("demmcsim.mom.R()>>hMomGen_min2(100,0,200)", "demmcsim.generation==-2", "PLC"); // "PLC" means new colors will be picked automatically - trkana->Draw("demmcsim.mom.R()>>hMomGen_min2_pi(100,0,200)", "demmcsim.generation==-2 && demmcsim.pdg==-211", "PFC SAME"); + trkana->Draw("demmcsim.mom.R()>>hMomGen_min2(100,0,200)", "demmcsim.trkrel._rem==2", "PLC"); // "PLC" means new colors will be picked automatically + trkana->Draw("demmcsim.mom.R()>>hMomGen_min2_pi(100,0,200)", "demmcsim.trkrel._rem==2 && demmcsim.pdg==-211", "PFC SAME"); TH1D* hMomGen_min2 = ((TH1D*)c1->GetPad(3)->GetPrimitive("hMomGen_min2")); hMomGen_min2->SetTitle(""); TH1D* hMomGen_min2_pi = ((TH1D*)c1->GetPad(3)->GetPrimitive("hMomGen_min2_pi")); @@ -48,14 +48,14 @@ void Genealogy() { leg_3->SetBorderSize(0); leg_3->SetTextSize(0.035); leg_3->SetFillColor(kWhite); - leg_3->AddEntry(hMomGen_min2, "generation = -2", "l"); - leg_3->AddEntry(hMomGen_min2_pi, "generation = -2 and pdg=pi-", "f"); + leg_3->AddEntry(hMomGen_min2, "generation = 2", "l"); + leg_3->AddEntry(hMomGen_min2_pi, "generation = 2 and pdg=pi-", "f"); leg_3->Draw(); c1->cd(4); - trkana->Draw("demmcsim.mom.R()>>hMomGen_min3(100,0,9000)", "demmcsim.generation==-3", "PLC"); // "PLC" means new colors will be picked automatically - trkana->Draw("demmcsim.mom.R()>>hMomGen_min3_pot(100,0,9000)", "demmcsim.generation==-3 && demmcsim.gen==56", "PFC SAME"); + trkana->Draw("demmcsim.mom.R()>>hMomGen_min3(100,0,9000)", "demmcsim.trkrel._rem==3", "PLC"); // "PLC" means new colors will be picked automatically + trkana->Draw("demmcsim.mom.R()>>hMomGen_min3_pot(100,0,9000)", "demmcsim.trkrel._rem==3 && demmcsim.gen==56", "PFC SAME"); TH1D* hMomGen_min3 = ((TH1D*)c1->GetPad(4)->GetPrimitive("hMomGen_min3")); hMomGen_min3->SetTitle(""); TH1D* hMomGen_min3_pot = ((TH1D*)c1->GetPad(4)->GetPrimitive("hMomGen_min3_pot")); @@ -63,13 +63,13 @@ void Genealogy() { leg_4->SetBorderSize(0); leg_4->SetTextSize(0.035); leg_4->SetFillColor(kWhite); - leg_4->AddEntry(hMomGen_min3, "generation = -3", "l"); - leg_4->AddEntry(hMomGen_min3_pot, "generation = -3 and proc=POT", "f"); + leg_4->AddEntry(hMomGen_min3, "generation = 3", "l"); + leg_4->AddEntry(hMomGen_min3_pot, "generation = 3 and proc=POT", "f"); leg_4->Draw(); c1->cd(5); - trkana->Draw("demmcsim.mom.R()>>hMomGen_min4(100,0,9000)", "demmcsim.generation==-4", "PLC"); // "PLC" means new colors will be picked automatically - trkana->Draw("demmcsim.mom.R()>>hMomGen_min4_pot(100,0,9000)", "demmcsim.generation==-4 && demmcsim.gen==56", "PFC SAME"); + trkana->Draw("demmcsim.mom.R()>>hMomGen_min4(100,0,9000)", "demmcsim.trkrel._rem==4", "PLC"); // "PLC" means new colors will be picked automatically + trkana->Draw("demmcsim.mom.R()>>hMomGen_min4_pot(100,0,9000)", "demmcsim.trkrel._rem==4 && demmcsim.gen==56", "PFC SAME"); TH1D* hMomGen_min4 = ((TH1D*)c1->GetPad(5)->GetPrimitive("hMomGen_min4")); hMomGen_min4->SetTitle(""); TH1D* hMomGen_min4_pot = ((TH1D*)c1->GetPad(5)->GetPrimitive("hMomGen_min4_pot")); @@ -77,8 +77,8 @@ void Genealogy() { leg_5->SetBorderSize(0); leg_5->SetTextSize(0.035); leg_5->SetFillColor(kWhite); - leg_5->AddEntry(hMomGen_min4, "generation = -4", "l"); - leg_5->AddEntry(hMomGen_min4_pot, "generation = -4 and proc=POT", "f"); + leg_5->AddEntry(hMomGen_min4, "generation = 4", "l"); + leg_5->AddEntry(hMomGen_min4_pot, "generation = 4 and proc=POT", "f"); leg_5->Draw(); c1->cd(6); diff --git a/example-analysis-scripts/Genealogy.py b/example-analysis-scripts/Genealogy.py index 587f677..3a6550f 100644 --- a/example-analysis-scripts/Genealogy.py +++ b/example-analysis-scripts/Genealogy.py @@ -16,7 +16,7 @@ def add_useful_columns(batch): # will calculate some useful values batch['demmcsim.mom'] = (batch['demmcsim.mom.fCoordinates.fX']**2 + batch['demmcsim.mom.fCoordinates.fY']**2 + batch['demmcsim.mom.fCoordinates.fZ']**2)**0.5 def zip_branches(batch): # will return a reduced awkward array with the branches we want - cols = ['demmcsim.mom', 'demmcsim.pdg', 'demmcsim.generation', 'demmcsim.gen'] + cols = ['demmcsim.mom', 'demmcsim.pdg', 'demmcsim.trkrel._rem', 'demmcsim.gen'] return ak.zip({col : batch[col] for col in cols}) # The script will start here @@ -35,11 +35,11 @@ def zip_branches(batch): # will return a reduced awkward array with the branches # Now define some useful masks # - for the different generations -gen_0_mask = reduced_trkana['demmcsim.generation']==0 -gen_min1_mask = reduced_trkana['demmcsim.generation']==-1 -gen_min2_mask = reduced_trkana['demmcsim.generation']==-2 -gen_min3_mask = reduced_trkana['demmcsim.generation']==-3 -gen_min4_mask = reduced_trkana['demmcsim.generation']==-4 +gen_0_mask = reduced_trkana['demmcsim.trkrel._rem']==0 +gen_min1_mask = reduced_trkana['demmcsim.trkrel._rem']==1 +gen_min2_mask = reduced_trkana['demmcsim.trkrel._rem']==2 +gen_min3_mask = reduced_trkana['demmcsim.trkrel._rem']==3 +gen_min4_mask = reduced_trkana['demmcsim.trkrel._rem']==4 # - for different particle creation processes process = MCDataProducts.ProcessCode() proc_ce_mask = reduced_trkana['demmcsim.gen']==process.mu2eCeMinusEndpoint @@ -52,16 +52,16 @@ def zip_branches(batch): # will return a reduced awkward array with the branches # - we will show that everything at gen=0 is a Ce gen_0_mom=ak.flatten(reduced_trkana['demmcsim.mom'][(gen_0_mask)]).to_numpy() gen_0_ce_mom=ak.flatten(batch['demmcsim.mom'][(gen_0_mask) & (proc_ce_mask)]).to_numpy() -# - we will show that everything at gen=-1 is a muon +# - we will show that everything at gen=1 is a muon gen_min1_mom=ak.flatten(reduced_trkana['demmcsim.mom'][(gen_min1_mask)]).to_numpy() gen_min1_mu_mom=ak.flatten(reduced_trkana['demmcsim.mom'][(gen_min1_mask) & (pdg_mu_mask)]).to_numpy() -# - we will show that everything at gen=-2 is a pion +# - we will show that everything at gen=2 is a pion gen_min2_mom=ak.flatten(reduced_trkana['demmcsim.mom'][(gen_min2_mask)]).to_numpy() gen_min2_pi_mom=ak.flatten(reduced_trkana['demmcsim.mom'][(gen_min2_mask) & (pdg_pi_mask)]).to_numpy() -# - we will show that *not* everything at gen=-3 is a POT +# - we will show that *not* everything at gen=3 is a POT gen_min3_mom=ak.flatten(reduced_trkana['demmcsim.mom'][(gen_min3_mask)]).to_numpy() gen_min3_pot_mom=ak.flatten(reduced_trkana['demmcsim.mom'][(gen_min3_mask) & (proc_pot_mask)]).to_numpy() -# - likewise, we will show that *not* everything at gen=-4 is a POT +# - likewise, we will show that *not* everything at gen=4 is a POT gen_min4_mom=ak.flatten(reduced_trkana['demmcsim.mom'][(gen_min4_mask)]).to_numpy() gen_min4_pot_mom=ak.flatten(reduced_trkana['demmcsim.mom'][(gen_min4_mask) & (proc_pot_mask)]).to_numpy() # - finally, we will show all the POTs regardless of generation @@ -72,17 +72,17 @@ def zip_branches(batch): # will return a reduced awkward array with the branches bins, counts, patches = axs.ravel()[0].hist(gen_0_mom, bins=n_mom_bins, range=(100,110), log=False, histtype='step', label='generation=0 (count = {})'.format(len(gen_0_mom))) bins, counts, patches = axs.ravel()[0].hist(gen_0_ce_mom, bins=n_mom_bins, range=(100,110), log=False, histtype='stepfilled', label='generation=0 and proc=mu2eCeMinusEndpoint (count = {})'.format(len(gen_0_ce_mom))) -bins, counts, patches = axs.ravel()[1].hist(gen_min1_mom, bins=n_mom_bins, range=(0,100), log=False, histtype='step', label='generation=-1 (count = {})'.format(len(gen_min1_mom))) -bins, counts, patches = axs.ravel()[1].hist(gen_min1_mu_mom, bins=n_mom_bins, range=(0,100), log=False, histtype='stepfilled', label='generation=-1 and pdg==mu- (count = {})'.format(len(gen_min1_mu_mom))) +bins, counts, patches = axs.ravel()[1].hist(gen_min1_mom, bins=n_mom_bins, range=(0,100), log=False, histtype='step', label='generation=1 (count = {})'.format(len(gen_min1_mom))) +bins, counts, patches = axs.ravel()[1].hist(gen_min1_mu_mom, bins=n_mom_bins, range=(0,100), log=False, histtype='stepfilled', label='generation=1 and pdg==mu- (count = {})'.format(len(gen_min1_mu_mom))) -bins, counts, patches = axs.ravel()[2].hist(gen_min2_mom, bins=n_mom_bins, range=(0,200), log=False, histtype='step', label='generation=-2 (count = {})'.format(len(gen_min2_mom))) -bins, counts, patches = axs.ravel()[2].hist(gen_min2_pi_mom, bins=n_mom_bins, range=(0,200), log=False, histtype='stepfilled', label='generation=-2 and pdg==pi- (count = {})'.format(len(gen_min2_pi_mom))) +bins, counts, patches = axs.ravel()[2].hist(gen_min2_mom, bins=n_mom_bins, range=(0,200), log=False, histtype='step', label='generation=2 (count = {})'.format(len(gen_min2_mom))) +bins, counts, patches = axs.ravel()[2].hist(gen_min2_pi_mom, bins=n_mom_bins, range=(0,200), log=False, histtype='stepfilled', label='generation=2 and pdg==pi- (count = {})'.format(len(gen_min2_pi_mom))) -bins, counts, patches = axs.ravel()[3].hist(gen_min3_mom, bins=n_mom_bins, range=(0, 9000), log=False, histtype='step', label='generation=-3 (count = {})'.format(len(gen_min3_mom))) -bins, counts, patches = axs.ravel()[3].hist(gen_min3_pot_mom, bins=n_mom_bins, range=(0, 9000), log=False, histtype='stepfilled', label='generation=-3 and proc=POT (count = {})'.format(len(gen_min3_pot_mom))) +bins, counts, patches = axs.ravel()[3].hist(gen_min3_mom, bins=n_mom_bins, range=(0, 9000), log=False, histtype='step', label='generation=3 (count = {})'.format(len(gen_min3_mom))) +bins, counts, patches = axs.ravel()[3].hist(gen_min3_pot_mom, bins=n_mom_bins, range=(0, 9000), log=False, histtype='stepfilled', label='generation=3 and proc=POT (count = {})'.format(len(gen_min3_pot_mom))) -bins, counts, patches = axs.ravel()[4].hist(gen_min4_mom, bins=n_mom_bins, range=(0, 9000), log=False, histtype='step', label='generation=-4 (count = {})'.format(len(gen_min4_mom))) -bins, counts, patches = axs.ravel()[4].hist(gen_min4_pot_mom, bins=n_mom_bins, range=(0, 9000), log=False, histtype='stepfilled', label='generation=-4 and proc=POT (count = {})'.format(len(gen_min4_pot_mom))) +bins, counts, patches = axs.ravel()[4].hist(gen_min4_mom, bins=n_mom_bins, range=(0, 9000), log=False, histtype='step', label='generation=4 (count = {})'.format(len(gen_min4_mom))) +bins, counts, patches = axs.ravel()[4].hist(gen_min4_pot_mom, bins=n_mom_bins, range=(0, 9000), log=False, histtype='stepfilled', label='generation=4 and proc=POT (count = {})'.format(len(gen_min4_pot_mom))) bins, counts, patches = axs.ravel()[5].hist(gen_all_pot_mom, bins=n_mom_bins, range=(0, 9000), log=False, histtype='stepfilled', label='all generations, proc=POT (count = {})'.format(len(gen_all_pot_mom))) for axi in axs.ravel(): From 478e1354f95fcf713ce5644b69064665e8bab478 Mon Sep 17 00:00:00 2001 From: Andrew Edmonds Date: Wed, 13 Sep 2023 12:30:10 -0500 Subject: [PATCH 3/3] Fix two issues: one where the primary particle was duplicated, the other where the particles were passed to MCRelationship in the wrong order --- src/InfoMCStructHelper.cc | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/InfoMCStructHelper.cc b/src/InfoMCStructHelper.cc index 118eadb..fb38ebb 100644 --- a/src/InfoMCStructHelper.cc +++ b/src/InfoMCStructHelper.cc @@ -138,12 +138,12 @@ namespace mu2e { for (int i_generation = 0; i_generation < n_generations; ++i_generation) { SimInfo sim_info; fillSimInfo(current_sim_particle, sim_info); - sim_info.trkrel = MCRelationship(trkprimaryptr, current_sim_particle_ptr); + sim_info.trkrel = MCRelationship(current_sim_particle_ptr, trkprimaryptr); auto bestprimarysp = primary.primarySimParticles().front(); MCRelationship bestrel; for(auto const& spp : primary.primarySimParticles()){ - MCRelationship mcrel(spp,current_sim_particle_ptr); + MCRelationship mcrel(current_sim_particle_ptr, spp); if(mcrel > bestrel){ bestrel = mcrel; bestprimarysp = spp; @@ -164,11 +164,19 @@ namespace mu2e { // Now add all the primary particles SimInfo sim_info; for(auto const& spp : primary.primarySimParticles()){ - // check whether we already put this primary in fillSimInfo(spp, sim_info); - sim_info.trkrel = MCRelationship(spp, trkprimaryptr); - sim_info.prirel = MCRelationship(spp, spp); - if (sim_info.trkrel != MCRelationship::same) { + + // check whether we already put this primary in + bool already_added = false; + for (const auto& i_sim_info : siminfos) { + if (i_sim_info.prirel == MCRelationship::same) { + already_added = true; + break; + } + } + if (!already_added) { + sim_info.trkrel = MCRelationship(spp, trkprimaryptr); + sim_info.prirel = MCRelationship(spp, spp); siminfos.push_back(sim_info); } }