Reference Manual

Go to the documentation of this file.
00001 #ifndef MuScleFitUtils_H
00002 #define MuScleFitUtils_H
00003 
00013 #include <CLHEP/Vector/LorentzVector.h>
00014 #include "DataFormats/MuonReco/interface/Muon.h"
00015 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00016 // #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
00017 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00018 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
00019 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00021 #include "TGraphErrors.h"
00022 #include "TH2F.h"
00023 #include "TMinuit.h"
00024 
00025 #include "MuonAnalysis/MomentumScaleCalibration/interface/CrossSectionHandler.h"
00026 #include "MuonAnalysis/MomentumScaleCalibration/interface/BackgroundHandler.h"
00027 #include "MuonAnalysis/MomentumScaleCalibration/interface/ResolutionFunction.h"
00028 
00029 #include <vector>
00030 
00031 // #include "Functions.h"
00032 // class biasFunctionBase<std::vector<double> >;
00033 // class scaleFunctionBase<double*>;
00034 template <class T> class biasFunctionBase;
00035 template <class T> class scaleFunctionBase;
00036 class smearFunctionBase;
00037 template <class T> class resolutionFunctionBase;
00038 class backgroundFunctionBase;
00039 class BackgroundHandler;
00040 
00041 class SimTrack;
00042 class TString;
00043 class TTree;
00044 
00045 typedef reco::Particle::LorentzVector lorentzVector;
00046 
00047 class MuScleFitUtils {
00048 public:
00049 
00050   // Constructor
00051   // ----------
00052   MuScleFitUtils() {};
00053 
00054   // Destructor
00055   // ----------
00056   virtual ~MuScleFitUtils() {};
00057 
00058   // Operations
00059   // ----------
00060   static std::pair<SimTrack, SimTrack> findBestSimuRes( const std::vector<SimTrack>& simMuons );
00061   static std::pair<lorentzVector, lorentzVector> findBestRecoRes( const std::vector<reco::LeafCandidate>& muons );
00062   static std::pair <lorentzVector, lorentzVector> findGenMuFromRes( const reco::GenParticleCollection* genParticles);
00063   static std::pair<lorentzVector, lorentzVector> findGenMuFromRes( const edm::HepMCProduct* evtMC );
00064   static std::pair<lorentzVector, lorentzVector> findSimMuFromRes( const edm::Handle<edm::HepMCProduct> & evtMC,
00065                                                                    const edm::Handle<edm::SimTrackContainer> & simTracks);
00066 
00067   static std::vector<TGraphErrors*> fitMass (TH2F* histo);
00068   static std::vector<TGraphErrors*> fitReso (TH2F* histo);
00069 
00070   static lorentzVector applyScale( const lorentzVector & muon, const std::vector<double> & parval, const int charge );
00071   static lorentzVector applyScale( const lorentzVector & muon, double* parval, const int charge );
00072   static lorentzVector applyBias( const lorentzVector & muon, const int charge );
00073   static lorentzVector applySmearing( const lorentzVector & muon );
00074   static lorentzVector fromPtEtaPhiToPxPyPz( const double* ptEtaPhiE );
00075 
00076   static void minimizeLikelihood();
00077 
00078   static double invDimuonMass( const lorentzVector & mu1, const lorentzVector & mu2 );
00079   static double massResolution( const lorentzVector & mu1, const lorentzVector & mu2 );
00080   static double massResolution( const lorentzVector & mu1, const lorentzVector & mu2, const std::vector<double> & parval );
00081   static double massResolution( const lorentzVector & mu1, const lorentzVector & mu2, std::auto_ptr<double> parval );
00082   static double massResolution( const lorentzVector & mu1, const lorentzVector & mu2, double* parval );
00083   static double massResolution( const lorentzVector& mu1, const lorentzVector& mu2, const ResolutionFunction & resolFunc );
00084 
00085   static double massProb( const double & mass, const double & rapidity, const int ires, const double & massResol );
00086   /* static double massProb( const double & mass, const double & resEta, const double & rapidity, const double & massResol, const std::vector<double> & parval, const bool doUseBkgrWindow = false ); */
00087   /* static double massProb( const double & mass, const double & resEta, const double & rapidity, const double & massResol, double * parval, const bool doUseBkgrWindow = false ); */
00088   static double massProb( const double & mass, const double & resEta, const double & rapidity, const double & massResol, const std::vector<double> & parval, const bool doUseBkgrWindow, const double & eta1, const double & eta2 );
00089   static double massProb( const double & mass, const double & resEta, const double & rapidity, const double & massResol, double * parval, const bool doUseBkgrWindow, const double & eta1, const double & eta2 );
00090   static double computeWeight( const double & mass, const int iev, const bool doUseBkgrWindow = false );
00091 
00092   static double deltaPhi( const double & phi1, const double & phi2 )
00093   {
00094     double deltaPhi = phi1 - phi2;
00095     while(deltaPhi >= TMath::Pi()) deltaPhi -= 2*TMath::Pi();
00096     while(deltaPhi < -TMath::Pi()) deltaPhi += 2*TMath::Pi();
00097     return fabs(deltaPhi);
00098   }
00100   static double deltaPhiNoFabs( const double & phi1, const double & phi2 )
00101   {
00102     double deltaPhi = phi1 - phi2;
00103     while(deltaPhi >= TMath::Pi()) deltaPhi -= 2*TMath::Pi();
00104     while(deltaPhi < -TMath::Pi()) deltaPhi += 2*TMath::Pi();
00105     return deltaPhi;
00106   }
00107   static double deltaR(const double & eta1, const double & eta2, const double & phi1, const double & phi2)
00108   {
00109     return sqrt( std::pow( eta1-eta2, 2 ) + std::pow( deltaPhi(phi1, phi2), 2 ) );
00110   }
00111 
00112   static int debug;       // debug option set by MuScleFit
00113   static bool ResFound;   // bool flag true if best resonance found (cuts on pt and eta)
00114 
00115   static const int totalResNum; // Total number of resonance: 6
00116   static double massWindowHalfWidth[3][6]; // parameter set by MuScleFitUtils
00117   static double ResGamma[6];     // parameter set by MuScleFitUtils
00118   static double ResMass[6];      // parameter set by MuScleFitUtils
00119   static double ResMinMass[6];      // parameter set by MuScleFitBase
00120   static double crossSection[6];
00121   static const double mMu2;
00122   static const double muMass;
00123 
00124   // Array of the pdgId of resonances
00125   static const unsigned int motherPdgIdArray[6];
00126 
00127   static unsigned int loopCounter; // parameter set by MuScleFit
00128 
00129   static int SmearType;
00130   static smearFunctionBase * smearFunction;
00131   static int BiasType;
00132   // No error, we take functions from the same group for scale and bias.
00133   static scaleFunctionBase<std::vector<double> > * biasFunction;
00134   static int ResolFitType;
00135   static resolutionFunctionBase<double *> * resolutionFunction;
00136   static resolutionFunctionBase<std::vector<double> > * resolutionFunctionForVec;
00137   static int ScaleFitType;
00138   static scaleFunctionBase<double*> * scaleFunction;
00139   static scaleFunctionBase<std::vector<double> > * scaleFunctionForVec;
00140   static int BgrFitType;
00141   // Three background regions:
00142   // - one for the Z
00143   // - one for the Upsilons
00144   // - one for J/Psi and Psi2S
00145   static const int backgroundFunctionsRegions;
00146   // static backgroundFunctionBase * backgroundFunctionForRegion[];
00147   // A background function for each resonance
00148   // static backgroundFunctionBase * backgroundFunction[];
00149 
00150   // The Cross section handler takes care of computing the relative cross
00151   // sections to be used depending on the resonances that are being fitted.
00152   // This corresponds to a normalization of the signal pdf.
00153   static CrossSectionHandler * crossSectionHandler;
00154 
00155   // The background handler takes care of using the correct function in each
00156   // window, use regions or resonance windows and rescale the fractions when needed
00157   static BackgroundHandler * backgroundHandler;
00158 
00159   // Parameters used to select whether to do a fit
00160   static std::vector<int> doResolFit;
00161   static std::vector<int> doScaleFit;
00162   static std::vector<int> doCrossSectionFit;
00163   static std::vector<int> doBackgroundFit;
00164 
00165   static int minuitLoop_;
00166   static TH1D* likelihoodInLoop_;
00167   static TH1D* signalProb_;
00168   static TH1D* backgroundProb_;
00169 
00170   static bool duringMinos_;
00171 
00172   static std::vector<double> parSmear;
00173   static std::vector<double> parBias;
00174   static std::vector<double> parResol;
00175   static std::vector<double> parResolStep;
00176   static std::vector<double> parResolMin;
00177   static std::vector<double> parResolMax;
00178   static std::vector<double> parScale;
00179   static std::vector<double> parScaleStep;
00180   static std::vector<double> parScaleMin;
00181   static std::vector<double> parScaleMax;
00182   static std::vector<double> parCrossSection;
00183   static std::vector<double> parBgr;
00184   static std::vector<int> parResolFix;
00185   static std::vector<int> parScaleFix;
00186   static std::vector<int> parCrossSectionFix;
00187   static std::vector<int> parBgrFix;
00188   static std::vector<int> parResolOrder;
00189   static std::vector<int> parScaleOrder;
00190   static std::vector<int> parCrossSectionOrder;
00191   static std::vector<int> parBgrOrder;
00192   static std::vector<int> resfind;
00193   static int FitStrategy;
00194   static bool speedup;       // parameter set by MuScleFit - whether to speedup processing
00195   static double x[7][10000]; // smearing values set by MuScleFit constructor
00196   static int goodmuon;       // number of events with a usable resonance
00197   static int counter_resprob;// number of times there are resolution problems
00198   static double GLZValue[40][1001][1001]; // matrix with integral values of Lorentz * Gaussian
00199   static double GLZNorm[40][1001];        // normalization values per each sigma
00200   static double GLValue[6][1001][1001]; // matrix with integral values of Lorentz * Gaussian
00201   static double GLNorm[6][1001];        // normalization values per each sigma
00202   static double ResMaxSigma[6];         // max sigma of matrix
00203   static double ResHalfWidth[6];        // halfwidth in matrix
00204   static int nbins;                     // number of bins in matrix
00205   static int MuonType; // 0, 1, 2 - 0 is GM, 1 is SM, 2 is track
00206   static int MuonTypeForCheckMassWindow; // Reduced to be 0, 1 or 2. It is = MuonType when MuonType < 3, = 2 otherwise.
00207 
00208   static std::vector<std::vector<double> > parvalue;
00209   // static std::map<unsigned int,std::vector<double> > parvalue;
00210   static std::vector<int> parfix;
00211   static std::vector<int> parorder;
00212 
00213   static std::vector<std::pair<lorentzVector,lorentzVector> > SavedPair;
00214   static std::vector<std::pair<lorentzVector,lorentzVector> > ReducedSavedPair;
00215   static std::vector<std::pair<lorentzVector,lorentzVector> > genPair;
00216   static std::vector<std::pair<lorentzVector,lorentzVector> > simPair;
00217 
00218   static bool scaleFitNotDone_;
00219 
00220   static bool normalizeLikelihoodByEventNumber_;
00221   // Pointer to the minuit object
00222   static TMinuit * rminPtr_;
00223   // Value stored to check whether to apply a new normalization to the likelihood
00224   static double oldNormalization_;
00225   static unsigned int normalizationChanged_;
00226 
00227   // This must be set to true if using events generated with Sherpa
00228   static bool sherpa_;
00229 
00230   // Decide whether to use the rapidity bins for the Z
00231   static bool rapidityBinsForZ_;
00232 
00233   static int iev_;
00234 
00235   static bool useProbsFile_;
00236 
00237   // Cuts on the muons to use in the fit
00238   static bool separateRanges_;
00239   static double minMuonPt_;
00240   static double maxMuonPt_;
00241   static double minMuonEtaFirstRange_;
00242   static double maxMuonEtaFirstRange_;
00243   static double minMuonEtaSecondRange_;
00244   static double maxMuonEtaSecondRange_;
00245   static double deltaPhiMinCut_;
00246   static double deltaPhiMaxCut_;
00247 
00248   static bool debugMassResol_;
00249   static struct massResolComponentsStruct
00250   {
00251     double dmdpt1;
00252     double dmdpt2;
00253     double dmdphi1;
00254     double dmdphi2;
00255     double dmdcotgth1;
00256     double dmdcotgth2;
00257   } massResolComponents;
00258 
00259   // Fit accuracy and debug parameters
00260   static bool startWithSimplex_;
00261   static bool computeMinosErrors_;
00262   static bool minimumShapePlots_;
00263 
00265   // static bool checkMassWindow( const double & mass, const int ires, const double & resMass, const double & leftFactor = 1., const double & rightFactor = 1. );
00266   static bool checkMassWindow( const double & mass, const double & leftBorder, const double & rightBorder );
00267 
00269   static double probability( const double & mass, const double & massResol,
00270                              const double GLvalue[][1001][1001], const double GLnorm[][1001],
00271                              const int iRes, const int iY );
00272 
00273 protected:
00274 
00275 private:
00276 
00277    struct byPt {
00278      bool operator() (const reco::Muon &a, const reco::Muon &b) const {
00279        return a.pt() > b.pt();
00280      }
00281    };
00282 
00283 };
00284 
00285 extern "C" void likelihood (int& npar, double* grad, double& fval, double* xval, int flag);
00286 
00287 #endif
/data/refman/pasoursint/CMSSW_5_3_10_patch2/src/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFitUtils.h
