delphes/include/DelphesConfig.h

0001
0002 #include <string>
0003
0004 #include <TObject.h>
0005 #include <TVector3.h>
0006
0007 class TDatabasePDG;
0008 #include <TParticlePDG.h>
0009 #include <TRandom.h>
0010
0011 #ifndef _DELPHES_CONFIG_
0012 #define _DELPHES_CONFIG_
0013
0014 class MassHypothesis: public TObject {
0015  public:
0016   MassHypothesis(TParticlePDG *pdg = 0, double max_contamination_left = 0.0,
0017          double max_contamination_right = 0.0): m_Diff(0.0), m_ChiSquare(0.0), m_PDG(pdg),
0018     m_MaxContaminationLeft(max_contamination_left),
0019     m_MaxContaminationRight(max_contamination_right), m_Threshold(0.0) {};
0020   ~ MassHypothesis() {};
0021
0022   double Mass( void ) const { return m_PDG->Mass(); };
0023   int PdgCode( void ) const { return m_PDG->PdgCode(); };
0024
0025   void SetThreshold(double value) { m_Threshold = value; };
0026   double GetThreshold( void ) const { return m_Threshold; };
0027
0028   double GetMeasurementOffset( void ) const { return m_Diff; };
0029
0030   double m_Diff;       //!
0031   double m_ChiSquare;  //!
0032
0033  private:
0034   TParticlePDG *m_PDG;
0035
0036   double m_MaxContaminationLeft, m_MaxContaminationRight;
0037
0038   double m_Threshold;
0039
0040   ClassDef(MassHypothesis, 1)
0041 };
0042
0043 class MomentumRange: public TObject {
0044  public:
0045   enum range {undefined, below_pion_threshold, below_kaon_threshold, below_proton_threshold};
0046
0047  MomentumRange(double min = 0.0, double max = 0.0): m_Min(min), m_Max(max), m_MatrixDim(0), m_Matrix(0),
0048     m_Range(MomentumRange::undefined) {
0049     if (min > max) std::swap(m_Min, m_Max);
0050   };
0051
0052   // Just fill the buffer; FIXME: eventually will mess up with these two calls;
0053   void AddSigmaValue(double value) { m_SigmaValues.push_back(value); };
0054   void AddSigmaValue(double measurement, double sigma) {
0055     m_MeasurementValues.push_back(measurement);
0056     m_SigmaValues.push_back(sigma);
0057   };
0058
0059   double GetMin( void )             const { return m_Min; };
0060   double GetMax( void )             const { return m_Max; };
0061   unsigned GetSigmaCount( void )    const { return m_SigmaValues.size(); }
0062   double GetSigma(unsigned ih) const {
0063     return ih < m_SigmaValues.size() ? m_SigmaValues[ih] : 0.0;
0064   };
0065   double GetMeasurement(unsigned ih) const {
0066     return ih < m_MeasurementValues.size() ? m_MeasurementValues[ih] : 0.0;
0067   };
0068
0069   unsigned m_MatrixDim;
0070   double *m_Matrix; //[m_MatrixDim]
0071
0072  private:
0073   double m_Min, m_Max;
0074
0075   std::vector<double> m_MeasurementValues;
0076   std::vector<double> m_SigmaValues;
0077   MomentumRange::range m_Range;
0078
0079   ClassDef(MomentumRange, 2)
0080 };
0081
0082 class EtaRange: public TObject {
0083  public:
0084  EtaRange(double min = 0.0, double max = 0.0): m_Min(min), m_Max(max) {
0085     if (min > max) std::swap(m_Min, m_Max);
0086   };
0087
0088   // A hack-like method, to populate PDG entries one by one;
0089   MomentumRange *GetMomentumRange(double min, double max);
0090
0091   // This will be the userland AddMomentumRange() wrapper;
0092   template<typename... Args> MomentumRange *AddMomentumRange(double min, double max, Args... args) {
0093     MomentumRange *prev = m_MomentumRanges.size() ? m_MomentumRanges.back() : 0;
0094     if (prev && min != prev->GetMax()) {
0095       printf("Momentum ranges sequence issue: %6.2f is followed by %6.2f\n",
0096          prev->GetMax(), min);
0097       exit(1);
0098     } //for if
0099     auto mrange = new MomentumRange(min, max);
0100     m_MomentumRanges.push_back(mrange);
0101
0102     ArgumentRecursion(mrange, args...);
0103
0104     return mrange;
0105   };
0106
0107   double GetMin( void ) const { return m_Min; };
0108   double GetMax( void ) const { return m_Max; };
0109
0110   unsigned GetMomentumRangeCount( void ) const { return m_MomentumRanges.size(); }
0111   const MomentumRange *FirstMomentumRange( void ) const {
0112     return m_MomentumRanges.empty() ? 0 : m_MomentumRanges.front();
0113   };
0114   const MomentumRange *LastMomentumRange( void ) const {
0115     return m_MomentumRanges.empty() ? 0 : m_MomentumRanges.back();
0116   };
0117
0118   std::vector<MomentumRange*> m_MomentumRanges;
0119
0120   MomentumRange *GetMomentumRange(double p) const {
0121     for(auto mrange: m_MomentumRanges)
0122       // Prefer to use '<=' here (help max value?);
0123       if (mrange->GetMin() <= p && p <= mrange->GetMax())
0124     return mrange;
0125
0126     return 0;
0127   };
0128
0129  private:
0130   template <typename T, typename... Args>
0131     void ArgumentRecursion(MomentumRange *mrange, const T& firstArg, const Args&... args) {
0132     // NB: c++11 static_assert() here guarantees that analysis.C does not *compile*
0133     // rather than fail at run-time;
0134     static_assert(std::is_same<T, double>::value, "'const char *' parameter(s) expected!");
0135
0136     mrange->AddSigmaValue(firstArg);
0137
0138     // Recursively process 'fnum' and the rest of detector name arguments;
0139     ArgumentRecursion(mrange, args...);
0140   };
0141   // Once no more arguments left, call this dummy;
0142   void ArgumentRecursion(MomentumRange *mrange) {};
0143
0144   double m_Min, m_Max;
0145
0146   ClassDef(EtaRange, 1)
0147 };
0148
0149 class DelphesConfig: public TObject {
0150  public:
0151   DelphesConfig( void );
0152   DelphesConfig(const char *dname);
0153   virtual ~DelphesConfig() {};
0154
0155   MassHypothesis *AddMassHypothesis(int pdg, double max_contamination_left = 1.0,
0156                     double max_contamination_right = 1.0);
0157   MassHypothesis *AddMassHypothesis(const char *pname, double max_contamination_left = 1.0,
0158                     double max_contamination_right = 1.0);
0159
0160   EtaRange *AddEtaRange(double min, double max);
0161
0162   bool StoreSigmaEntry(MomentumRange *mrange, int pdg, double sigma);
0163   bool StoreSigmaEntry(MomentumRange *mrange, int pdg, double measurement, double sigma);
0164
0165   void UsePtMode( void ) { m_PtMode = true; };
0166   void AddZeroSigmaEntries( void );
0167   void Print();
0168   int  Check(bool rigorous = true);
0169   virtual int  Calculate()               = 0;
0170   //virtual bool ApplyThresholdModeLogic() = 0;
0171   void WriteTcl(bool check = true);
0172
0173   unsigned GetMassHypothesisCount( void ) const { return m_MassHypotheses.size(); };
0174   void GetMassHypotheses(int harr[]) const {
0175     unsigned id = 0;
0176
0177     for(auto hypo: m_MassHypotheses)
0178       harr[id++] = hypo->PdgCode();
0179   };
0180   MassHypothesis *GetMassHypothesis(int pdg, bool ignore_sign = true);
0181
0182   double GetEtaMin( void ) {
0183     if (!m_EtaRanges.size()) return -99;
0184
0185     return m_EtaRanges[0]->GetMin();
0186   };
0187   double GetEtaMax( void ) {
0188     if (!m_EtaRanges.size()) return -99;
0189
0190     return m_EtaRanges[m_EtaRanges.size()-1]->GetMax();
0191   };
0192
0193   EtaRange *GetEtaRange(double eta) const {
0194     for(auto erange: m_EtaRanges)
0195       // Prefer to use '<=' here (help max value?);
0196       if (erange->GetMin() <= eta && eta <= erange->GetMax())
0197     return erange;
0198
0199     return 0;
0200   };
0201
0202   int GetSmearingMatrix(double eta, double p, double hmtx[]) const {
0203     for(unsigned iq=0; iq<m_MassHypotheses.size()*m_MassHypotheses.size(); iq++)
0204       hmtx[iq] = 0.0;
0205
0206     auto erange = GetEtaRange(eta);
0207     if (!erange) return -1;
0208
0209     auto mrange = erange->GetMomentumRange(p);
0210     if (!mrange) return -2;
0211
0212     for(unsigned iq=0; iq<mrange->m_MatrixDim; iq++)
0213       hmtx[iq] = mrange->m_Matrix[iq];
0214
0215     return 0;
0216   };
0217   MomentumRange *GetEtaMomentumRange(double eta, double p) const {
0218     auto erange = GetEtaRange(eta);
0219     if (!erange) return 0;
0220
0221     return erange->GetMomentumRange(p);
0222   };
0223
0224   int GetSmearingMatrix(const TVector3 p, double hmtx[]) const {
0225     return GetSmearingMatrix(p.Eta(), p.Mag(), hmtx);
0226   };
0227
0228  protected:
0229   std::string m_Name;
0230
0231   TDatabasePDG *m_DatabasePDG;
0232
0233   // Eta ranges in ascending order;
0234   std::vector<EtaRange*> m_EtaRanges;
0235
0236   // Assume all masses are different and provided in ascending order;
0237   std::vector<MassHypothesis*> m_MassHypotheses;
0238
0239   // FIXME: need a std::set, eventually?;
0240   unsigned GetHypoIndex(MassHypothesis *hypo) {
0241     for(unsigned ih=0; ih<m_MassHypotheses.size(); ih++)
0242       if (m_MassHypotheses[ih] == hypo)
0243     return ih;
0244
0245     return 0;
0246   };
0247
0248  private:
0249   MassHypothesis *AddMassHypothesisCore(TParticlePDG *pdg, double max_contamination_left,
0250                     double max_contamination_right);
0251
0252   void DetermineThresholds( void );
0253   void WriteMassHypothesis(FILE *fout, unsigned ih);
0254
0255   double m_EtaMin, m_EtaMax;
0256   double m_MomentumMin, m_MomentumMax;
0257   // These only make sense for pi/K/p case;
0258   //double m_PionThreshold, m_KaonThreshold, m_ProtonThreshold;
0259
0260  protected:
0261   // Yes, this is a global parameter, so that all entries in the output table have
0262   // a consistent meaning;
0263   bool m_EfficiencyContaminationMode;
0264
0265   bool m_PtMode;
0266
0267   // FIXME: need seed setting method;
0268   TRandom m_rndm; //!
0269
0270   ClassDef(DelphesConfig, 1)
0271 };
0272
0273 #endif