From 7ee981e05df254d403a1cf1f1a8cbcc1088029e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Sat, 4 Sep 2021 10:20:37 +0200 Subject: [PATCH] Add ECAL to AOD creation --- examples/aod/createO2tables.C | 18 ++++++++ examples/aod/createO2tables.h | 34 +++++++++++++-- examples/scripts/default_configfile.ini | 3 ++ rpythia8/rpythia8-box.cc | 4 +- src/CMakeLists.txt | 2 + src/ECALdetector.cc | 56 +++++++++++++++++++++++++ src/ECALdetector.hh | 38 +++++++++++++++++ 7 files changed, 151 insertions(+), 4 deletions(-) create mode 100644 src/ECALdetector.cc create mode 100644 src/ECALdetector.hh diff --git a/examples/aod/createO2tables.C b/examples/aod/createO2tables.C index 5a98014..cfc8b96 100644 --- a/examples/aod/createO2tables.C +++ b/examples/aod/createO2tables.C @@ -31,6 +31,7 @@ R__LOAD_LIBRARY(libDelphesO2) #include "TrackSmearer.hh" #include "TOFLayer.hh" #include "RICHdetector.hh" +#include "ECALdetector.hh" #include "MIDdetector.hh" #include "TrackUtils.hh" @@ -156,6 +157,9 @@ int createO2tables(const char* inputFile = "delphes.root", forward_rich_detector.setType(o2::delphes::RICHdetector::kForward); forward_rich_detector.setRadiusIn(forward_rich_radius_in); + // ECAL detector + o2::delphes::ECALdetector ecal_detector; + // MID detector o2::delphes::MIDdetector mid_detector; const bool isMID = mid_detector.setup(inputFileAccMuonPID); @@ -173,6 +177,7 @@ int createO2tables(const char* inputFile = "delphes.root", MakeTreeO2ftof(); MakeTreeO2rich(); MakeTreeO2frich(); + MakeTreeO2ecal(); MakeTreeO2mid(); MakeTreeO2collision(); MakeTreeO2collisionExtra(); @@ -259,6 +264,19 @@ int createO2tables(const char* inputFile = "delphes.root", dNdEta += 1.f; } FillTree(kMcParticle); + + // info for the ECAL + std::array pos; + float energy = 0; + if (ecal_detector.makeSignal(*particle, pos, energy)) { + ecal.fIndexCollisions = ientry + eventOffset; + ecal.fIndexMcParticles = TMath::Abs(iparticle + fOffsetLabel); + ecal.fEnergy = energy; + ecal.fPosX = pos[0]; + ecal.fPosY = pos[1]; + ecal.fPosZ = pos[2]; + FillTree(kA3ECAL); + } } dNdEta = 0.5f * dNdEta / multEtaRange; fOffsetLabel += particles->GetEntries(); diff --git a/examples/aod/createO2tables.h b/examples/aod/createO2tables.h index c0a0a7f..d14a995 100644 --- a/examples/aod/createO2tables.h +++ b/examples/aod/createO2tables.h @@ -32,6 +32,7 @@ enum TreeIndex { // Index of the output trees kFRICH, kMID, kFTOF, + kA3ECAL, kTrees }; @@ -68,7 +69,8 @@ const TString TreeName[kTrees] = {"O2collision", "O2rich", "O2frich", "O2mid", - "O2ftof"}; + "O2ftof", + "O2a3ecal"}; const TString TreeTitle[kTrees] = {"Collision tree", "Collision extra", @@ -100,7 +102,8 @@ const TString TreeTitle[kTrees] = {"Collision tree", "RICH info", "Forward RICH info", "MID info", - "Forward TOF info"}; + "Forward TOF info", + "ALICE3 ECAL"}; TTree* Trees[kTrees] = {nullptr}; // Array of created TTrees TTree* CreateTree(TreeIndex t) @@ -373,7 +376,7 @@ struct { Float_t fRICHNsigmaPi = -999.f; /// Nsigma for Pi Float_t fRICHNsigmaKa = -999.f; /// Nsigma for Ka Float_t fRICHNsigmaPr = -999.f; /// Nsigma for Pr -} rich, frich; //! structure to keep RICH info +} rich, frich; //! structure to keep RICH info void MakeTreeO2rich() { @@ -470,6 +473,31 @@ void MakeTreeO2ftof() tFTOF->SetBasketSize("*", fBasketSizeTracks); } +struct { + // ALICE3 ECAL data + Int_t fIndexCollisions = -1; /// Collision ID + Int_t fIndexMcParticles = -1; /// Particle ID + Int_t fIndexTracks = -1; /// Track ID + + Float_t fEnergy = -999.f; /// Energy + Float_t fPosX = -999.f; /// Position in X + Float_t fPosY = -999.f; /// Position in Y + Float_t fPosZ = -999.f; /// Position in Z +} ecal; //! structure to keep ECAL info + +void MakeTreeO2ecal() +{ + TTree* tECAL = CreateTree(kA3ECAL); + tECAL->Branch("fIndexCollisions", &ecal.fIndexCollisions, "fIndexCollisions/I"); + tECAL->Branch("fIndexMcParticles", &ecal.fIndexMcParticles, "fIndexMcParticles/I"); + tECAL->Branch("fIndexTracks", &ecal.fIndexTracks, "fIndexTracks/I"); + tECAL->Branch("fEnergy", &ecal.fEnergy, "fEnergy/F"); + tECAL->Branch("fPosX", &ecal.fPosX, "fPosX/F"); + tECAL->Branch("fPosY", &ecal.fPosY, "fPosY/F"); + tECAL->Branch("fPosZ", &ecal.fPosZ, "fPosZ/F"); + tECAL->SetBasketSize("*", fBasketSizeTracks); +} + struct { // MC particle diff --git a/examples/scripts/default_configfile.ini b/examples/scripts/default_configfile.ini index de78eb8..528e71c 100644 --- a/examples/scripts/default_configfile.ini +++ b/examples/scripts/default_configfile.ini @@ -79,6 +79,9 @@ custom_gen = rpythia8-gun --pdg 421 --px 1. --py 0. --pz 0. --xProd 0. --yProd 0 [BOX_pion] custom_gen = rpythia8-box --pdg 211 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 +[BOX_photons] +custom_gen = rpythia8-box --pdg 22 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 + [BOX_proton] custom_gen = rpythia8-box --pdg 2212 --etamin -2. --etamax 2. --phimin 0. --phimax 6.28 --pmin 0 --pmax 10 --xProd 0. --yProd 0. --zProd 0. --decay --npart 100 --config $DELPHESO2_ROOT/examples/pythia8/pythia_nuclei.cfg --background-config $O2_ROOT/share/Generators/egconfig/pythia8_hi.cfg diff --git a/rpythia8/rpythia8-box.cc b/rpythia8/rpythia8-box.cc index 3c1d155..859c7e1 100644 --- a/rpythia8/rpythia8-box.cc +++ b/rpythia8/rpythia8-box.cc @@ -81,9 +81,11 @@ int main(int argc, char** argv) if (!pythia.particleData.isLepton(pdg) && !pythia.particleData.isHadron(pdg) && !pythia.particleData.isResonance(pdg)) { - if (abs(pdg) < 1000000000) { + if (abs(pdg) < 1000000000 && pdg != 22) { std::cout << "Error: invalid PDG code \"" << pdg << "\"" << std::endl; return 1; + } else if (pdg == 22) { + std::cout << "PDG code \"" << pdg << "\" stands for a gamma" << std::endl; } else { std::cout << "PDG code \"" << pdg << "\" stands for a nucleous" << std::endl; } diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index eef564c..6e40024 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,6 +6,7 @@ set(SOURCES TrackUtils.cc TOFLayer.cc RICHdetector.cc + ECALdetector.cc MIDdetector.cc PreShower.cc ) @@ -16,6 +17,7 @@ set(HEADERS TrackUtils.hh TOFLayer.hh RICHdetector.hh + ECALdetector.hh MIDdetector.hh PreShower.hh ) diff --git a/src/ECALdetector.cc b/src/ECALdetector.cc new file mode 100644 index 0000000..43c2784 --- /dev/null +++ b/src/ECALdetector.cc @@ -0,0 +1,56 @@ +/// @author: Yuri Kharlov +/// @author: Nicolo' Jacazio +/// @since 04/09/2021 + +#include "ECALdetector.hh" +#include "TDatabasePDG.h" +#include "TRandom.h" +#include "TLorentzVector.h" + +namespace o2 +{ +namespace delphes +{ + +/*****************************************************************/ + +void ECALdetector::setup(float resoEA, float resoEB, float resoEC, float resoPos) +{ +} + +/*****************************************************************/ + +bool ECALdetector::hasECAL(const Track& track) const +{ + auto x = track.XOuter * 0.1; // [cm] + auto y = track.YOuter * 0.1; // [cm] + auto z = track.ZOuter * 0.1; // [cm] + /** check if hit **/ + bool ishit = false; + auto r = hypot(x, y); + ishit = (fabs(r - mRadius) < 0.001 && fabs(z) < mLength); + if (!ishit) + return false; + auto particle = (GenParticle*)track.Particle.GetObject(); + return true; +} + +/*****************************************************************/ +bool ECALdetector::makeSignal(const GenParticle& particle, std::array& pos, float& energy) const +{ + TLorentzVector vec; + const int pid = particle.PID; + if (pid != 22) { + return false; + } + pos[0] = 0.f; + pos[1] = 0.f; + pos[2] = 0.f; + energy = 0.f; + return true; +} + +/*****************************************************************/ + +} // namespace delphes +} // namespace o2 diff --git a/src/ECALdetector.hh b/src/ECALdetector.hh new file mode 100644 index 0000000..b9572aa --- /dev/null +++ b/src/ECALdetector.hh @@ -0,0 +1,38 @@ +/// @author: Roberto Preghenella +/// @email: preghenella@bo.infn.it + +#ifndef _DelphesO2_ECALdetector_h_ +#define _DelphesO2_ECALdetector_h_ + +#include "classes/DelphesClasses.h" + +namespace o2 +{ +namespace delphes +{ + +class ECALdetector +{ + + public: + ECALdetector() = default; + ~ECALdetector() = default; + + void setup(float resoEA, float resoEB, float resoEC, float resoPos); + bool hasECAL(const Track& track) const; + bool makeSignal(const GenParticle& particle, std::array& pos, float& energy) const; + + protected: + float mRadius = 120.; // [cm] + float mLength = 200.; // [cm] + + float mEnergyResolutionA = 0; // [] + float mEnergyResolutionB = 0; // [] + float mEnergyResolutionC = 0; // [] + float mPositionResolution = 0; // [] +}; + +} // namespace delphes +} // namespace o2 + +#endif /** _DelphesO2_ECALdetector_h_ **/