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