Skip to content
Snippets Groups Projects

adding codes to look at seeding performance

Closed adding codes to look at seeding performance
sazhang/calypso:alma9-tracking into alma9-dev
Closed Saimeng Zhang requested to merge sazhang/calypso:alma9-tracking into alma9-dev
Viewing commit 1fd2e445
Prev
Show latest version
14 files
+ 871
− 123
Compare changes
  • Side-by-side
  • Inline
Files
14
PhysicsAnalysis/NtupleDumper/src/NtupleDumperAlg.cxx
+ 99
− 96
@@ -47,32 +47,25 @@ void NtupleDumperAlg::addBranch(const std::string &name,
m_tree->Branch(name.c_str(),var,(name+"/I").c_str());
}
std::map<std::string, std::list<int>>&
NtupleDumperAlg::getChannelMap() const {
// Have to declare this cost to call during execute
void NtupleDumperAlg::defineWaveBranches() const {
ATH_MSG_DEBUG("defineWaveBranches called");
// Use waveform map to find all defined waveform channels
auto mapping = m_mappingTool->getCableMapping();
ATH_MSG_DEBUG("Cable mapping contains " << mapping.size() << " entries");
// mapping is std::map<int, std::pair<std::string, Identifier> >
// use this to fill my own map of channel lists keyed by type
static std::map<std::string, std::list<int>> wave_map;
std::map<std::string, std::list<int>> wave_map;
for (const auto& [key, value] : mapping) {
wave_map[value.first].push_back(key);
ATH_MSG_DEBUG("Found mapping " << value.first << " chan " << key);
}
return wave_map;
}
// Have to declare this cost to call during execute
void NtupleDumperAlg::defineWaveBranches() const {
ATH_MSG_DEBUG("defineWaveBranches called");
auto wave_map = getChannelMap();
// Now go through found types and define ntuple entries
// Keys are defined in cable map and used by RawWaveformDecoder
for (const auto& [key, value] : wave_map) {
@@ -294,7 +287,19 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("TBP", &m_tbp, "TBP/I");
m_tree->Branch("TAP", &m_tap, "TAP/I");
m_tree->Branch("inputBits", &m_inputBits, "inputBits/I");
m_tree->Branch("inputBitsNext", &m_inputBitsNext, "inputBitsNext/I");
m_tree->Branch("inputBitsNext", &m_inputBitsNext, "inputBitsNext/I");
//WAVEFORMS
// Need to put an option in here for pre/post 2024 data
// Define these at run time so we can use cable map
//addWaveBranches("VetoNu",2,4);
//addWaveBranches("VetoSt1",2,14);
//addWaveBranches("VetoSt2",2,6);
//addWaveBranches("Timing",4,8);
//addWaveBranches("Preshower",2,12);
//addWaveBranches("Calo",4,0);
//addWaveBranches("CaloHi", 4,16);
m_tree->Branch("ScintHit", &m_scintHit);
addCalibratedBranches("Calo",4,0);
@@ -471,7 +476,14 @@ StatusCode NtupleDumperAlg::initialize()
m_tree->Branch("truth_py", &m_truth_py);
m_tree->Branch("truth_pz", &m_truth_pz);
m_tree->Branch("truth_m", &m_truth_m);
m_tree->Branch("truth_theta", &m_truth_theta);
m_tree->Branch("truth_phi", &m_truth_phi);
m_tree->Branch("truth_pT", &m_truth_pT);
m_tree->Branch("truth_eta", &m_truth_eta);
m_tree->Branch("truth_barcode", &m_truth_barcode);
m_tree->Branch("truth_pdg", &m_truth_pdg);
m_tree->Branch("truth_isFiducial", &m_truth_isFiducial);
m_tree->Branch("truth_prod_x", &m_truth_prod_x);
m_tree->Branch("truth_prod_y", &m_truth_prod_y);
@@ -522,6 +534,51 @@ StatusCode NtupleDumperAlg::initialize()
ATH_CHECK(histSvc()->regTree("/HIST2/tree", m_tree));
// Register histograms
/*
m_HistRandomCharge[0] = new TH1F("hRandomCharge0", "Calo ch0 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[1] = new TH1F("hRandomCharge1", "Calo ch1 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[2] = new TH1F("hRandomCharge2", "Calo ch2 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[3] = new TH1F("hRandomCharge3", "Calo ch3 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[4] = new TH1F("hRandomCharge4", "VetoNu ch4 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[5] = new TH1F("hRandomCharge5", "VetoNu ch5 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[6] = new TH1F("hRandomCharge6", "Veto ch6 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[7] = new TH1F("hRandomCharge7", "Veto ch7 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[8] = new TH1F("hRandomCharge8", "Trig ch8 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[9] = new TH1F("hRandomCharge9", "Trig ch9 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[10] = new TH1F("hRandomCharge10", "Trig ch10 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[11] = new TH1F("hRandomCharge11", "Trig ch11 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[12] = new TH1F("hRandomCharge12", "Preshower ch12 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[13] = new TH1F("hRandomCharge13", "Preshower ch13 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[14] = new TH1F("hRandomCharge14", "Veto ch14 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[15] = new TH1F("hRandomCharge15", "Veto ch15 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[16] = new TH1F("hRandomCharge16", "CaloHi ch0 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[17] = new TH1F("hRandomCharge17", "CaloHi ch1 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[18] = new TH1F("hRandomCharge18", "CaloHi ch2 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
m_HistRandomCharge[19] = new TH1F("hRandomCharge19", "CaloHi ch3 Charge from Random Events;charge (pC);Events/bin", 100, -1.0, 1.0);
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge0", m_HistRandomCharge[0]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge1", m_HistRandomCharge[1]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge2", m_HistRandomCharge[2]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge3", m_HistRandomCharge[3]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge4", m_HistRandomCharge[4]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge5", m_HistRandomCharge[5]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge6", m_HistRandomCharge[6]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge7", m_HistRandomCharge[7]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge8", m_HistRandomCharge[8]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge9", m_HistRandomCharge[9]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge10", m_HistRandomCharge[10]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge11", m_HistRandomCharge[11]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge12", m_HistRandomCharge[12]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge13", m_HistRandomCharge[13]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge14", m_HistRandomCharge[14]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge15", m_HistRandomCharge[15]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge16", m_HistRandomCharge[16]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge17", m_HistRandomCharge[17]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge18", m_HistRandomCharge[18]));
ATH_CHECK(histSvc()->regHist("/HIST2/RandomCharge19", m_HistRandomCharge[19]));
*/
if (m_onlyBlinded){
ATH_MSG_INFO("Only events that would be blinded are saved in ntuple");
m_doBlinding = false;
@@ -531,10 +588,6 @@ StatusCode NtupleDumperAlg::initialize()
ATH_MSG_INFO("Blinding will NOT be enforced for real data.");
}
// Make sure charge histogram pointers are empty
for (unsigned int chan = 0; chan<max_chan; chan++)
m_HistRandomCharge[chan] = NULL;
return StatusCode::SUCCESS;
}
@@ -554,6 +607,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
if (truthEventContainer.isValid() && truthEventContainer->size() > 0)
{
isMC = true;
}
SG::ReadHandle<McEventCollection> mcEventContainer {m_mcEventContainer, ctx};
@@ -692,10 +746,8 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// Fill histograms
for (unsigned int chan = 0; chan<max_chan; chan++) {
// Only fill histograms that have been defined
if (m_HistRandomCharge[chan]) {
ATH_MSG_DEBUG("Filling random charge histogram chan= " << chan);
if (m_HistRandomCharge[chan])
m_HistRandomCharge[chan]->Fill(m_wave_raw_charge[chan]);
}
}
}
@@ -713,81 +765,18 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
}
}
if (m_doScintFilter) {
// Get channel mapping
auto wave_map = getChannelMap();
// filter events, but use waveform peak cuts instead of triggers, as triggers could miss signals slightly out of time
// digitizer channels described here: https://twiki.cern.ch/twiki/bin/viewauth/FASER/CaloScintillator
// Make calorimeter trigger, use CaloHi if it exists
bool calo_trig = false;
if (wave_map.contains(std::string("calo2"))) {
for (auto chan : wave_map[std::string("calo2")]) {
calo_trig |= m_wave_raw_peak[chan] > 3.0;
}
}
else if (wave_map.contains(std::string("calo"))) {
for (auto chan : wave_map[std::string("calo")]) {
calo_trig |= m_wave_raw_peak[chan] > 3.0;
}
}
bool calo_trig = (m_wave_raw_peak[0] > 3.0 || m_wave_raw_peak[1] > 3.0 || m_wave_raw_peak[2] > 3.0 || m_wave_raw_peak[3] > 3.0);
bool vetoNu_trig = (m_wave_raw_peak[4] > 25.0 && m_wave_raw_peak[5] > 25.0);
bool vetoSt2_trig = (m_wave_raw_peak[6] > 25.0 && m_wave_raw_peak[7] > 25.0);
bool timingScint_trig = ((m_wave_raw_peak[8] > 25.0 && m_wave_raw_peak[9] > 25.0) || (m_wave_raw_peak[10] > 25.0 && m_wave_raw_peak[11] > 25.0));
bool preshower_trig = (m_wave_raw_peak[12] > 3.0 && m_wave_raw_peak[13] > 3.0);
bool vetoSt1_trig = (m_wave_raw_peak[14] > 25.0);
// VetoNu trigger
bool vetoNu_trig = true;
if (wave_map.contains(std::string("vetonu"))) {
for (auto chan : wave_map[std::string("vetonu")]) {
vetoNu_trig &= m_wave_raw_peak[chan] > 25.0;
}
} else {
vetoNu_trig = false;
}
bool veto_OR_trig = (vetoNu_trig || vetoSt1_trig || vetoSt2_trig);
// Veto Station 2
bool vetoSt1_trig = true;
bool vetoSt2_trig = true;
if (wave_map.contains(std::string("veto"))) {
int ncount = 0;
for (auto chan : wave_map[std::string("veto")]) {
if (++ncount < 3) {
vetoSt2_trig &= m_wave_raw_peak[chan] > 25.0;
} else {
vetoSt1_trig &= m_wave_raw_peak[chan] > 25.0;
}
}
} else {
vetoSt1_trig = false;
vetoSt2_trig = false;
}
bool veto_OR_trig = vetoNu_trig || vetoSt1_trig || vetoSt2_trig;
// Timing layer (either top or bottom)
bool timingScint_trig = true;
bool timingScint_trig1 = true;
bool timingScint_trig2 = true;
if (wave_map.contains(std::string("trigger"))) {
int ncount = 0;
for (auto chan : wave_map[std::string("trigger")]) {
if (++ncount < 3) {
timingScint_trig1 &= m_wave_raw_peak[chan] > 25.0;
} else {
timingScint_trig2 &= m_wave_raw_peak[chan] > 25.0;
}
}
timingScint_trig = timingScint_trig1 || timingScint_trig2;
} else {
timingScint_trig = false;
}
bool preshower_trig = true;
if (wave_map.contains(std::string("preshower"))) {
for (auto chan : wave_map[std::string("preshower")]) {
preshower_trig &= m_wave_raw_peak[chan] > 3.0;
}
} else {
preshower_trig = false;
}
bool passes_ScintFilter = false;
if (calo_trig) {
passes_ScintFilter = true;
@@ -899,15 +888,21 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
// Find truth particle information
SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainer { m_truthParticleContainer, ctx };
if (truthParticleContainer.isValid() && truthParticleContainer->size() > 0) {
int ipart(0);
//int ipart(0);
for (auto particle : *truthParticleContainer) { // loop over first 10 truth particles (for non A' samples)
if (ipart++ > 9) break;
//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_pT.push_back(particle->p4().Pt());
m_truth_theta.push_back(particle->p4().Theta());
m_truth_eta.push_back(particle->p4().Eta());
m_truth_phi.push_back(particle->p4().Phi());
m_truth_barcode.push_back(particle->barcode());
m_truth_isFiducial.push_back(m_fiducialParticleTool->isFiducial(particle->barcode()));
m_truth_pdg.push_back(particle->pdgId());
if ( particle->hasProdVtx()) {
@@ -1041,6 +1036,7 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
m_calo_total_E_EM_fudged = m_calo_total_E_EM_fudged * m_caloMC_FudgeFactor;
}
// load in calibrated preshower container
SG::ReadHandle<xAOD::CalorimeterHitContainer> preshowerCalibratedContainer { m_preshowerCalibratedContainer, ctx };
ATH_CHECK(preshowerCalibratedContainer.isValid());
@@ -1054,7 +1050,8 @@ StatusCode NtupleDumperAlg::execute(const EventContext &ctx) const
ATH_MSG_DEBUG("Calibrated preshower: 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());
}
// enforce blinding such that events that miss the veto layers and have a large calo signal are skipped and not in the output root file
// Only blind colliding BCIDs (+/- 1)
bool blinded_event = false;
@@ -1590,8 +1587,6 @@ NtupleDumperAlg::clearTree() const
for(unsigned int ii=0;ii<max_chan;ii++) {
m_wave_localtime[ii]=0;
m_wave_triggertime[ii]=0;
m_wave_bcidtime[ii]=0;
m_wave_peak[ii]=0;
m_wave_width[ii]=0;
m_wave_charge[ii]=0;
@@ -1758,6 +1753,14 @@ NtupleDumperAlg::clearTree() const
m_truth_py.clear();
m_truth_pz.clear();
m_truth_m.clear();
m_truth_isFiducial.clear();
m_truth_pT.clear();
m_truth_theta.clear();
m_truth_phi.clear();
m_truth_eta.clear();
m_truth_barcode.clear();
m_truth_pdg.clear();
m_truth_prod_x.clear();
m_truth_prod_y.clear();