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