diff --git a/delphes/include/DelphesConfig.h b/delphes/include/DelphesConfig.h index ae6104f86533619afa41e12c91d19ab53e2dab09..e4e65fd84ef883facd4ef25d1aa62dba4cb93ee5 100644 --- a/delphes/include/DelphesConfig.h +++ b/delphes/include/DelphesConfig.h @@ -86,10 +86,10 @@ class EtaRange { double GetMax( void ) const { return m_Max; }; unsigned GetMomentumRangeCount( void ) const { return m_MomentumRanges.size(); } - const MomentumRange *FirstRange( void ) const { + const MomentumRange *FirstMomentumRange( void ) const { return m_MomentumRanges.empty() ? 0 : m_MomentumRanges.front(); }; - const MomentumRange *LastRange( void ) const { + const MomentumRange *LastMomentumRange( void ) const { return m_MomentumRanges.empty() ? 0 : m_MomentumRanges.back(); }; @@ -127,6 +127,7 @@ class DelphesConfig { bool StoreSigmaEntry(MomentumRange *mrange, int pdg, double sigma); + void UsePtMode( void ) { m_PtMode = true; }; void AddZeroSigmaEntries( void ); void Print(); int Check(); @@ -160,6 +161,8 @@ class DelphesConfig { // Yes, this is a global parameter, so that all entries in the output table have // a consistent meaning; bool m_EfficiencyContaminationMode; + + bool m_PtMode; }; #endif diff --git a/delphes/include/DelphesConfigTOF.h b/delphes/include/DelphesConfigTOF.h index 096eb50e13007ab4e48ec5081c8099e943627afc..9c59ffcf6dd8cef1ef085cd2b29c87f106e4b2ac 100644 --- a/delphes/include/DelphesConfigTOF.h +++ b/delphes/include/DelphesConfigTOF.h @@ -9,13 +9,20 @@ class DelphesConfigTOF: public DelphesConfig { DelphesConfigTOF(const char *dname): DelphesConfig(dname), m_T0Resolution(0.0), m_DetectorResolution(0.0), m_InstallationRadius(0.0), m_EtaMin(0.0), m_EtaMax(0.0), m_MomentumMin(0.0), m_MomentumMax(0.0), - m_EtaBinCount(0), m_MomentumBinCount(0), m_MagneticField(0.0) {}; + m_EtaBinCount(0), m_MomentumBinCount(0), m_MagneticField(0.0), + m_MomentumResolutionA(0.0), m_MomentumResolutionB(0.0), m_PathLengthResolution(0.0) {}; ~DelphesConfigTOF() {}; - void SetT0Resolution (double value) { m_T0Resolution = value; } - void SetDetectorResolution(double value) { m_DetectorResolution = value; } - void SetInstallationRadius(double value) { m_InstallationRadius = value; } - void SetMagneticField (double value) { m_MagneticField = value; } + void SetT0Resolution (double value) { m_T0Resolution = value; } + // Assume a constant one; [mm]; + void SetPathLengthResolution(double value) { m_PathLengthResolution = value; } + // dp/p ~ a*p + b; either in momentum or in Pt; + void SetMomentumResolution(double a, double b) { + m_MomentumResolutionA = a; m_MomentumResolutionB = b; + }; + void SetDetectorResolution (double value) { m_DetectorResolution = value; } + void SetInstallationRadius (double value) { m_InstallationRadius = value; } + void SetMagneticField (double value) { m_MagneticField = value; } // Well, would be more natural to give Z-range along the beam line; void SetEtaRange(double min, double max, unsigned ebins) { m_EtaMin = min; @@ -26,7 +33,7 @@ class DelphesConfigTOF: public DelphesConfig { m_MomentumMin = min; m_MomentumMax = max; m_MomentumBinCount = ebins; - }; + }; // TOF-specific call (input parameters -> sigma counts); int DoSigmaCalculations( void ); @@ -40,6 +47,8 @@ class DelphesConfigTOF: public DelphesConfig { double m_MomentumMin, m_MomentumMax; unsigned m_EtaBinCount, m_MomentumBinCount; double m_MagneticField; + + double m_MomentumResolutionA, m_MomentumResolutionB, m_PathLengthResolution; }; #endif diff --git a/delphes/scripts/delphes-btof.C b/delphes/scripts/delphes-btof.C index c33ed79e551fe31f7f172ecd9c1b59b7dd95b791..c23e641af556753360cdde0e90cdce4653ce2b4f 100644 --- a/delphes/scripts/delphes-btof.C +++ b/delphes/scripts/delphes-btof.C @@ -6,6 +6,7 @@ void delphes_btof( void ) { auto btof = new DelphesConfigTOF("BTOF"); + btof->UsePtMode(); // Define particle mass hypotheses in ascending mass order; yes, there is no // reason to overcomplicate things; @@ -15,16 +16,24 @@ void delphes_btof( void ) btof->AddMassHypothesis("proton"); // Define t0 and detector time resolution is [ns]; - btof->SetT0Resolution (0.030); - btof->SetDetectorResolution(0.020); + btof->SetT0Resolution (0.030); + btof->SetDetectorResolution (0.020); + btof->SetMomentumResolution (0.000, 0.000); + // [mm] here; + btof->SetPathLengthResolution(1.000); // Installation radius in [m]; constant magnetic field in [T]; - btof->SetInstallationRadius(0.500); - btof->SetMagneticField (3.000); + btof->SetInstallationRadius (0.500); + btof->SetMagneticField (3.000); + + // eta and momentum range and binning; + btof->SetEtaRange(-1.0, 1.0, 10); + // Do not mind to use Pt rather than 1/Pt bins; [GeV/c]; + btof->SetMomentumRange(1.0, 2.0, 10); // This input is sufficient to allocate the internal tables and calculate // time of flight for various mass hypotheses; - DoSigmaCalculations(); + btof->DoSigmaCalculations(); // This is again some generic stuff; btof->AddZeroSigmaEntries(); diff --git a/delphes/source/DelphesConfig.cc b/delphes/source/DelphesConfig.cc index e6ba6fef9c5e096029972dbbf57933d072f75286..1645068e9490ca7730377294bbe056d0dccba1a0 100644 --- a/delphes/source/DelphesConfig.cc +++ b/delphes/source/DelphesConfig.cc @@ -10,8 +10,9 @@ DelphesConfig::DelphesConfig(const char *dname): m_Name(dname), m_EtaMin(0.0), m_EtaMax(0.0), - m_MomentumMin(0.0), m_MomentumMax(0.0), m_EfficiencyContaminationMode(false)//, + m_MomentumMin(0.0), m_MomentumMax(0.0), m_EfficiencyContaminationMode(false), //m_PionThreshold(0.0), m_KaonThreshold(0.0), m_ProtonThreshold(0.0) + m_PtMode(false) { m_DatabasePDG = new TDatabasePDG(); } // DelphesConfig::DelphesConfig() @@ -186,15 +187,15 @@ int DelphesConfig::Check( void ) // First eta range: assign momentum range; if (erange == m_EtaRanges.front()) { - m_MomentumMin = erange->FirstRange()->GetMin(); - m_MomentumMax = erange->LastRange ()->GetMax(); + m_MomentumMin = erange->FirstMomentumRange()->GetMin(); + m_MomentumMax = erange->LastMomentumRange ()->GetMax(); continue; } //if // Subsequent eta ranges: check that momentum range is the same; - if (erange->FirstRange()->GetMin() != m_MomentumMin || - erange->LastRange() ->GetMax() != m_MomentumMax) + if (erange->FirstMomentumRange()->GetMin() != m_MomentumMin || + erange->LastMomentumRange() ->GetMax() != m_MomentumMax) _ERROR_("Full momentum range should be the same in all eta ranges!"); for(auto mrange: erange->m_MomentumRanges) @@ -210,7 +211,7 @@ int DelphesConfig::Check( void ) void DelphesConfig::WriteMassHypothesis(FILE *fout, unsigned ih) { unsigned dim = m_MassHypotheses.size(); - const char *pstr = "pt * cosh(eta)"; + const char *pstr = m_PtMode ? "pt" : "pt * cosh(eta)"; auto prev = ih ? m_MassHypotheses[ih-1] : 0; auto hypo = m_MassHypotheses[ih ]; diff --git a/delphes/source/DelphesConfigTOF.cc b/delphes/source/DelphesConfigTOF.cc index 0769e72a725259eed73d1610fb24dbddd138841a..30b3d8a450216f8800296b1aada47ac1cb80440f 100644 --- a/delphes/source/DelphesConfigTOF.cc +++ b/delphes/source/DelphesConfigTOF.cc @@ -8,6 +8,7 @@ int DelphesConfigTOF::DoSigmaCalculations( void ) // // FIXME: sanity check on the input parameters is needed; // + assert(m_PtMode); double etaStep = (m_EtaMax - m_EtaMin ) / m_EtaBinCount; double pStep = (m_MomentumMax - m_MomentumMin) / m_MomentumBinCount; @@ -16,13 +17,19 @@ int DelphesConfigTOF::DoSigmaCalculations( void ) for(unsigned ie=0; ie<m_EtaBinCount; ie++) { double eta = etaStep*(ie - (m_EtaBinCount-1)/2.); - auto erange = AddEtaRange(eta - etaStep/2, eta + etaStep/2); + // Assign eta range; avoid machine accuracy issues when checking continuity; + auto erange = AddEtaRange(ie ? m_EtaRanges.back()->GetMax() : eta - etaStep/2, eta + etaStep/2); // Loop though all of the requested momentum bins; for(unsigned ip=0; ip<m_MomentumBinCount; ip++) { - double p = pStep*(ip - (m_MomentumBinCount-1)/2.); + // So momentum has a mening of Pt in this code; see assert() above; + double pt = pStep*(ip - (m_MomentumBinCount-1)/2.); - auto mrange = erange->GetMomentumRange(p - pStep/2, p + pStep/2); + // Assign momentum range; avoid machine accuracy issues when checking continuity; + auto mrange = erange->GetMomentumRange(ip ? erange->LastMomentumRange()->GetMax() : + pt - pStep/2, pt + pStep/2); + // mnemonic rule: 30 MeV/c in 1kGs field -> curvature radius of 1m; convert to [mm]; + double r = 1000*(0.1/m_MagneticField)*(pt/0.03); // Loop through all of the mass hypotheses and populate entries one by one; for(unsigned ih=0; ih<m_MassHypotheses.size(); ih++) {