Skip to content
Snippets Groups Projects
Commit 4258974f authored by Eric Torrence's avatar Eric Torrence
Browse files

Merge branch 'ntupledumper' into 'master' 

update ntuple: blinding default to 100 GeV, add truth positions to DP daughters, add raw calibrated energy

See merge request faser/calypso!340
parents fad536c9 b5863048
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+ 110
− 108
View file @ 4258974f
......@@ -170,7 +170,9 @@ StatusCode NtupleDumperAlg::initialize()
addBranch("Calo_total_nMIP", &m_calo_total_nMIP);
addBranch("Calo_total_E_dep", &m_calo_total_E_dep);
addBranch("Calo_total_E_EM", &m_calo_total_E_EM);
addBranch("Calo_total_raw_E_EM", &m_calo_total_raw_E_EM);
addBranch("Calo_total_E_EM_fudged", &m_calo_total_E_EM_fudged);
addBranch("Calo_total_raw_E_EM_fudged", &m_calo_total_raw_E_EM_fudged);
addCalibratedBranches("Preshower",2,12);
addBranch("Preshower_total_nMIP", &m_preshower_total_nMIP);
......@@ -582,119 +584,114 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
if (isMC) { // if simulation find MC cross section and primary lepton
// Work out effective cross section for MC
if (m_useFlukaWeights)
{
double flukaWeight = truthEventContainer->at(0)->weights()[0];
ATH_MSG_ALWAYS("Found fluka weight = " << flukaWeight);
m_crossSection = m_baseEventCrossSection * flukaWeight;
}
else if (m_useGenieWeights)
{
m_crossSection = m_baseEventCrossSection;
}
else
{
if (m_useFlukaWeights) {
double flukaWeight = truthEventContainer->at(0)->weights()[0];
ATH_MSG_ALWAYS("Found fluka weight = " << flukaWeight);
m_crossSection = m_baseEventCrossSection * flukaWeight;
} else if (m_useGenieWeights) {
m_crossSection = m_baseEventCrossSection;
} else {
//ATH_MSG_WARNING("Monte carlo event with no weighting scheme specified. Setting crossSection (weight) to " << m_baseEventCrossSection << " fb.");
m_crossSection = m_baseEventCrossSection;
m_crossSection = m_baseEventCrossSection;
}
// Find the M d0 and d1 truth information
// Find truth particle information
SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx };
if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0)
{
if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0) {
int ipart(0);
for (auto particle : *truthParticleContainer)
{
// loop over first 10 truth particles (for non A' samples)
if (ipart++ < 10) {
m_truth_P.push_back(particle->p4().P());
m_truth_px.push_back(particle->p4().X());
m_truth_py.push_back(particle->p4().Y());
m_truth_pz.push_back(particle->p4().Z());
m_truth_m.push_back(particle->m());
m_truth_pdg.push_back(particle->pdgId());
if ( particle->hasProdVtx()) {
m_truth_prod_x.push_back(particle->prodVtx()->x());
m_truth_prod_y.push_back(particle->prodVtx()->y());
m_truth_prod_z.push_back(particle->prodVtx()->z());
} else {
m_truth_prod_x.push_back(999999);
m_truth_prod_y.push_back(999999);
m_truth_prod_z.push_back(999999);
}
if ( particle->hasDecayVtx()) {
m_truth_dec_x.push_back(particle->decayVtx()->x());
m_truth_dec_y.push_back(particle->decayVtx()->y());
m_truth_dec_z.push_back(particle->decayVtx()->z());
} else {
m_truth_dec_x.push_back(999999);
m_truth_dec_y.push_back(999999);
m_truth_dec_z.push_back(999999);
}
}
if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 ) {
if ( particle->pdgId() == 32) { // mother particle (A')
m_truthM_P.push_back(particle->p4().P());
m_truthM_px.push_back(particle->p4().X());
m_truthM_py.push_back(particle->p4().Y());
m_truthM_pz.push_back(particle->p4().Z());
if ( particle->hasDecayVtx()) { // decay vertex for A' particle
m_truthM_x.push_back(particle->decayVtx()->x());
m_truthM_y.push_back(particle->decayVtx()->y());
m_truthM_z.push_back(particle->decayVtx()->z());
} else {
m_truthM_x.push_back(999999);
m_truthM_y.push_back(999999);
m_truthM_z.push_back(999999);
}
}
if ( particle->pdgId() == 11) // daughter particle (positron)
{
m_truthd0_P.push_back(particle->p4().P());
m_truthd0_px.push_back(particle->p4().X());
m_truthd0_py.push_back(particle->p4().Y());
m_truthd0_pz.push_back(particle->p4().Z());
if ( particle->hasProdVtx()) {
m_truthd0_x.push_back(particle->prodVtx()->x());
m_truthd0_y.push_back(particle->prodVtx()->y());
m_truthd0_z.push_back(particle->prodVtx()->z());
} else {
m_truthd0_x.push_back(999999);
m_truthd0_y.push_back(999999);
m_truthd0_z.push_back(999999);
}
}
if ( particle->pdgId() == -11) // daughter particle (electron)
{
m_truthd1_P.push_back(particle->p4().P());
m_truthd1_px.push_back(particle->p4().X());
m_truthd1_py.push_back(particle->p4().Y());
m_truthd1_pz.push_back(particle->p4().Z());
if ( particle->hasProdVtx()) {
m_truthd1_x.push_back(particle->prodVtx()->x());
m_truthd1_y.push_back(particle->prodVtx()->y());
m_truthd1_z.push_back(particle->prodVtx()->z());
} else {
m_truthd1_x.push_back(999999);
m_truthd1_y.push_back(999999);
m_truthd1_z.push_back(999999);
}
}
}
for (auto particle : *truthParticleContainer) { // loop over first 10 truth particles (for non A' samples)
if (ipart++ > 9) break;
m_truth_P.push_back(particle->p4().P());
m_truth_px.push_back(particle->p4().X());
m_truth_py.push_back(particle->p4().Y());
m_truth_pz.push_back(particle->p4().Z());
m_truth_m.push_back(particle->m());
m_truth_pdg.push_back(particle->pdgId());
if ( particle->hasProdVtx()) {
m_truth_prod_x.push_back(particle->prodVtx()->x());
m_truth_prod_y.push_back(particle->prodVtx()->y());
m_truth_prod_z.push_back(particle->prodVtx()->z());
} else {
m_truth_prod_x.push_back(999999);
m_truth_prod_y.push_back(999999);
m_truth_prod_z.push_back(999999);
}
if ( particle->hasDecayVtx()) {
m_truth_dec_x.push_back(particle->decayVtx()->x());
m_truth_dec_y.push_back(particle->decayVtx()->y());
m_truth_dec_z.push_back(particle->decayVtx()->z());
} else {
m_truth_dec_x.push_back(999999);
m_truth_dec_y.push_back(999999);
m_truth_dec_z.push_back(999999);
}
// Find the M d0 and d1 truth information for dark photon
if ( particle->barcode() == 1 || particle->barcode() == 2 || particle->barcode() == 3 ) {
auto positions = m_fiducialParticleTool->getTruthPositions(particle->barcode());
if ( particle->pdgId() == 32) { // mother particle (A')
m_truthM_P.push_back(particle->p4().P());
m_truthM_px.push_back(particle->p4().X());
m_truthM_py.push_back(particle->p4().Y());
m_truthM_pz.push_back(particle->p4().Z());
if ( particle->hasDecayVtx()) { // decay vertex for A' particle
m_truthM_x.push_back(particle->decayVtx()->x());
m_truthM_y.push_back(particle->decayVtx()->y());
m_truthM_z.push_back(particle->decayVtx()->z());
} else {
m_truthM_x.push_back(999999);
m_truthM_y.push_back(999999);
m_truthM_z.push_back(999999);
}
}
if ( particle->pdgId() == 11) { // daughter particle (electron)
m_truthd0_P.push_back(particle->p4().P());
m_truthd0_px.push_back(particle->p4().X());
m_truthd0_py.push_back(particle->p4().Y());
m_truthd0_pz.push_back(particle->p4().Z());
if ( particle->hasProdVtx()) {
m_truthd0_x.push_back(particle->prodVtx()->x());
m_truthd0_y.push_back(particle->prodVtx()->y());
m_truthd0_z.push_back(particle->prodVtx()->z());
} else {
m_truthd0_x.push_back(999999);
m_truthd0_y.push_back(999999);
m_truthd0_z.push_back(999999);
}
for (int station = 1; station < 4; ++station) {
m_truthd0_x.push_back(positions[station].x());
m_truthd0_y.push_back(positions[station].y());
m_truthd0_z.push_back(positions[station].z());
}
}
if ( particle->pdgId() == -11) { // daughter particle (positron)
m_truthd1_P.push_back(particle->p4().P());
m_truthd1_px.push_back(particle->p4().X());
m_truthd1_py.push_back(particle->p4().Y());
m_truthd1_pz.push_back(particle->p4().Z());
if ( particle->hasProdVtx()) {
m_truthd1_x.push_back(particle->prodVtx()->x());
m_truthd1_y.push_back(particle->prodVtx()->y());
m_truthd1_z.push_back(particle->prodVtx()->z());
} else {
m_truthd1_x.push_back(999999);
m_truthd1_y.push_back(999999);
m_truthd1_z.push_back(999999);
}
for (int station = 1; station < 4; ++station) {
m_truthd1_x.push_back(positions[station].x());
m_truthd1_y.push_back(positions[station].y());
m_truthd1_z.push_back(positions[station].z());
}
}
}
}
}
}
......@@ -711,6 +708,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_calo_total_nMIP += hit->Nmip();
m_calo_total_E_dep += hit->E_dep();
m_calo_total_E_EM += hit->E_EM();
m_calo_total_raw_E_EM += (hit->E_EM()*hit->fit_to_raw_ratio());
ATH_MSG_DEBUG("Calibrated calo: ch is " << ch << ", edep is " << hit->E_dep() << ", E_EM is " << hit->E_EM() << ", Nmip is " << hit->Nmip() << ", fit_to_raw_ratio is " << hit->fit_to_raw_ratio());
......@@ -728,8 +726,10 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// add a fudged energy variable that corrects the MC to agree with the testbeam results
m_calo_total_E_EM_fudged = m_calo_total_E_EM;
m_calo_total_raw_E_EM_fudged = m_calo_total_raw_E_EM;
if (isMC) {
m_calo_total_E_EM_fudged = m_calo_total_E_EM_fudged * m_caloMC_FudgeFactor;
m_calo_total_raw_E_EM_fudged = m_calo_total_raw_E_EM_fudged * m_caloMC_FudgeFactor;
}
// load in calibrated preshower container
......@@ -1211,7 +1211,9 @@ NtupleDumperAlg::clearTree() const
m_calo_total_nMIP=0;
m_calo_total_E_dep=0;
m_calo_total_E_EM=0;
m_calo_total_raw_E_EM=0;
m_calo_total_E_EM_fudged=0;
m_calo_total_raw_E_EM_fudged=0;
m_preshower_total_nMIP=0;
m_preshower_total_E_dep=0;
......
PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.h
+ 3
− 1
View file @ 4258974f
......@@ -117,7 +117,7 @@ private:
DoubleProperty m_genieLuminosity { this, "GenieLuminosity", 150.0, "Genie luminosity in inverse fb." };
DoubleProperty m_minMomentum { this, "MinMomentum", 50000.0, "Write out all truth particles with a momentum larger MinMomentum"};
DoubleProperty m_blindingCaloThreshold {this, "BlindingCaloThreshold", 25000.0, "Blind events with Ecal energy above threshold (in MeV)"};
DoubleProperty m_blindingCaloThreshold {this, "BlindingCaloThreshold", 100000.0, "Blind events with Ecal energy above threshold (in MeV)"};
DoubleProperty m_caloMC_FudgeFactor {this, "CaloFudgeFactorMC", 1.088, "Correct the MC energy calibration by this fudge factor"};
double m_baseEventCrossSection {1.0};
......@@ -167,8 +167,10 @@ private:
mutable float m_calo_total_nMIP;
mutable float m_calo_total_E_dep;
mutable float m_calo_total_E_EM;
mutable float m_calo_total_raw_E_EM;
mutable float m_calo_total_E_EM_fudged;
mutable float m_calo_total_raw_E_EM_fudged;
mutable float m_preshower_total_nMIP;
mutable float m_preshower_total_E_dep;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment