From 54dfe81ee22e63875d9cf65124d8779837b7d5d7 Mon Sep 17 00:00:00 2001 From: Yuri Kharlov Date: Fri, 10 Sep 2021 18:50:46 +0200 Subject: [PATCH 1/5] First ECAL implementation from Nicolo Jacazio --- examples/aod/createO2tables.C | 18 ++++++++ examples/aod/createO2tables.h | 32 +++++++++++++- 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, 150 insertions(+), 3 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 3546be5..2459033 100644 --- a/examples/aod/createO2tables.C +++ b/examples/aod/createO2tables.C @@ -32,6 +32,7 @@ R__LOAD_LIBRARY(libDelphesO2) #include "TrackSmearer.hh" #include "TOFLayer.hh" #include "RICHdetector.hh" +#include "ECALdetector.hh" #include "MIDdetector.hh" #include "TrackUtils.hh" @@ -167,6 +168,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); @@ -183,6 +187,7 @@ int createO2tables(const char* inputFile = "delphes.root", MakeTreeO2trackExtra(); MakeTreeO2ftof(); MakeTreeO2rich(); + MakeTreeO2ecal(); MakeTreeO2frich(); MakeTreeO2mid(); MakeTreeO2collision(); @@ -270,6 +275,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); + } if constexpr (debug_qa) { if (!debugEffDenPart[particle->PID]) { debugEffDenPart[particle->PID] = new TH1F(Form("denPart%i", particle->PID), Form("denPart%i;#it{p}_{T} (GeV/#it{c})", particle->PID), 1000, 0, 10); diff --git a/examples/aod/createO2tables.h b/examples/aod/createO2tables.h index c0a0a7f..783d1e4 100644 --- a/examples/aod/createO2tables.h +++ b/examples/aod/createO2tables.h @@ -19,6 +19,7 @@ enum TreeIndex { // Index of the output trees kV0s, kCascades, kTOF, + kA3ECAL, kMcParticle, kMcCollision, kMcTrackLabel, @@ -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) @@ -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 256c669..bbf8ab6 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_ **/ From 87730861976a9c4dbffd8b1b1203eb92de0028ea Mon Sep 17 00:00:00 2001 From: Yuri Kharlov Date: Sat, 11 Sep 2021 23:29:10 +0200 Subject: [PATCH 2/5] Energy and position smearing is added --- src/ECALdetector.cc | 134 +++++++++++++++++++++++++++++++------------- src/ECALdetector.hh | 43 +++++++------- 2 files changed, 118 insertions(+), 59 deletions(-) diff --git a/src/ECALdetector.cc b/src/ECALdetector.cc index 43c2784..ae363a7 100644 --- a/src/ECALdetector.cc +++ b/src/ECALdetector.cc @@ -9,48 +9,102 @@ namespace o2 { -namespace delphes -{ + namespace delphes + { -/*****************************************************************/ + /*****************************************************************/ -void ECALdetector::setup(float resoEA, float resoEB, float resoEC, float resoPos) -{ -} + void ECALdetector::setup(float resoEA, float resoEB, float resoEC, float resoXA, float resoXB) + { + mEnergyResolutionA = resoEA; + mEnergyResolutionB = resoEB; + mEnergyResolutionC = resoEC; + mPositionResolutionA = resoXA; + mPositionResolutionB = resoXB; + } -/*****************************************************************/ + /*****************************************************************/ -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 + 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 + bool ECALdetector::makeSignal(const GenParticle& particle, TLorentzVector &pECAL, float & posZ, float & posPhi) + { + const int pid = particle.PID; + if (pid != 22) { + return false; + } + + TLorentzVector p4True = particle.P4(); // true 4-momentum + Float_t vtX = particle.X; // particle vertex X + Float_t vtY = particle.Y; // particle vertex Y + Float_t vtZ = particle.Z; // particle vertex Z + + posPhi = particle.Phi; // azimuth angle of a photon hit + posZ = mRadius * particle.CtgTheta; // z-coodrinate of a photon hit + + Double_t eSmeared = smearPhotonE(particle.E); // smeared photon energy + TLorentzVector p4Smeared = smearPhotonP4(p4True); + + return true; + } + + /*****************************************************************/ + TLorentzVector ECALdetector::smearPhotonP4(const TLorentzVector pTrue) + { + // This function smears the photon 4-momentum from the true one via applying + // parametrized energy and coordinate resolution + + // Get true energy from true 4-momentum and smear this energy + Double_t eTrue = pTrue.E(); + Double_t eSmeared = smearPhotonE(eTrue); + // Smear direction of 3-vector + Double_t phi = pTrue.Phi() + gRandom->Gaus(0.,sigmaX(eTrue)/mRadius); + Double_t theta = pTrue.Theta() + gRandom->Gaus(0.,sigmaX(eTrue)/mRadius); + // Calculate smeared components of 3-vector + Double_t pxSmeared = eSmeared*TMath::Cos(phi)*TMath::Sin(theta); + Double_t pySmeared = eSmeared*TMath::Sin(phi)*TMath::Sin(theta); + Double_t pzSmeared = eSmeared*TMath::Cos(theta); + // Construct new 4-momentum from smeared energy and 3-momentum + TLorentzVector pSmeared(pxSmeared,pySmeared,pzSmeared,eSmeared); + return pSmeared; + } + /*****************************************************************/ + Double_t ECALdetector::sigmaX(const Double_t eTrue) + { + // Calculate sigma of photon coordinate smearing [cm] + // E is the photon energy + Double_t dX = sqrt(mPositionResolutionA*mPositionResolutionA + mPositionResolutionB*mPositionResolutionB/eTrue); + return dX; + } + /*****************************************************************/ + Double_t ECALdetector::smearPhotonE(const Double_t eTrue) + { + // Smear a photon energy eTrue according to a Gaussian distribution with energy resolution parameters + // sigma of Gaussian smearing is calculated from parameters A,B,C and true energy + + Double_t sigmaE = eTrue * sqrt(mEnergyResolutionA*mEnergyResolutionA/eTrue/eTrue + + mEnergyResolutionB*mEnergyResolutionB/eTrue + + mEnergyResolutionC*mEnergyResolutionC); + Double_t eSmeared = gRandom->Gaus(eTrue,sigmaE); + if (eSmeared < 0) eSmeared = 0; + return eSmeared; + } + /*****************************************************************/ + + } // namespace delphes } // namespace o2 diff --git a/src/ECALdetector.hh b/src/ECALdetector.hh index b9572aa..5354981 100644 --- a/src/ECALdetector.hh +++ b/src/ECALdetector.hh @@ -8,31 +8,36 @@ namespace o2 { -namespace delphes -{ + namespace delphes + { -class ECALdetector -{ + class ECALdetector + { + + public: + ECALdetector() = default; + ~ECALdetector() = default; - public: - ECALdetector() = default; - ~ECALdetector() = default; + void setup(float resoEA, float resoEB, float resoEC, float resoPosA, float resoPosB); + bool hasECAL(const Track& track) const; + bool makeSignal(const GenParticle& particle, TLorentzVector &pECAL, float & posZ, float & posPhi); - 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: + Double_t smearPhotonE (const Double_t eTrue); + Double_t sigmaX (const Double_t eTrue); + TLorentzVector smearPhotonP4(const TLorentzVector pTrue); - protected: - float mRadius = 120.; // [cm] - float mLength = 200.; // [cm] + float mRadius = 120.; // ECAL barrel inner radius [cm] + float mLength = 200.; // ECAL half-length along beam axis [cm] - float mEnergyResolutionA = 0; // [] - float mEnergyResolutionB = 0; // [] - float mEnergyResolutionC = 0; // [] - float mPositionResolution = 0; // [] -}; + float mEnergyResolutionA = 0.002; // parameter A of energy resolution in GeV + float mEnergyResolutionB = 0.02; // parameter B of energy resolution in GeV^{1/2} + float mEnergyResolutionC = 0.10; // parameter C of energy resolution + float mPositionResolutionA = 0.15; // parameter A of coordinate resolution in cm + float mPositionResolutionB = 0.30; // parameter B of coordinate resolution in cm*GeV^{1/2} + }; -} // namespace delphes + } // namespace delphes } // namespace o2 #endif /** _DelphesO2_ECALdetector_h_ **/ From 7314886414b99bbf20b2c5f3b532842b52fe3fe2 Mon Sep 17 00:00:00 2001 From: Yuri Kharlov Date: Thu, 23 Sep 2021 13:18:20 +0200 Subject: [PATCH 3/5] ECAL response updated, bug fixed --- examples/aod/createO2tables.C | 34 ++++++++++++++++++++-------------- examples/aod/createO2tables.h | 22 +++++++++++++--------- src/ECALdetector.cc | 20 +++++++++++++++----- 3 files changed, 48 insertions(+), 28 deletions(-) diff --git a/examples/aod/createO2tables.C b/examples/aod/createO2tables.C index 2459033..96d6e57 100644 --- a/examples/aod/createO2tables.C +++ b/examples/aod/createO2tables.C @@ -39,12 +39,12 @@ R__LOAD_LIBRARY(libDelphesO2) #include "createO2tables.h" // Detector parameters -const double Bz = 0.2; // [T] +const double Bz = 5.e-1; // TOF -const double tof_radius = 100.; // [cm] Radius of the TOF detector (used to compute acceptance) -const double tof_length = 200.; // [cm] Length of the TOF detector (used to compute acceptance) -const double tof_sigmat = 0.02; // [ns] Resolution of the TOF detector -const double tof_sigmat0 = 0.2; // [ns] Time spread of the vertex +const double tof_radius = 100.; +const double tof_length = 200.; +const double tof_sigmat = 0.020; +const double tof_sigmat0 = 0.20; // Forward TOF const double forward_tof_radius = 100.; // [cm] Radius of the Forward TOF detector (used to compute acceptance) const double forward_tof_radius_in = 10.; // [cm] Inner radius of the Forward TOF detector (used to compute acceptance) @@ -71,7 +71,7 @@ const char* inputFileAccMuonPID = "muonAccEffPID.root"; // Simulation parameters constexpr bool do_vertexing = true; // Vertexing with the O2 -constexpr bool enable_nuclei = true; // Nuclei LUTs +constexpr bool enable_nuclei = true; constexpr bool debug_qa = false; // Debug QA histograms int createO2tables(const char* inputFile = "delphes.root", @@ -277,15 +277,21 @@ int createO2tables(const char* inputFile = "delphes.root", FillTree(kMcParticle); // info for the ECAL - std::array pos; - float energy = 0; - if (ecal_detector.makeSignal(*particle, pos, energy)) { - ecal.fIndexCollisions = ientry + eventOffset; + float posZ, posPhi; + TLorentzVector pECAL; + if (ecal_detector.makeSignal(*particle, pECAL, posZ, posPhi)) { // to be updated 13.09.2021 + printf("ECAL particle: pid=%d, p=(%g,%g,%g,%g)\n", + particle->PID,particle->Px,particle->Py,particle->Pz,particle->E); + printf("ECAL p=(%g,%g,%g,%g)\n", + pECAL.Px(),pECAL.Py(),pECAL.Pz(),pECAL.E()); + 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]; + ecal.fPx = pECAL.Px(); + ecal.fPy = pECAL.Py(); + ecal.fPz = pECAL.Pz(); + ecal.fE = pECAL.E(); + ecal.fPosZ = 0.f; + ecal.fPosPhi = 0.f; FillTree(kA3ECAL); } if constexpr (debug_qa) { diff --git a/examples/aod/createO2tables.h b/examples/aod/createO2tables.h index 783d1e4..ebec2c2 100644 --- a/examples/aod/createO2tables.h +++ b/examples/aod/createO2tables.h @@ -19,7 +19,6 @@ enum TreeIndex { // Index of the output trees kV0s, kCascades, kTOF, - kA3ECAL, kMcParticle, kMcCollision, kMcTrackLabel, @@ -33,6 +32,7 @@ enum TreeIndex { // Index of the output trees kFRICH, kMID, kFTOF, + kA3ECAL, kTrees }; @@ -479,10 +479,12 @@ struct { 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 + Double_t fPx = -1.e10; /// px + Double_t fPy = -1.e10; /// py + Double_t fPz = -1.e10; /// pz + Double_t fE = -1.e10; /// E + Float_t fPosZ = -999.f; /// Position in Z + Float_t fPosPhi = -999.f; /// Position in phi } ecal; //! structure to keep ECAL info void MakeTreeO2ecal() @@ -491,10 +493,12 @@ void MakeTreeO2ecal() 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->Branch("fPx" , &ecal.fPx , "fPx/D"); + tECAL->Branch("fPy" , &ecal.fPy , "fPy/D"); + tECAL->Branch("fPz" , &ecal.fPz , "fPz/D"); + tECAL->Branch("fE" , &ecal.fE , "fE/D" ); + tECAL->Branch("fPosZ" , &ecal.fPosZ , "fPosZ/F"); + tECAL->Branch("fPosPhi", &ecal.fPosPhi, "fPosPhi/F"); tECAL->SetBasketSize("*", fBasketSizeTracks); } diff --git a/src/ECALdetector.cc b/src/ECALdetector.cc index ae363a7..a3aec31 100644 --- a/src/ECALdetector.cc +++ b/src/ECALdetector.cc @@ -41,15 +41,25 @@ namespace o2 } /*****************************************************************/ - // bool ECALdetector::makeSignal(const GenParticle& particle, std::array& pos, float& energy) const - bool ECALdetector::makeSignal(const GenParticle& particle, TLorentzVector &pECAL, float & posZ, float & posPhi) + bool ECALdetector::makeSignal(const GenParticle& particle, + TLorentzVector &p4ECAL, + float & posZ, + float & posPhi) { + // Simulate fast response of ECAL to photons: + // take generated particle as input and calculate its smeared 4-momentum p4ECAL + // and hit coordinates posZ, posPhi + const int pid = particle.PID; - if (pid != 22) { + if (pid != 22) { // only photons are treated so far. e+- and MIPs will be added later. return false; } TLorentzVector p4True = particle.P4(); // true 4-momentum + if (TMath::Abs(p4True.Eta())>4) { // ECAL acceptance is rougly limited by |eta|<4 + return false; + } + Float_t vtX = particle.X; // particle vertex X Float_t vtY = particle.Y; // particle vertex Y Float_t vtZ = particle.Z; // particle vertex Z @@ -58,8 +68,8 @@ namespace o2 posZ = mRadius * particle.CtgTheta; // z-coodrinate of a photon hit Double_t eSmeared = smearPhotonE(particle.E); // smeared photon energy - TLorentzVector p4Smeared = smearPhotonP4(p4True); - + p4ECAL = smearPhotonP4(p4True); + return true; } From 6a384e2ae0569e750c8ae5362315024e5dc6e5b4 Mon Sep 17 00:00:00 2001 From: Yuri Kharlov Date: Thu, 23 Sep 2021 18:56:16 +0200 Subject: [PATCH 4/5] Merging with master --- examples/aod/createO2tables.C | 28 ++--- examples/aod/createO2tables.h | 22 ++-- src/ECALdetector.cc | 217 +++++++++++++++++----------------- src/ECALdetector.hh | 46 +++---- 4 files changed, 157 insertions(+), 156 deletions(-) diff --git a/examples/aod/createO2tables.C b/examples/aod/createO2tables.C index 8e65306..21b912a 100644 --- a/examples/aod/createO2tables.C +++ b/examples/aod/createO2tables.C @@ -71,10 +71,10 @@ const double forward_rich_sigma = 7.e-3; // Resolution of the Forward RIC const char* inputFileAccMuonPID = "muonAccEffPID.root"; // Simulation parameters -constexpr bool do_vertexing = true; // Vertexing with the O2 +constexpr bool do_vertexing = true; // Vertexing with the O2 constexpr bool enable_nuclei = true; -constexpr bool debug_qa = false; // Debug QA histograms -constexpr int tof_mismatch = 0; // Flag to configure the TOF mismatch running mode: 0 off, 1 create, 2 use +constexpr bool debug_qa = false; // Debug QA histograms +constexpr int tof_mismatch = 0; // Flag to configure the TOF mismatch running mode: 0 off, 1 create, 2 use int createO2tables(const char* inputFile = "delphes.root", const char* outputFile = "AODRun5.root", @@ -295,18 +295,18 @@ int createO2tables(const char* inputFile = "delphes.root", float posZ, posPhi; TLorentzVector pECAL; if (ecal_detector.makeSignal(*particle, pECAL, posZ, posPhi)) { // to be updated 13.09.2021 - printf("ECAL particle: pid=%d, p=(%g,%g,%g,%g)\n", - particle->PID,particle->Px,particle->Py,particle->Pz,particle->E); - printf("ECAL p=(%g,%g,%g,%g)\n", - pECAL.Px(),pECAL.Py(),pECAL.Pz(),pECAL.E()); - ecal.fIndexCollisions = ientry + eventOffset; + printf("ECAL particle: pid=%d, p=(%g,%g,%g,%g)\n", + particle->PID, particle->Px, particle->Py, particle->Pz, particle->E); + printf("ECAL p=(%g,%g,%g,%g)\n", + pECAL.Px(), pECAL.Py(), pECAL.Pz(), pECAL.E()); + ecal.fIndexCollisions = ientry + eventOffset; ecal.fIndexMcParticles = TMath::Abs(iparticle + fOffsetLabel); - ecal.fPx = pECAL.Px(); - ecal.fPy = pECAL.Py(); - ecal.fPz = pECAL.Pz(); - ecal.fE = pECAL.E(); - ecal.fPosZ = 0.f; - ecal.fPosPhi = 0.f; + ecal.fPx = pECAL.Px(); + ecal.fPy = pECAL.Py(); + ecal.fPz = pECAL.Pz(); + ecal.fE = pECAL.E(); + ecal.fPosZ = 0.f; + ecal.fPosPhi = 0.f; FillTree(kA3ECAL); } if constexpr (debug_qa) { diff --git a/examples/aod/createO2tables.h b/examples/aod/createO2tables.h index 422e5fc..0b8d88b 100644 --- a/examples/aod/createO2tables.h +++ b/examples/aod/createO2tables.h @@ -483,12 +483,12 @@ struct { Int_t fIndexMcParticles = -1; /// Particle ID Int_t fIndexTracks = -1; /// Track ID - Double_t fPx = -1.e10; /// px - Double_t fPy = -1.e10; /// py - Double_t fPz = -1.e10; /// pz - Double_t fE = -1.e10; /// E - Float_t fPosZ = -999.f; /// Position in Z - Float_t fPosPhi = -999.f; /// Position in phi + Double_t fPx = -1.e10; /// px + Double_t fPy = -1.e10; /// py + Double_t fPz = -1.e10; /// pz + Double_t fE = -1.e10; /// E + Float_t fPosZ = -999.f; /// Position in Z + Float_t fPosPhi = -999.f; /// Position in phi } ecal; //! structure to keep ECAL info void MakeTreeO2ecal() @@ -497,11 +497,11 @@ void MakeTreeO2ecal() tECAL->Branch("fIndexCollisions", &ecal.fIndexCollisions, "fIndexCollisions/I"); tECAL->Branch("fIndexMcParticles", &ecal.fIndexMcParticles, "fIndexMcParticles/I"); tECAL->Branch("fIndexTracks", &ecal.fIndexTracks, "fIndexTracks/I"); - tECAL->Branch("fPx" , &ecal.fPx , "fPx/D"); - tECAL->Branch("fPy" , &ecal.fPy , "fPy/D"); - tECAL->Branch("fPz" , &ecal.fPz , "fPz/D"); - tECAL->Branch("fE" , &ecal.fE , "fE/D" ); - tECAL->Branch("fPosZ" , &ecal.fPosZ , "fPosZ/F"); + tECAL->Branch("fPx", &ecal.fPx, "fPx/D"); + tECAL->Branch("fPy", &ecal.fPy, "fPy/D"); + tECAL->Branch("fPz", &ecal.fPz, "fPz/D"); + tECAL->Branch("fE", &ecal.fE, "fE/D"); + tECAL->Branch("fPosZ", &ecal.fPosZ, "fPosZ/F"); tECAL->Branch("fPosPhi", &ecal.fPosPhi, "fPosPhi/F"); tECAL->SetBasketSize("*", fBasketSizeTracks); } diff --git a/src/ECALdetector.cc b/src/ECALdetector.cc index a3aec31..ffe098d 100644 --- a/src/ECALdetector.cc +++ b/src/ECALdetector.cc @@ -9,112 +9,113 @@ namespace o2 { - namespace delphes - { - - /*****************************************************************/ - - void ECALdetector::setup(float resoEA, float resoEB, float resoEC, float resoXA, float resoXB) - { - mEnergyResolutionA = resoEA; - mEnergyResolutionB = resoEB; - mEnergyResolutionC = resoEC; - mPositionResolutionA = resoXA; - mPositionResolutionB = resoXB; - } - - /*****************************************************************/ - - 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, - TLorentzVector &p4ECAL, - float & posZ, - float & posPhi) - { - // Simulate fast response of ECAL to photons: - // take generated particle as input and calculate its smeared 4-momentum p4ECAL - // and hit coordinates posZ, posPhi - - const int pid = particle.PID; - if (pid != 22) { // only photons are treated so far. e+- and MIPs will be added later. - return false; - } - - TLorentzVector p4True = particle.P4(); // true 4-momentum - if (TMath::Abs(p4True.Eta())>4) { // ECAL acceptance is rougly limited by |eta|<4 - return false; - } - - Float_t vtX = particle.X; // particle vertex X - Float_t vtY = particle.Y; // particle vertex Y - Float_t vtZ = particle.Z; // particle vertex Z - - posPhi = particle.Phi; // azimuth angle of a photon hit - posZ = mRadius * particle.CtgTheta; // z-coodrinate of a photon hit - - Double_t eSmeared = smearPhotonE(particle.E); // smeared photon energy - p4ECAL = smearPhotonP4(p4True); - - return true; - } - - /*****************************************************************/ - TLorentzVector ECALdetector::smearPhotonP4(const TLorentzVector pTrue) - { - // This function smears the photon 4-momentum from the true one via applying - // parametrized energy and coordinate resolution - - // Get true energy from true 4-momentum and smear this energy - Double_t eTrue = pTrue.E(); - Double_t eSmeared = smearPhotonE(eTrue); - // Smear direction of 3-vector - Double_t phi = pTrue.Phi() + gRandom->Gaus(0.,sigmaX(eTrue)/mRadius); - Double_t theta = pTrue.Theta() + gRandom->Gaus(0.,sigmaX(eTrue)/mRadius); - // Calculate smeared components of 3-vector - Double_t pxSmeared = eSmeared*TMath::Cos(phi)*TMath::Sin(theta); - Double_t pySmeared = eSmeared*TMath::Sin(phi)*TMath::Sin(theta); - Double_t pzSmeared = eSmeared*TMath::Cos(theta); - // Construct new 4-momentum from smeared energy and 3-momentum - TLorentzVector pSmeared(pxSmeared,pySmeared,pzSmeared,eSmeared); - return pSmeared; - } - /*****************************************************************/ - Double_t ECALdetector::sigmaX(const Double_t eTrue) - { - // Calculate sigma of photon coordinate smearing [cm] - // E is the photon energy - Double_t dX = sqrt(mPositionResolutionA*mPositionResolutionA + mPositionResolutionB*mPositionResolutionB/eTrue); - return dX; - } - /*****************************************************************/ - Double_t ECALdetector::smearPhotonE(const Double_t eTrue) - { - // Smear a photon energy eTrue according to a Gaussian distribution with energy resolution parameters - // sigma of Gaussian smearing is calculated from parameters A,B,C and true energy - - Double_t sigmaE = eTrue * sqrt(mEnergyResolutionA*mEnergyResolutionA/eTrue/eTrue + - mEnergyResolutionB*mEnergyResolutionB/eTrue + - mEnergyResolutionC*mEnergyResolutionC); - Double_t eSmeared = gRandom->Gaus(eTrue,sigmaE); - if (eSmeared < 0) eSmeared = 0; - return eSmeared; - } - /*****************************************************************/ - - } // namespace delphes +namespace delphes +{ + +/*****************************************************************/ + +void ECALdetector::setup(float resoEA, float resoEB, float resoEC, float resoXA, float resoXB) +{ + mEnergyResolutionA = resoEA; + mEnergyResolutionB = resoEB; + mEnergyResolutionC = resoEC; + mPositionResolutionA = resoXA; + mPositionResolutionB = resoXB; +} + +/*****************************************************************/ + +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, + TLorentzVector& p4ECAL, + float& posZ, + float& posPhi) +{ + // Simulate fast response of ECAL to photons: + // take generated particle as input and calculate its smeared 4-momentum p4ECAL + // and hit coordinates posZ, posPhi + + const int pid = particle.PID; + if (pid != 22) { // only photons are treated so far. e+- and MIPs will be added later. + return false; + } + + TLorentzVector p4True = particle.P4(); // true 4-momentum + if (TMath::Abs(p4True.Eta()) > 4) { // ECAL acceptance is rougly limited by |eta|<4 + return false; + } + + Float_t vtX = particle.X; // particle vertex X + Float_t vtY = particle.Y; // particle vertex Y + Float_t vtZ = particle.Z; // particle vertex Z + + posPhi = particle.Phi; // azimuth angle of a photon hit + posZ = mRadius * particle.CtgTheta; // z-coodrinate of a photon hit + + Double_t eSmeared = smearPhotonE(particle.E); // smeared photon energy + p4ECAL = smearPhotonP4(p4True); + + return true; +} + +/*****************************************************************/ +TLorentzVector ECALdetector::smearPhotonP4(const TLorentzVector pTrue) +{ + // This function smears the photon 4-momentum from the true one via applying + // parametrized energy and coordinate resolution + + // Get true energy from true 4-momentum and smear this energy + Double_t eTrue = pTrue.E(); + Double_t eSmeared = smearPhotonE(eTrue); + // Smear direction of 3-vector + Double_t phi = pTrue.Phi() + gRandom->Gaus(0., sigmaX(eTrue) / mRadius); + Double_t theta = pTrue.Theta() + gRandom->Gaus(0., sigmaX(eTrue) / mRadius); + // Calculate smeared components of 3-vector + Double_t pxSmeared = eSmeared * TMath::Cos(phi) * TMath::Sin(theta); + Double_t pySmeared = eSmeared * TMath::Sin(phi) * TMath::Sin(theta); + Double_t pzSmeared = eSmeared * TMath::Cos(theta); + // Construct new 4-momentum from smeared energy and 3-momentum + TLorentzVector pSmeared(pxSmeared, pySmeared, pzSmeared, eSmeared); + return pSmeared; +} +/*****************************************************************/ +Double_t ECALdetector::sigmaX(const Double_t eTrue) +{ + // Calculate sigma of photon coordinate smearing [cm] + // E is the photon energy + Double_t dX = sqrt(mPositionResolutionA * mPositionResolutionA + mPositionResolutionB * mPositionResolutionB / eTrue); + return dX; +} +/*****************************************************************/ +Double_t ECALdetector::smearPhotonE(const Double_t eTrue) +{ + // Smear a photon energy eTrue according to a Gaussian distribution with energy resolution parameters + // sigma of Gaussian smearing is calculated from parameters A,B,C and true energy + + Double_t sigmaE = eTrue * sqrt(mEnergyResolutionA * mEnergyResolutionA / eTrue / eTrue + + mEnergyResolutionB * mEnergyResolutionB / eTrue + + mEnergyResolutionC * mEnergyResolutionC); + Double_t eSmeared = gRandom->Gaus(eTrue, sigmaE); + if (eSmeared < 0) + eSmeared = 0; + return eSmeared; +} +/*****************************************************************/ + +} // namespace delphes } // namespace o2 diff --git a/src/ECALdetector.hh b/src/ECALdetector.hh index 5354981..6c305ac 100644 --- a/src/ECALdetector.hh +++ b/src/ECALdetector.hh @@ -8,36 +8,36 @@ namespace o2 { - namespace delphes - { +namespace delphes +{ - class ECALdetector - { +class ECALdetector +{ - public: - ECALdetector() = default; - ~ECALdetector() = default; + public: + ECALdetector() = default; + ~ECALdetector() = default; - void setup(float resoEA, float resoEB, float resoEC, float resoPosA, float resoPosB); - bool hasECAL(const Track& track) const; - bool makeSignal(const GenParticle& particle, TLorentzVector &pECAL, float & posZ, float & posPhi); + void setup(float resoEA, float resoEB, float resoEC, float resoPosA, float resoPosB); + bool hasECAL(const Track& track) const; + bool makeSignal(const GenParticle& particle, TLorentzVector& pECAL, float& posZ, float& posPhi); - protected: - Double_t smearPhotonE (const Double_t eTrue); - Double_t sigmaX (const Double_t eTrue); - TLorentzVector smearPhotonP4(const TLorentzVector pTrue); + protected: + Double_t smearPhotonE(const Double_t eTrue); + Double_t sigmaX(const Double_t eTrue); + TLorentzVector smearPhotonP4(const TLorentzVector pTrue); - float mRadius = 120.; // ECAL barrel inner radius [cm] - float mLength = 200.; // ECAL half-length along beam axis [cm] + float mRadius = 120.; // ECAL barrel inner radius [cm] + float mLength = 200.; // ECAL half-length along beam axis [cm] - float mEnergyResolutionA = 0.002; // parameter A of energy resolution in GeV - float mEnergyResolutionB = 0.02; // parameter B of energy resolution in GeV^{1/2} - float mEnergyResolutionC = 0.10; // parameter C of energy resolution - float mPositionResolutionA = 0.15; // parameter A of coordinate resolution in cm - float mPositionResolutionB = 0.30; // parameter B of coordinate resolution in cm*GeV^{1/2} - }; + float mEnergyResolutionA = 0.002; // parameter A of energy resolution in GeV + float mEnergyResolutionB = 0.02; // parameter B of energy resolution in GeV^{1/2} + float mEnergyResolutionC = 0.10; // parameter C of energy resolution + float mPositionResolutionA = 0.15; // parameter A of coordinate resolution in cm + float mPositionResolutionB = 0.30; // parameter B of coordinate resolution in cm*GeV^{1/2} +}; - } // namespace delphes +} // namespace delphes } // namespace o2 #endif /** _DelphesO2_ECALdetector_h_ **/ From c8e40e1b1291146e68463354f56a21e1fc7c35d5 Mon Sep 17 00:00:00 2001 From: Yuri Kharlov Date: Mon, 27 Sep 2021 12:42:35 +0200 Subject: [PATCH 5/5] Fix bug in ECal hit coordinate calculation --- examples/aod/createO2tables.C | 12 ++++++------ examples/aod/createO2tables.h | 2 +- src/ECALdetector.cc | 8 ++++++-- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/examples/aod/createO2tables.C b/examples/aod/createO2tables.C index 21b912a..d967044 100644 --- a/examples/aod/createO2tables.C +++ b/examples/aod/createO2tables.C @@ -249,6 +249,7 @@ int createO2tables(const char* inputFile = "delphes.root", treeReader->ReadEntry(ientry); constexpr float multEtaRange = 2.f; // Range in eta to count the charged particles float dNdEta = 0.f; // Charged particle multiplicity to use in the efficiency evaluation + TLorentzVector pECAL; // 4-momentum of photon in ECAL for (Int_t iparticle = 0; iparticle < particles->GetEntries(); ++iparticle) { // Loop over particles auto particle = (GenParticle*)particles->At(iparticle); @@ -293,10 +294,9 @@ int createO2tables(const char* inputFile = "delphes.root", // info for the ECAL float posZ, posPhi; - TLorentzVector pECAL; if (ecal_detector.makeSignal(*particle, pECAL, posZ, posPhi)) { // to be updated 13.09.2021 - printf("ECAL particle: pid=%d, p=(%g,%g,%g,%g)\n", - particle->PID, particle->Px, particle->Py, particle->Pz, particle->E); + printf("ECAL particle: pid=%d, p=(%g,%g,%g,%g), posZ=%f, posPhi=%f\n", + particle->PID, particle->Px, particle->Py, particle->Pz, particle->E,posZ,posPhi); printf("ECAL p=(%g,%g,%g,%g)\n", pECAL.Px(), pECAL.Py(), pECAL.Pz(), pECAL.E()); ecal.fIndexCollisions = ientry + eventOffset; @@ -304,9 +304,9 @@ int createO2tables(const char* inputFile = "delphes.root", ecal.fPx = pECAL.Px(); ecal.fPy = pECAL.Py(); ecal.fPz = pECAL.Pz(); - ecal.fE = pECAL.E(); - ecal.fPosZ = 0.f; - ecal.fPosPhi = 0.f; + ecal.fE = pECAL.E(); + ecal.fPosZ = posZ; + ecal.fPosPhi = posPhi; FillTree(kA3ECAL); } if constexpr (debug_qa) { diff --git a/examples/aod/createO2tables.h b/examples/aod/createO2tables.h index 0b8d88b..4301c43 100644 --- a/examples/aod/createO2tables.h +++ b/examples/aod/createO2tables.h @@ -486,7 +486,7 @@ struct { Double_t fPx = -1.e10; /// px Double_t fPy = -1.e10; /// py Double_t fPz = -1.e10; /// pz - Double_t fE = -1.e10; /// E + Double_t fE = -1.e10; /// E Float_t fPosZ = -999.f; /// Position in Z Float_t fPosPhi = -999.f; /// Position in phi } ecal; //! structure to keep ECAL info diff --git a/src/ECALdetector.cc b/src/ECALdetector.cc index ffe098d..656878c 100644 --- a/src/ECALdetector.cc +++ b/src/ECALdetector.cc @@ -64,8 +64,12 @@ bool ECALdetector::makeSignal(const GenParticle& particle, Float_t vtY = particle.Y; // particle vertex Y Float_t vtZ = particle.Z; // particle vertex Z - posPhi = particle.Phi; // azimuth angle of a photon hit - posZ = mRadius * particle.CtgTheta; // z-coodrinate of a photon hit + posPhi = p4True.Phi(); // azimuth angle of a photon hit + posZ = -1e6; + Double_t tanTheta = TMath::Tan(p4True.Theta()); + if (tanTheta != 0.) { + posZ = mRadius / tanTheta; // z-coodrinate of a photon hit + } Double_t eSmeared = smearPhotonE(particle.E); // smeared photon energy p4ECAL = smearPhotonP4(p4True);