Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:22:41 

0001 #ifndef MuScleFitUtils_H
0002 #define MuScleFitUtils_H
0003 
0004 /** \class DTHitQualityUtils
0005  *
0006  *  Provide basic functionalities useful for MuScleFit
0007  *
0008  *  \author S. Bolognesi - INFN Torino / T. Dorigo - INFN Padova
0009  * Revised S. Casasso, E. Migliore - UniTo & INFN Torino
0010  */
0011 
0012 #include <CLHEP/Vector/LorentzVector.h>
0013 #include "DataFormats/MuonReco/interface/Muon.h"
0014 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0015 // #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
0016 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0017 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0018 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "TGraphErrors.h"
0021 #include "TH2F.h"
0022 #include "TMinuit.h"
0023 
0024 #include "MuonAnalysis/MomentumScaleCalibration/interface/CrossSectionHandler.h"
0025 #include "MuonAnalysis/MomentumScaleCalibration/interface/BackgroundHandler.h"
0026 #include "MuonAnalysis/MomentumScaleCalibration/interface/ResolutionFunction.h"
0027 #include "MuonAnalysis/MomentumScaleCalibration/interface/MuonPair.h"
0028 #include "MuonAnalysis/MomentumScaleCalibration/interface/GenMuonPair.h"
0029 
0030 #include <vector>
0031 
0032 // #include "Functions.h"
0033 // class biasFunctionBase<std::vector<double> >;
0034 // class scaleFunctionBase<double*>;
0035 template <class T>
0036 class biasFunctionBase;
0037 template <class T>
0038 class scaleFunctionBase;
0039 class smearFunctionBase;
0040 template <class T>
0041 class resolutionFunctionBase;
0042 class backgroundFunctionBase;
0043 class BackgroundHandler;
0044 
0045 class SimTrack;
0046 class TString;
0047 class TTree;
0048 
0049 typedef reco::Particle::LorentzVector lorentzVector;
0050 
0051 class MuScleFitUtils {
0052 public:
0053   // Constructor
0054   // ----------
0055   MuScleFitUtils(){};
0056 
0057   // Destructor
0058   // ----------
0059   virtual ~MuScleFitUtils(){};
0060 
0061   // Operations
0062   // ----------
0063   static std::pair<SimTrack, SimTrack> findBestSimuRes(const std::vector<SimTrack>& simMuons);
0064   static std::pair<MuScleFitMuon, MuScleFitMuon> findBestRecoRes(const std::vector<MuScleFitMuon>& muons);
0065   static std::pair<lorentzVector, lorentzVector> findGenMuFromRes(const reco::GenParticleCollection* genParticles);
0066   static std::pair<lorentzVector, lorentzVector> findGenMuFromRes(const edm::HepMCProduct* evtMC);
0067   static std::pair<lorentzVector, lorentzVector> findSimMuFromRes(const edm::Handle<edm::HepMCProduct>& evtMC,
0068                                                                   const edm::Handle<edm::SimTrackContainer>& simTracks);
0069 
0070   static std::vector<TGraphErrors*> fitMass(TH2F* histo);
0071   static std::vector<TGraphErrors*> fitReso(TH2F* histo);
0072 
0073   static lorentzVector applyScale(const lorentzVector& muon, const std::vector<double>& parval, const int charge);
0074   static lorentzVector applyScale(const lorentzVector& muon, double* parval, const int charge);
0075   static lorentzVector applyBias(const lorentzVector& muon, const int charge);
0076   static lorentzVector applySmearing(const lorentzVector& muon);
0077   static lorentzVector fromPtEtaPhiToPxPyPz(const double* ptEtaPhiE);
0078 
0079   static void minimizeLikelihood();
0080 
0081   static double invDimuonMass(const lorentzVector& mu1, const lorentzVector& mu2);
0082   static double massResolution(const lorentzVector& mu1, const lorentzVector& mu2);
0083   static double massResolution(const lorentzVector& mu1, const lorentzVector& mu2, const std::vector<double>& parval);
0084   static double massResolution(const lorentzVector& mu1, const lorentzVector& mu2, std::unique_ptr<double> parval);
0085   static double massResolution(const lorentzVector& mu1, const lorentzVector& mu2, double* parval);
0086   static double massResolution(const lorentzVector& mu1, const lorentzVector& mu2, const ResolutionFunction& resolFunc);
0087 
0088   static double massProb(const double& mass, const double& rapidity, const int ires, const double& massResol);
0089   /* static double massProb( const double & mass, const double & resEta, const double & rapidity, const double & massResol, const std::vector<double> & parval, const bool doUseBkgrWindow = false ); */
0090   /* static double massProb( const double & mass, const double & resEta, const double & rapidity, const double & massResol, double * parval, const bool doUseBkgrWindow = false ); */
0091   static double massProb(const double& mass,
0092                          const double& resEta,
0093                          const double& rapidity,
0094                          const double& massResol,
0095                          const std::vector<double>& parval,
0096                          const bool doUseBkgrWindow,
0097                          const double& eta1,
0098                          const double& eta2);
0099   static double massProb(const double& mass,
0100                          const double& resEta,
0101                          const double& rapidity,
0102                          const double& massResol,
0103                          double* parval,
0104                          const bool doUseBkgrWindow,
0105                          const double& eta1,
0106                          const double& eta2);
0107   static double computeWeight(const double& mass, const int iev, const bool doUseBkgrWindow = false);
0108 
0109   static double deltaPhi(const double& phi1, const double& phi2) {
0110     double deltaPhi = phi1 - phi2;
0111     while (deltaPhi >= TMath::Pi())
0112       deltaPhi -= 2 * TMath::Pi();
0113     while (deltaPhi < -TMath::Pi())
0114       deltaPhi += 2 * TMath::Pi();
0115     return fabs(deltaPhi);
0116   }
0117   /// Without fabs at the end, used to have a symmetric distribution for the resolution fits and variance computations
0118   static double deltaPhiNoFabs(const double& phi1, const double& phi2) {
0119     double deltaPhi = phi1 - phi2;
0120     while (deltaPhi >= TMath::Pi())
0121       deltaPhi -= 2 * TMath::Pi();
0122     while (deltaPhi < -TMath::Pi())
0123       deltaPhi += 2 * TMath::Pi();
0124     return deltaPhi;
0125   }
0126   static double deltaR(const double& eta1, const double& eta2, const double& phi1, const double& phi2) {
0127     return sqrt(std::pow(eta1 - eta2, 2) + std::pow(deltaPhi(phi1, phi2), 2));
0128   }
0129 
0130   static int debug;      // debug option set by MuScleFit
0131   static bool ResFound;  // bool flag true if best resonance found (cuts on pt and eta)
0132 
0133   static const int totalResNum;             // Total number of resonance: 6
0134   static double massWindowHalfWidth[3][6];  // parameter set by MuScleFitUtils
0135   static double ResGamma[6];                // parameter set by MuScleFitUtils
0136   static double ResMass[6];                 // parameter set by MuScleFitUtils
0137   static double ResMinMass[6];              // parameter set by MuScleFitBase
0138   static double crossSection[6];
0139   static const double mMu2;
0140   static const double muMass;
0141 
0142   // Array of the pdgId of resonances
0143   static const unsigned int motherPdgIdArray[6];
0144 
0145   static unsigned int loopCounter;  // parameter set by MuScleFit
0146 
0147   static int SmearType;
0148   static smearFunctionBase* smearFunction;
0149   static int BiasType;
0150   // No error, we take functions from the same group for scale and bias.
0151   static scaleFunctionBase<std::vector<double> >* biasFunction;
0152   static int ResolFitType;
0153   static resolutionFunctionBase<double*>* resolutionFunction;
0154   static resolutionFunctionBase<std::vector<double> >* resolutionFunctionForVec;
0155   static int ScaleFitType;
0156   static scaleFunctionBase<double*>* scaleFunction;
0157   static scaleFunctionBase<std::vector<double> >* scaleFunctionForVec;
0158   static int BgrFitType;
0159   // Three background regions:
0160   // - one for the Z
0161   // - one for the Upsilons
0162   // - one for J/Psi and Psi2S
0163   static const int backgroundFunctionsRegions;
0164   // static backgroundFunctionBase * backgroundFunctionForRegion[];
0165   // A background function for each resonance
0166   // static backgroundFunctionBase * backgroundFunction[];
0167 
0168   // The Cross section handler takes care of computing the relative cross
0169   // sections to be used depending on the resonances that are being fitted.
0170   // This corresponds to a normalization of the signal pdf.
0171   static CrossSectionHandler* crossSectionHandler;
0172 
0173   // The background handler takes care of using the correct function in each
0174   // window, use regions or resonance windows and rescale the fractions when needed
0175   static BackgroundHandler* backgroundHandler;
0176 
0177   // Parameters used to select whether to do a fit
0178   static std::vector<int> doResolFit;
0179   static std::vector<int> doScaleFit;
0180   static std::vector<int> doCrossSectionFit;
0181   static std::vector<int> doBackgroundFit;
0182 
0183   static int minuitLoop_;
0184   static TH1D* likelihoodInLoop_;
0185   static TH1D* signalProb_;
0186   static TH1D* backgroundProb_;
0187 
0188   static bool duringMinos_;
0189 
0190   static std::vector<double> parSmear;
0191   static std::vector<double> parBias;
0192   static std::vector<double> parResol;
0193   static std::vector<double> parResolStep;
0194   static std::vector<double> parResolMin;
0195   static std::vector<double> parResolMax;
0196   static std::vector<double> parScale;
0197   static std::vector<double> parScaleStep;
0198   static std::vector<double> parScaleMin;
0199   static std::vector<double> parScaleMax;
0200   static std::vector<double> parCrossSection;
0201   static std::vector<double> parBgr;
0202   static std::vector<int> parResolFix;
0203   static std::vector<int> parScaleFix;
0204   static std::vector<int> parCrossSectionFix;
0205   static std::vector<int> parBgrFix;
0206   static std::vector<int> parResolOrder;
0207   static std::vector<int> parScaleOrder;
0208   static std::vector<int> parCrossSectionOrder;
0209   static std::vector<int> parBgrOrder;
0210   static std::vector<int> resfind;
0211   static int FitStrategy;
0212   static bool speedup;                     // parameter set by MuScleFit - whether to speedup processing
0213   static double x[7][10000];               // smearing values set by MuScleFit constructor
0214   static int goodmuon;                     // number of events with a usable resonance
0215   static int counter_resprob;              // number of times there are resolution problems
0216   static double GLZValue[40][1001][1001];  // matrix with integral values of Lorentz * Gaussian
0217   static double GLZNorm[40][1001];         // normalization values per each sigma
0218   static double GLValue[6][1001][1001];    // matrix with integral values of Lorentz * Gaussian
0219   static double GLNorm[6][1001];           // normalization values per each sigma
0220   static double ResMaxSigma[6];            // max sigma of matrix
0221   static double ResHalfWidth[6];           // halfwidth in matrix
0222   static int nbins;                        // number of bins in matrix
0223   static int MuonType;                     // 0, 1, 2 - 0 is GM, 1 is SM, 2 is track
0224   static int MuonTypeForCheckMassWindow;  // Reduced to be 0, 1 or 2. It is = MuonType when MuonType < 3, = 2 otherwise.
0225 
0226   static std::vector<std::vector<double> > parvalue;
0227   // static std::map<unsigned int,std::vector<double> > parvalue;
0228   static std::vector<int> parfix;
0229   static std::vector<int> parorder;
0230 
0231   static std::vector<std::pair<lorentzVector, lorentzVector> > SavedPair;
0232   static std::vector<std::pair<lorentzVector, lorentzVector> > ReducedSavedPair;
0233   static std::vector<std::pair<lorentzVector, lorentzVector> > genPair;
0234   static std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> > SavedPairMuScleFitMuons;
0235   static std::vector<std::pair<MuScleFitMuon, MuScleFitMuon> > genMuscleFitPair;
0236   static std::vector<std::pair<lorentzVector, lorentzVector> > simPair;
0237 
0238   static bool scaleFitNotDone_;
0239 
0240   static bool normalizeLikelihoodByEventNumber_;
0241   // Pointer to the minuit object
0242   static TMinuit* rminPtr_;
0243   // Value stored to check whether to apply a new normalization to the likelihood
0244   static double oldNormalization_;
0245   static unsigned int normalizationChanged_;
0246 
0247   // This must be set to true if using events generated with Sherpa
0248   static bool sherpa_;
0249 
0250   // Decide whether to use the rapidity bins for the Z
0251   static bool rapidityBinsForZ_;
0252 
0253   static int iev_;
0254 
0255   static bool useProbsFile_;
0256 
0257   // Cuts on the muons to use in the fit
0258   static bool separateRanges_;
0259   static double minMuonPt_;
0260   static double maxMuonPt_;
0261   static double minMuonEtaFirstRange_;
0262   static double maxMuonEtaFirstRange_;
0263   static double minMuonEtaSecondRange_;
0264   static double maxMuonEtaSecondRange_;
0265   static double deltaPhiMinCut_;
0266   static double deltaPhiMaxCut_;
0267 
0268   static bool debugMassResol_;
0269   static struct massResolComponentsStruct {
0270     double dmdpt1;
0271     double dmdpt2;
0272     double dmdphi1;
0273     double dmdphi2;
0274     double dmdcotgth1;
0275     double dmdcotgth2;
0276   } massResolComponents;
0277 
0278   // Fit accuracy and debug parameters
0279   static bool startWithSimplex_;
0280   static bool computeMinosErrors_;
0281   static bool minimumShapePlots_;
0282 
0283   /// Method to check if the mass value is within the mass window of the i-th resonance.
0284   // static bool checkMassWindow( const double & mass, const int ires, const double & resMass, const double & leftFactor = 1., const double & rightFactor = 1. );
0285   static bool checkMassWindow(const double& mass, const double& leftBorder, const double& rightBorder);
0286 
0287   /// Computes the probability given the mass, mass resolution and the arrays with the probabilities and the normalizations.
0288   static double probability(const double& mass,
0289                             const double& massResol,
0290                             const double GLvalue[][1001][1001],
0291                             const double GLnorm[][1001],
0292                             const int iRes,
0293                             const int iY);
0294 
0295 protected:
0296 private:
0297   struct byPt {
0298     bool operator()(const reco::Muon& a, const reco::Muon& b) const { return a.pt() > b.pt(); }
0299   };
0300 };
0301 
0302 extern "C" void likelihood(int& npar, double* grad, double& fval, double* xval, int flag);
0303 
0304 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!