00001 #ifndef RoccoR_h
00002 #define RoccoR_h
00003
00004 #include <iostream>
00005 #include <fstream>
00006 #include "TMath.h"
00007 #include <cmath>
00008
00009 struct CrystalBall{
00010 private:
00011 static constexpr double pi = M_PI;
00012 static constexpr double SPiO2 = sqrt(M_PI/2);
00013 static constexpr double S2 = sqrt(2);
00014
00015 public:
00016 double m;
00017 double s;
00018 double a;
00019 double n;
00020
00021 double B;
00022 double C;
00023 double D;
00024 double N;
00025
00026 double NA;
00027 double Ns;
00028 double NC;
00029 double F;
00030 double G;
00031 double k;
00032
00033 double cdfMa;
00034 double cdfPa;
00035
00036 CrystalBall(double m_=0, double s_=1, double a_=10, double n_=10){
00037 init(m_, s_, a_, n_);
00038 }
00039
00040 void init(double m_, double s_, double a_, double n_){
00041 m=m_;
00042 s=s_;
00043 a=a_;
00044 n=n_;
00045
00046 double fa = fabs(a);
00047 double expa = exp(-fa*fa/2);
00048 double A = pow(n/fa, n)*expa;
00049 double C1 = n/fa/(n-1)*expa;
00050 double D1 = 2*SPiO2*erf(fa/S2);
00051
00052 B = n/fa-fa;
00053 C = (D1+2*C1)/C1;
00054 D = (D1+2*C1)/2;
00055
00056 N = 1.0/s/(D1+2*C1);
00057 k = 1.0/(n-1);
00058
00059 NA = N*A;
00060 Ns = N*s;
00061 NC = Ns*C1;
00062 F = 1-fa*fa/n;
00063 G = s*n/fa;
00064
00065 cdfMa=cdf(m-a*s);
00066 cdfPa=cdf(m+a*s);
00067 }
00068
00069 double pdf(double x) const{
00070 double d=(x-m)/s;
00071 if(d<-a) return NA*pow(B-d, -n);
00072 if(d> a) return NA*pow(B+d, -n);
00073 return N*exp(-d*d/2);
00074 }
00075
00076 double pdf(double x, double ks, double dm) const{
00077 double d=(x-m-dm)/(s*ks);
00078 if(d<-a) return NA/ks*pow(B-d, -n);
00079 if(d> a) return NA/ks*pow(B+d, -n);
00080 return N/ks*exp(-d*d/2);
00081 }
00082
00083 double cdf(double x) const{
00084 double d = (x-m)/s;
00085 if(d<-a) return NC / pow(F-s*d/G, n-1);
00086 if(d> a) return NC * (C - pow(F+s*d/G, 1-n) );
00087 return Ns*(D-SPiO2*erf(-d/S2));
00088 }
00089
00090 double invcdf(double u) const{
00091 if(u<cdfMa) return m + G*(F - pow(NC/u, k) );
00092 if(u>cdfPa) return m - G*(F - pow(C-u/NC, -k) );
00093 return m - S2*s*TMath::ErfInverse((D - u/Ns ) / SPiO2);
00094 }
00095 };
00096
00097
00098
00099
00100
00101 class RocRes{
00102 private:
00103 static const int NMAXETA=12;
00104 static const int NMAXTRK=12;
00105
00106 int NETA;
00107 int NTRK;
00108 int NMIN;
00109
00110 double BETA[NMAXETA+1];
00111 double ntrk[NMAXETA][NMAXTRK+1];
00112 double dtrk[NMAXETA][NMAXTRK+1];
00113
00114 double width[NMAXETA][NMAXTRK];
00115 double alpha[NMAXETA][NMAXTRK];
00116 double power[NMAXETA][NMAXTRK];
00117
00118 double rmsA[NMAXETA][NMAXTRK];
00119 double rmsB[NMAXETA][NMAXTRK];
00120 double rmsC[NMAXETA][NMAXTRK];
00121
00122 double kDat[NMAXETA];
00123 double kRes[NMAXETA];
00124
00125 int getBin(double x, const int NN, const double *b) const;
00126
00127
00128 public:
00129 enum TYPE {MC, Data, Extra};
00130
00131 CrystalBall cb[NMAXETA][NMAXTRK];
00132
00133 RocRes();
00134 int getEtaBin(double feta) const;
00135 int getNBinDT(double v, int H) const;
00136 int getNBinMC(double v, int H) const;
00137 double getUrnd(int H, int F, double v) const;
00138 void dumpParams();
00139 void init(std::string filename);
00140
00141 void reset();
00142
00143 ~RocRes(){}
00144
00145 double Sigma(double pt, int H, int F) const;
00146 double kSpread(double gpt, double rpt, double eta, int nlayers, double w) const;
00147 double kSmear(double pt, double eta, TYPE type, double v, double u) const;
00148 double kSmear(double pt, double eta, TYPE type, double v, double u, int n) const;
00149 double kExtra(double pt, double eta, int nlayers, double u, double w) const;
00150 double getkDat(int H) const{return kDat[H];}
00151 double getkRes(int H) const{return kRes[H];}
00152 };
00153
00154
00155 class RocOne{
00156 private:
00157 static const int NMAXETA=22;
00158 static const int NMAXPHI=16;
00159 static const double MPHI;
00160
00161 int NETA;
00162 int NPHI;
00163
00164 double BETA[NMAXETA+1];
00165 double DPHI;
00166
00167 double M[2][NMAXETA][NMAXPHI];
00168 double A[2][NMAXETA][NMAXPHI];
00169 double D[2][NMAXETA];
00170
00171 RocRes RR;
00172
00173 int getBin(double x, const int NN, const double *b) const;
00174 int getBin(double x, const int nmax, const double xmin, const double dx) const;
00175
00176 public:
00177 enum TYPE{MC, DT};
00178
00179 RocOne();
00180 ~RocOne(){}
00181
00182 RocOne(std::string filename, int iTYPE=0, int iSYS=0, int iMEM=0);
00183 bool checkSYS(int iSYS, int iMEM, int kSYS=0, int kMEM=0);
00184 bool checkTIGHT(int iTYPE, int iSYS, int iMEM, int kTYPE=0, int kSYS=0, int kMEM=0);
00185 void reset();
00186 void init(std::string filename, int iTYPE=0, int iSYS=0, int iMEM=0);
00187
00188 double kScaleDT(int Q, double pt, double eta, double phi) const;
00189 double kScaleMC(int Q, double pt, double eta, double phi, double kSMR=1) const;
00190 double kScaleAndSmearMC(int Q, double pt, double eta, double phi, int n, double u, double w) const;
00191 double kScaleFromGenMC(int Q, double pt, double eta, double phi, int n, double gt, double w) const;
00192 double kGenSmear(double pt, double eta, double v, double u, RocRes::TYPE TT=RocRes::Data) const;
00193
00194 double getM(int T, int H, int F) const{return M[T][H][F];}
00195 double getA(int T, int H, int F) const{return A[T][H][F];}
00196 double getK(int T, int H) const{return T==DT?RR.getkDat(H):RR.getkRes(H);}
00197 RocRes& getR() {return RR;}
00198 };
00199
00200
00201 class RoccoR{
00202 public:
00203 RoccoR();
00204 RoccoR(std::string dirname);
00205 ~RoccoR();
00206
00207 void init(std::string dirname);
00208
00209 double kGenSmear(double pt, double eta, double v, double u, RocRes::TYPE TT=RocRes::Data, int s=0, int m=0) const;
00210 double kScaleDT(int Q, double pt, double eta, double phi, int s=0, int m=0) const;
00211
00212 double kScaleAndSmearMC(int Q, double pt, double eta, double phi, int n, double u, double w, int s=0, int m=0) const;
00213 double kScaleFromGenMC(int Q, double pt, double eta, double phi, int n, double gt, double w, int s=0, int m=0) const;
00214
00215
00216 double getM(int T, int H, int F, int E=0, int m=0) const{return RC[E][m].getM(T,H,F);}
00217 double getA(int T, int H, int F, int E=0, int m=0) const{return RC[E][m].getA(T,H,F);}
00218 double getK(int T, int H, int E=0, int m=0) const{return RC[E][m].getK(T,H);}
00219
00220 int Nset() const{return RC.size();}
00221 int Nmem(int s=0) const{return RC[s].size();}
00222
00223 private:
00224 std::vector<std::vector<RocOne> > RC;
00225 };
00226
00227 #endif
00228
