diff --git a/delphes/include/DelphesConfig.h b/delphes/include/DelphesConfig.h
index e4e65fd84ef883facd4ef25d1aa62dba4cb93ee5..f441addf9596390c5ca55270455b8f182e94f13c 100644
--- a/delphes/include/DelphesConfig.h
+++ b/delphes/include/DelphesConfig.h
@@ -38,8 +38,12 @@ class MomentumRange {
if (min > max) std::swap(m_Min, m_Max);
};
- // Just fill the buffer;
+ // Just fill the buffer; FIXME: eventually will mess up with these two calls;
void AddSigmaValue(double value) { m_SigmaValues.push_back(value); };
+ void AddSigmaValue(double measurement, double sigma) {
+ m_MeasurementValues.push_back(measurement);
+ m_SigmaValues.push_back(sigma);
+ };
double GetMin( void ) const { return m_Min; };
double GetMax( void ) const { return m_Max; };
@@ -53,6 +57,7 @@ class MomentumRange {
private:
double m_Min, m_Max;
+ std::vector<double> m_MeasurementValues;
std::vector<double> m_SigmaValues;
MomentumRange::range m_Range;
};
@@ -126,6 +131,7 @@ class DelphesConfig {
EtaRange *AddEtaRange(double min, double max);
bool StoreSigmaEntry(MomentumRange *mrange, int pdg, double sigma);
+ bool StoreSigmaEntry(MomentumRange *mrange, int pdg, double measurement, double sigma);
void UsePtMode( void ) { m_PtMode = true; };
void AddZeroSigmaEntries( void );
diff --git a/delphes/include/DelphesConfigTOF.h b/delphes/include/DelphesConfigTOF.h
index 9c59ffcf6dd8cef1ef085cd2b29c87f106e4b2ac..340cba93709c5c101506266dfce4ee2431fcb567 100644
--- a/delphes/include/DelphesConfigTOF.h
+++ b/delphes/include/DelphesConfigTOF.h
@@ -49,6 +49,8 @@ class DelphesConfigTOF: public DelphesConfig {
double m_MagneticField;
double m_MomentumResolutionA, m_MomentumResolutionB, m_PathLengthResolution;
+
+ double tof(double m, double p, double l);
};
#endif
diff --git a/delphes/scripts/delphes-btof.C b/delphes/scripts/delphes-btof.C
index c23e641af556753360cdde0e90cdce4653ce2b4f..ea0aeb383b6d945355f8ab983d6e741b46705999 100644
--- a/delphes/scripts/delphes-btof.C
+++ b/delphes/scripts/delphes-btof.C
@@ -1,4 +1,6 @@
+//
+// Units are [mm], [GeV], [T], [ps];
//
// root -l delphes_btof.C
//
@@ -12,32 +14,32 @@ void delphes_btof( void )
// reason to overcomplicate things;
//btof->AddMassHypothesis(-11);
btof->AddMassHypothesis("pi+");
- btof->AddMassHypothesis("K+");
- btof->AddMassHypothesis("proton");
+ //btof->AddMassHypothesis("K+");
+ //btof->AddMassHypothesis("proton");
- // Define t0 and detector time resolution is [ns];
- btof->SetT0Resolution (0.030);
- btof->SetDetectorResolution (0.020);
+ // Define t0 and detector time resolution is [ps];
+ btof->SetT0Resolution (30.00);
+ btof->SetDetectorResolution (30.00);
btof->SetMomentumResolution (0.000, 0.000);
- // [mm] here;
+ // Units are [mm] throughout the code;
btof->SetPathLengthResolution(1.000);
- // Installation radius in [m]; constant magnetic field in [T];
- btof->SetInstallationRadius (0.500);
+ // Installation radius in [mm]; constant magnetic field in [T];
+ btof->SetInstallationRadius (500.0);
btof->SetMagneticField (3.000);
// eta and momentum range and binning;
- btof->SetEtaRange(-1.0, 1.0, 10);
+ btof->SetEtaRange (0.0, 0.0, 1);
// Do not mind to use Pt rather than 1/Pt bins; [GeV/c];
- btof->SetMomentumRange(1.0, 2.0, 10);
+ btof->SetMomentumRange(0.1, 1.1, 20);
// This input is sufficient to allocate the internal tables and calculate
// time of flight for various mass hypotheses;
btof->DoSigmaCalculations();
// This is again some generic stuff;
- btof->AddZeroSigmaEntries();
- btof->Print();
+ //btof->AddZeroSigmaEntries();
+ //btof->Print();
//btof->Write();
exit(0);
} // delphes_btof()
diff --git a/delphes/source/DelphesConfig.cc b/delphes/source/DelphesConfig.cc
index 1645068e9490ca7730377294bbe056d0dccba1a0..01b520b6f35b65b7edae2dc32d49dbef885af429 100644
--- a/delphes/source/DelphesConfig.cc
+++ b/delphes/source/DelphesConfig.cc
@@ -80,7 +80,8 @@ MomentumRange *EtaRange::GetMomentumRange(double min, double max)
// -------------------------------------------------------------------------------------
-bool DelphesConfig::StoreSigmaEntry(MomentumRange *mrange, int pdg, double sigma)
+bool DelphesConfig::StoreSigmaEntry(MomentumRange *mrange, int pdg, double measurement,
+ double sigma)
{
// Check that at least 'pdg' code is in order;
unsigned mdim = mrange->GetSigmaCount();
@@ -90,13 +91,20 @@ bool DelphesConfig::StoreSigmaEntry(MomentumRange *mrange, int pdg, double sigma
if (pdg != m_MassHypotheses[mdim]->PdgCode()) return false;
- mrange->AddSigmaValue(sigma);
+ mrange->AddSigmaValue(measurement, sigma);
return true;
} // DelphesConfig::StoreSigmaEntry()
// -------------------------------------------------------------------------------------
+bool DelphesConfig::StoreSigmaEntry(MomentumRange *mrange, int pdg, double sigma)
+{
+ return StoreSigmaEntry(mrange, pdg, 0.0, sigma);
+} // DelphesConfig::StoreSigmaEntry()
+
+// -------------------------------------------------------------------------------------
+
void DelphesConfig::AddZeroSigmaEntries( void )
{
for(auto erange: m_EtaRanges)
diff --git a/delphes/source/DelphesConfigTOF.cc b/delphes/source/DelphesConfigTOF.cc
index 30b3d8a450216f8800296b1aada47ac1cb80440f..6760435c6590bed4333fef6d231ffee5dd9bc3bf 100644
--- a/delphes/source/DelphesConfigTOF.cc
+++ b/delphes/source/DelphesConfigTOF.cc
@@ -1,6 +1,18 @@
#include <DelphesConfigTOF.h>
+// Speed of light; [mm/ps];
+#define _SoL_ ( 0.29979)
+
+// -------------------------------------------------------------------------------------
+
+double DelphesConfigTOF::tof(double m, double p, double l)
+{
+ double e = sqrt(p*p + m*m), beta = p/e;//, c = 299.79, t = l / (beta*_SoL_);
+
+ return l / (beta*_SoL_);
+} // DelphesConfigTOF::tof()
+
// -------------------------------------------------------------------------------------
int DelphesConfigTOF::DoSigmaCalculations( void )
@@ -15,32 +27,52 @@ int DelphesConfigTOF::DoSigmaCalculations( void )
// Loop though all of the requested eta bins;
for(unsigned ie=0; ie<m_EtaBinCount; ie++) {
- double eta = etaStep*(ie - (m_EtaBinCount-1)/2.);
+ double eta = m_EtaMin + etaStep*(ie + 0.5);
+ double theta = 2*std::atan(exp(-eta));
// 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++) {
- // So momentum has a mening of Pt in this code; see assert() above;
- double pt = pStep*(ip - (m_MomentumBinCount-1)/2.);
+ // So momentum has a meaning of Pt in this code; see assert() above;
+ double pt = m_MomentumMin + pStep*(ip +0.5), p = pt / sin(theta);
// 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);
+ // Particle will not reach the barrel layer; FIXME: may want to give few cm up?;
+ if (2*r < m_InstallationRadius) continue;
// Loop through all of the mass hypotheses and populate entries one by one;
for(unsigned ih=0; ih<m_MassHypotheses.size(); ih++) {
auto hypo = m_MassHypotheses[ih];
- //mass[ih] = m_MassHypotheses[ih]->Mass();
-
- bool ret = StoreSigmaEntry(mrange, hypo->PdgCode(), 0.030);
- if (!ret) {
- printf("Failed to insert a sigma entry\n");
- exit(0);
- } //if
+ double m = m_MassHypotheses[ih]->Mass();
+
+ {
+ double alfa = 2*asin(m_InstallationRadius/(2.0*r));
+ double lxy = r*alfa, l = lxy/sin(theta);
+ // Use speed of light in [mm/ps] units;
+ double e = sqrt(p*p + m*m), beta = p/e;//, c = 0.29979;
+
+ double sp = 0.01, t = tof(m, p, l);
+ double p1 = p - 0.1*sp, p2 = p + 0.1*p, dt = tof(m, p2, l) - tof(m, p1, l);
+
+ printf("%3d %3d: eta= %5.2f th=%5.1f pt= %5.2f, r= %7.2f, m= %5.3f alfa=%5.1f, beta= %5.3f, l = %5.1f t = %7.4f\n",
+ ie, ip, eta, theta*180/M_PI, pt, r, m, alfa*180/M_PI, beta, l, t);
+
+ // Now calculate the uncertainty;
+ double s1 = m_T0Resolution, s2 = m_DetectorResolution, s3 = m_PathLengthResolution / _SoL_, s4 = sp*fabs(dt)/(0.2*sp);
+ printf(" -> %7.2f %7.2f %7.2f %7.2f\n", s1, s2, s3, s4);
+
+ bool ret = StoreSigmaEntry(mrange, hypo->PdgCode(), sqrt(s1*s1+s2*s2+s3*s3+s4*s4));
+ if (!ret) {
+ printf("Failed to insert a sigma entry\n");
+ exit(0);
+ } //if
+ }
} //for ih
} //for ip
} //for ie