Heppy/interface/MuScleFitCorrector.h

0001 #ifndef PhysicsTools_Heppy_MuScleFitCorrector_h
0002 #define PhysicsTools_Heppy_MuScleFitCorrector_h
0003
0004 /**
0005  * MuScleFitCorrector class
0006  * Author M. De Mattia - 18/11/2008
0007  * Author S. Casasso   - 25/10/2012
0008  * Author E. Migliore  - 25/10/2012
0009  */
0010
0011 #include <iostream>
0012 #include <fstream>
0013 #include <sstream>
0014 #include "PhysicsTools/Heppy/interface/MuScleFitCorrector_Functions.h"
0015 #include "TLorentzVector.h"
0016 #include "TRandom3.h"
0017
0018 /**
0019  * This is used to have a common set of functions for the specialized templates to use.
0020  * The constructor receives the name identifying the parameters for the correction function.
0021  * It reads the parameters from a txt file in data/.
0022  */
0023
0024 namespace heppy {
0025
0026   class MuScleFitCorrector {
0027   public:
0028     /**
0029    * The constructor takes a string identifying the parameters to read. It
0030    * parses the txt file containing the parameters, extracts the index of the
0031    * correction function and saves the corresponding pointer. It then fills the
0032    * vector of parameters.
0033    */
0034     MuScleFitCorrector(const TString& identifier) {
0035       fileName_ = identifier;
0036       readParameters(identifier);
0037       gRandom_ = new TRandom3();
0038     }
0039     // dummy empty constructor to allow dictionaries
0040     MuScleFitCorrector() {}
0041     ~MuScleFitCorrector() {}
0042
0043     // Returns a pointer to the selected function
0044     scaleFunctionBase<double*>* getScaleFunction() { return scaleFunction_; }
0045     resolutionFunctionBase<double*>* getResolFunction() { return resolutionFunction_; }
0046
0047     // Method to get correct pt
0048     double getCorrectPt(const TLorentzVector& lorentzVector, const int& chg) {
0049       // Loop on all the functions and apply them iteratively on the pt corrected by the previous function.
0050       double pt = lorentzVector.Pt();
0051       pt = (scaleFunction_->scale(pt, lorentzVector.Eta(), lorentzVector.Phi(), chg, scaleParArray_));
0052       return pt;
0053     };
0054
0055     // Return the squared difference of relative pT resolutions data-MC
0056     double getSigmaPtDiffSquared(const double& pt, const double& eta) {
0057       double sigmaPtData = resolutionFunction_->sigmaPt(pt, eta, resolDataParArray_);
0058       double sigmaPtMC = resolutionFunction_->sigmaPt(pt, eta, resolMCParArray_);
0059       return (sigmaPtData * sigmaPtData - sigmaPtMC * sigmaPtMC);
0060     };
0061
0062     // Method to get correct pt (baseline)
0063     double getSmearedPt(const TLorentzVector& lorentzVector, const int& chg, bool fake) {
0064       double pt = lorentzVector.Pt();
0065       double eta = lorentzVector.Eta();
0066       double squaredDiff = getSigmaPtDiffSquared(pt, eta);
0067       if (squaredDiff < 0)
0068         return pt;
0069       double Cfact = 0.;
0070       if (fileName_.Contains("smearPrompt"))
0071         Cfact = 0.85;  // smear Summer12 against 2012 Prompt data
0072       else if (fileName_.Contains("smearReReco"))
0073         Cfact = 0.75;  // smear Summer12 against 2012 ReReco data
0074       else if (fileName_.Contains("2011_MC"))
0075         Cfact = 0.8;  // smear Fall11 against 2011 ReReco data
0076       else {
0077         std::cout << "Are you sure you want to smear data??" << std::endl;
0078         exit(1);
0079       }
0080       double sPar = Cfact * sqrt(squaredDiff);
0081       double curv = ((double)chg / pt);
0082       double normSmearFact = gRandom_->Gaus(0, sPar);
0083       if (fake)
0084         normSmearFact = sPar;
0085       double smearedCurv = curv + fabs(curv) * normSmearFact;
0086
0087       return 1. / ((double)chg * smearedCurv);
0088     };
0089
0090     // Method to apply correction to a TLorentzVector
0091     void applyPtCorrection(TLorentzVector& lorentzVector, const int& chg) {
0092       double eta = lorentzVector.Eta();
0093       double phi = lorentzVector.Phi();
0094       double m = lorentzVector.M();
0095       double corrPt = getCorrectPt(lorentzVector, chg);
0096       lorentzVector.SetPtEtaPhiM(corrPt, eta, phi, m);
0097     };
0098
0099     // Method to apply smearing to a TLorentzVector
0100     void applyPtSmearing(TLorentzVector& lorentzVector, const int& chg, bool fake = false) {
0101       double eta = lorentzVector.Eta();
0102       double phi = lorentzVector.Phi();
0103       double m = lorentzVector.M();
0104       double smearedPt = getSmearedPt(lorentzVector, chg, fake);
0105       lorentzVector.SetPtEtaPhiM(smearedPt, eta, phi, m);
0106     };
0107
0108     std::vector<double> getResolMCParVec() { return resolMCParVec_; }
0109
0110   protected:
0111     // File name
0112     TString fileName_;
0113
0114     // Identification numbers
0115     int scaleFunctionId_;
0116     int resolutionFunctionId_;
0117
0118     // Vectors of parameters
0119     std::vector<double> scaleParVec_;
0120     std::vector<double> resolDataParVec_;
0121     std::vector<double> resolMCParVec_;
0122
0123     // We will use the array for the function calls because it is faster than the vector for random access.
0124     double* scaleParArray_;
0125     double* resolDataParArray_;
0126     double* resolMCParArray_;
0127
0128     // Convert vectors to arrays for faster random access. The first pointer is replaced, thus it is taken by reference.
0129     void convertToArrays() {
0130       int resolParNum = resolutionFunction_->parNum();
0131       int scaleParNum = scaleFunction_->parNum();
0132
0133       int resolDataParVecSize = resolDataParVec_.size();
0134       int resolMCParVecSize = resolMCParVec_.size();
0135       int scaleParVecSize = scaleParVec_.size();
0136
0137       useResol_ = false;
0138       if (resolDataParVecSize != 0 && resolMCParVecSize != 0)
0139         useResol_ = true;
0140
0141       if (resolParNum != resolDataParVecSize && useResol_) {
0142         std::cout << "Error: inconsistent number of parameters: resolutionFunction has " << resolParNum
0143                   << " parameters but " << resolDataParVecSize << " have been read from file" << std::endl;
0144         exit(1);
0145       }
0146
0147       if (resolParNum != resolMCParVecSize && useResol_) {
0148         std::cout << "Error: inconsistent number of parameters: resolutionFunction has " << resolParNum
0149                   << " parameters but " << resolMCParVecSize << " have been read from file" << std::endl;
0150         exit(1);
0151       }
0152
0153       if (scaleParNum != scaleParVecSize) {
0154         std::cout << "Error: inconsistent number of parameters: scaleFunction has " << scaleParNum << " parameters but "
0155                   << scaleParVecSize << " have been read from file" << std::endl;
0156         exit(1);
0157       }
0158
0159       std::vector<double>::const_iterator scaleParIt = scaleParVec_.begin();
0160       std::vector<double>::const_iterator resolDataParIt = resolDataParVec_.begin();
0161       std::vector<double>::const_iterator resolMCParIt = resolMCParVec_.begin();
0162
0163       scaleParArray_ = new double[scaleParNum];
0164       resolDataParArray_ = new double[resolParNum];
0165       resolMCParArray_ = new double[resolParNum];
0166
0167       for (int iPar = 0; iPar < scaleParNum; ++iPar) {
0168         double parameter = *scaleParIt;
0169         scaleParArray_[iPar] = parameter;
0170         ++scaleParIt;
0171       }
0172
0173       if (useResol_) {
0174         for (int iPar = 0; iPar < resolParNum; ++iPar) {
0175           double parameter = *resolDataParIt;
0176           resolDataParArray_[iPar] = parameter;
0177           ++resolDataParIt;
0178         }
0179
0180         for (int iPar = 0; iPar < resolParNum; ++iPar) {
0181           double parameter = *resolMCParIt;
0182           resolMCParArray_[iPar] = parameter;
0183           ++resolMCParIt;
0184         }
0185       }
0186     }
0187     //----------------------------//
0188
0189     // Parser of the parameters file
0190     void readParameters(const TString& fileName) {
0191       scaleParArray_ = 0;
0192       resolDataParArray_ = 0;
0193       resolMCParArray_ = 0;
0194
0195       // Read the parameters file
0196       std::ifstream parametersFile(fileName.Data());
0197
0198       if (!parametersFile.is_open()) {
0199         std::cout << "Error: file " << fileName << " not found. Aborting." << std::endl;
0200         abort();
0201       }
0202
0203       std::string line;
0204       int nReadPar = 0;
0205       std::string iteration("Iteration ");
0206       // Loop on the file lines
0207       while (parametersFile) {
0208         getline(parametersFile, line);
0209         size_t lineInt = line.find("value");
0210         size_t iterationSubStr = line.find(iteration);
0211
0212         if (iterationSubStr != std::string::npos) {
0213           std::stringstream sLine(line);
0214           std::string num;
0215           int scaleFunctionNum = 0;
0216           int resolutionFunctionNum = 0;
0217           int wordCounter = 0;
0218
0219           // Warning: this strongly depends on the parameters file structure.
0220           while (sLine >> num) {
0221             ++wordCounter;
0222             if (wordCounter == 8) {
0223               std::stringstream in(num);
0224               in >> resolutionFunctionNum;
0225             }
0226
0227             if (wordCounter == 9) {
0228               std::stringstream in(num);
0229               in >> scaleFunctionNum;
0230             }
0231           }
0232
0233           scaleFunctionId_ = scaleFunctionNum;
0234           scaleFunction_ = scaleFunctionService(scaleFunctionNum);
0235           resolutionFunctionId_ = resolutionFunctionNum;
0236           resolutionFunction_ = resolutionFunctionService(resolutionFunctionNum);
0237
0238           /*       std::cout<<"Function IDs: "<<std::endl; */
0239           /*       std::cout<<"     scale function number "<<scaleFunctionId_<<std::endl; */
0240           /*       std::cout<<"     resolution function number "<<resolutionFunctionId_<<std::endl; */
0241         }
0242
0243         int nScalePar = scaleFunction_->parNum();
0244         int nResolPar = resolutionFunction_->parNum();
0245
0246         if ((lineInt != std::string::npos)) {
0247           size_t subStr1 = line.find("value");
0248           std::stringstream paramStr;
0249           double param = 0;
0250           paramStr << line.substr(subStr1 + 5);
0251           paramStr >> param;
0252           // Fill the last vector of parameters, which corresponds to this iteration.
0253           if (nReadPar < nScalePar) {
0254             scaleParVec_.push_back(param);
0255           } else if (nReadPar < (nScalePar + nResolPar)) {
0256             resolDataParVec_.push_back(param);
0257           } else if (nReadPar < (nScalePar + 2 * nResolPar)) {
0258             resolMCParVec_.push_back(param);
0259           }
0260           ++nReadPar;
0261         }
0262       }
0263       convertToArrays();
0264
0265       /*   std::cout<<"Scale function n. "<<scaleFunctionId_<<" has "<<scaleFunction_->parNum()<<"parameters:"<<std::endl; */
0266       /*   for (int ii=0; ii<scaleFunction_->parNum(); ++ii){ */
0267       /*     std::cout<<"par["<<ii<<"] = "<<scaleParArray_[ii]<<std::endl; */
0268       /*   } */
0269
0270       if (useResol_) {
0271         /*     std::cout<<"Resolution data function n. "<<resolutionFunctionId_<<" has "<<resolutionFunction_->parNum()<<"parameters:"<<std::endl; */
0272         /*     for (int ii=0; ii<resolutionFunction_->parNum(); ++ii){ */
0273         /*       std::cout<<"par["<<ii<<"] = "<<resolDataParArray_[ii]<<std::endl; */
0274         /*     } */
0275
0276         /*     std::cout<<"Resolution MC function n. "<<resolutionFunctionId_<<" has "<<resolutionFunction_->parNum()<<"parameters:"<<std::endl; */
0277         /*     for (int ii=0; ii<resolutionFunction_->parNum(); ++ii){ */
0278         /*       std::cout<<"par["<<ii<<"] = "<<resolMCParArray_[ii]<<std::endl; */
0279         /*     } */
0280       }
0281     }
0282
0283     // Functions
0284     scaleFunctionBase<double*>* scaleFunction_;
0285     resolutionFunctionBase<double*>* resolutionFunction_;
0286
0287     // Pointer for TRandom3 access
0288     TRandom3* gRandom_;
0289
0290     // Bool for using resolution function or not (value depends from the information on the parameters txt file)
0291     bool useResol_;
0292   };
0293 }  // namespace heppy
0294
0295 #endif