Skip to content

Commit

Permalink
Merge pull request #101 from AndrewEdmonds11/genealogy-update
Browse files Browse the repository at this point in the history 
Genealogy update
  • Loading branch information
@AndrewEdmonds11
AndrewEdmonds11 authored Sep 13, 2023
2 parents 474e5dd + 478e135 commit b8acbea
Show file tree
Hide file tree
Showing 3 changed files with 62 additions and 55 deletions.
40 changes: 20 additions & 20 deletions example-analysis-scripts/Genealogy.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,74 +11,74 @@ 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"));
TLegend* leg_1 = new TLegend(0.2, 0.5, 0.4, 0.8); // defined in NDC coordinates (i.e. fraction of canvas)
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"));
TLegend* leg_2 = new TLegend(0.2, 0.5, 0.4, 0.8); // defined in NDC coordinates (i.e. fraction of canvas)
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"));
TLegend* leg_3 = new TLegend(0.2, 0.5, 0.4, 0.8); // defined in NDC coordinates (i.e. fraction of canvas)
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"));
TLegend* leg_4 = new TLegend(0.2, 0.5, 0.4, 0.8); // defined in NDC coordinates (i.e. fraction of canvas)
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"));
TLegend* leg_5 = new TLegend(0.2, 0.5, 0.4, 0.8); // defined in NDC coordinates (i.e. fraction of canvas)
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);
Expand Down
36 changes: 18 additions & 18 deletions example-analysis-scripts/Genealogy.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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():
Expand Down
41 changes: 24 additions & 17 deletions src/InfoMCStructHelper.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -149,23 +138,20 @@ 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;
}
}
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();
Expand All @@ -174,6 +160,27 @@ 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()){
fillSimInfo(spp, sim_info);

// 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);
}
}

}


Expand Down

0 comments on commit b8acbea

Please sign in to comment.