Second commit of BTOF Delphes files

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

delphes/include/DelphesConfigTOF.h

@@ -9,10 +9,17 @@ 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; }
+  // 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; }
@@ -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

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;
@@ -17,14 +18,22 @@ void delphes_btof( void )
   // Define t0 and detector time resolution is [ns];
   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);
 
+  // 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();

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  ];

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.);
-    
-      auto mrange = erange->GetMomentumRange(p - pStep/2, p + pStep/2);
+      // So momentum has a mening of Pt in this code; see assert() above;
+      double pt = pStep*(ip - (m_MomentumBinCount-1)/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++) {