From f5500461a1ff746b768f2e1008b96afb78f6c5da Mon Sep 17 00:00:00 2001 From: Marco De Mattia Date: Mon, 5 Oct 2009 11:48:40 +0000 Subject: [PATCH] --- yaml --- r: 74719 b: "refs/heads/CMSSW_7_1_X" c: 7437c0c5a85104af36f281189eb38fa410d5881b h: "refs/heads/CMSSW_7_1_X" i: 74717: 9b6c8806fbd024115a3ea8043d04d4ed514443bb 74715: 44e759c55f8e1ffd49f8f0458e9383bd0366378e 74711: ddfe56e3ac1291d9a155d66fcaaa321f3dd1da52 74703: 4012db231ec23e0e23b23296113dd8038263f7b5 74687: 583ddbf69d3c8c7f43ffff29a17c824f54f4bbce v: v3 --- [refs] | 2 +- .../interface/Functions.h | 98 +++++---------- .../interface/MuScleFitBase.h | 5 + .../interface/MuScleFitUtils.h | 19 ++- .../plugins/MuScleFit.cc | 56 ++++++--- .../src/MuScleFitBase.cc | 43 +++++-- .../src/MuScleFitPlotter.cc | 60 ++++++--- .../src/MuScleFitUtils.cc | 115 +++++++++++------- 8 files changed, 240 insertions(+), 158 deletions(-) diff --git a/[refs] b/[refs] index 41d68365e1c28..2e7ac0576307e 100644 --- a/[refs] +++ b/[refs] @@ -1,3 +1,3 @@ --- refs/heads/gh-pages: 09c786f70121f131b3715aaf3464996502bbeb7e -"refs/heads/CMSSW_7_1_X": d4aa8aecaf0b3374aa37a33991d54d223e9df430 +"refs/heads/CMSSW_7_1_X": 7437c0c5a85104af36f281189eb38fa410d5881b diff --git a/trunk/MuonAnalysis/MomentumScaleCalibration/interface/Functions.h b/trunk/MuonAnalysis/MomentumScaleCalibration/interface/Functions.h index 5f6d05b17961b..0c0efa011e9a7 100644 --- a/trunk/MuonAnalysis/MomentumScaleCalibration/interface/Functions.h +++ b/trunk/MuonAnalysis/MomentumScaleCalibration/interface/Functions.h @@ -1121,16 +1121,17 @@ class resolutionFunctionType11 : public resolutionFunctionBase { }; /** - * Same as type12 but with free parameters for transition region and center of second parabola. + * Same as type11 but with free parameters for transition region and center of second parabola. * It also imposes continuity of the two fuctions. * Adds also two additional parameters to allow a linear and a quadratic dependence from pt (the * resolution vs Pt has been seen to grow with Pt for misaligned samples. + * Also replaced the sigmaCotgTh and sigmaPhi with those from type8. */ // Resolution Type 12 template class resolutionFunctionType12 : public resolutionFunctionBase { public: - resolutionFunctionType12() { this->parNum_ = 11; } + resolutionFunctionType12() { this->parNum_ = 19; } // linear in pt and by points in eta virtual double sigmaPt(const double & pt, const double & eta, const T & parval) { double fabsEta = fabs(eta); @@ -1148,32 +1149,51 @@ class resolutionFunctionType12 : public resolutionFunctionBase { } // 1/pt in pt and quadratic in eta virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) { - return( 0.004 ); + return( parval[10]+parval[11]/pt + parval[12]*fabs(eta)+parval[13]*eta*eta ); } // 1/pt in pt and quadratic in eta virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) { - return( 0.001 ); + return( parval[14]+parval[15]/pt + parval[16]*fabs(eta)+parval[17]*eta*eta ); } +// // 1/pt in pt and quadratic in eta +// virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) { +// return( 0.004 ); +// } +// // 1/pt in pt and quadratic in eta +// virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) { +// return( 0.001 ); +// } virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const vector & parResolOrder, const int muonType) { double thisStep[] = { 0.001, 0.00001, 0.0000001, 0.00000001, 0.00000001, 0.00000001, 0.00000001, 0.00000001, - 0.001, 0.0001, 0.000001 }; + 0.001, 0.0001, 0.000001, + 0.00002, 0.0002, 0.0000002, 0.00002, + 0.00002, 0.0002, 0.00000002, 0.000002 }; TString thisParName[] = { "etaTransition", "offsetEtaHigh", "coeffOverPt", "coeffHighPt", "linaerEtaCentral", "parabEtaCentral", "linaerEtaHigh", "parabEtaHigh", - "secondParabolaCenter", "linearPt", "quadraticPt" }; - double thisMini[] = { -1.1, -1.1, -0.1, -0.1, + "secondParabolaCenter", "linearPt", "quadraticPt", + "Cth res. sc.", "Cth res. 1/Pt sc.", "Cth res. Eta sc.", "Cth res. Eta^2 sc.", + "Phi res. sc.", "Phi res. 1/Pt sc.", "Phi res. Eta sc.", "Phi res. Eta^2 sc." }; + double thisMini[] = { -1.1, -1.1, -1.1, -1.1, 0.0001, 0.0005, 0.0005, 0.001, - -1.0, -1.0, -1.0 }; + -1.0, -1.0, -1.0, + -0.001, 0.002, -0.0001, -0.0001, + -0.0001, 0.0005, -0.0001, -0.00001 }; if( muonType == 1 ) { double thisMaxi[] = { 1., 1., 1., 1., 1., 1., 1., 1., - 1., 1. ,1. }; + 1., 1., 1., + 1., 1., 1., 1., + 1., 1., 1., 1. }; + this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName ); } else { - double thisMaxi[] = { 1.8, -0.8, -0.001, -0.001, + double thisMaxi[] = { 1.8, -0.8, 0.1, 0.1, 0.005, 0.05, 0.05, 0.05, - 2.4, 2.0, 2.0 }; + 2.4, 2.0, 2.0, + 0.001, 0.005, 0.00004, 0.0007, + 0.001, 0.01, -0.0000015, 0.0004 }; this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName ); } } @@ -1235,62 +1255,6 @@ class resolutionFunctionType12 : public resolutionFunctionBase { // } // }; -/// This resolution function uses the same function as type11, but one for Pt < 4GeV and one for Pt >= 4 GeV. -// Resolution Type 13 -template -class resolutionFunctionType13 : public resolutionFunctionBase { - public: - resolutionFunctionType13() { this->parNum_ = 16; } - // linear in pt and by points in eta - virtual double sigmaPt(const double & pt, const double & eta, const T & parval) { - double fabsEta = fabs(eta); - int shift = 0; - if( pt >= 4 ) { - shift = 8; - } - if(fabsEta<1.2) - return (parval[0+shift]+ parval[2+shift]*1./pt + pt/(pt+parval[3+shift]) + parval[4+shift]*fabsEta + parval[5+shift]*eta*eta); - else - return (parval[1+shift]+ parval[2+shift]*1./pt + pt/(pt+parval[3+shift]) + parval[6+shift]*fabs((fabsEta-1.6)) + parval[7+shift]*(fabsEta-1.6)*(fabsEta-1.6)); - } - // 1/pt in pt and quadratic in eta - virtual double sigmaCotgTh(const double & pt, const double & eta, const T & parval) { - return( 0.004 ); - } - // 1/pt in pt and quadratic in eta - virtual double sigmaPhi(const double & pt, const double & eta, const T & parval) { - return( 0.001 ); - } - virtual void setParameters(double* Start, double* Step, double* Mini, double* Maxi, int* ind, TString* parname, const T & parResol, const vector & parResolOrder, const int muonType) { - - double thisStep[] = { 0.00001, 0.00001, 0.0000001, 0.00000001, - 0.00000001, 0.00000001, 0.00000001, 0.00001, - 0.00001, 0.00001, 0.0000001, 0.00000001, - 0.00000001, 0.00000001, 0.00000001, 0.00001 }; - TString thisParName[] = { "offsetEtaCentral", "offsetEtaHigh", "coeffOverPt", "coeffHighPt", - "linaerEtaCentral", "parabEtaCentral", "linaerEtaHigh", "parabEtaHigh", - "offsetEtaCentral_2", "offsetEtaHigh_2", "coeffOverPt_2", "coeffHighPt_2", - "linaerEtaCentral_2", "parabEtaCentral_2", "linaerEtaHigh_2", "parabEtaHigh_2" }; - double thisMini[] = { -1.1, -1.1, -0.1, -0.1, - 0.0001, 0.0005, 0.0005, 0.001, - -1.1, -1.1, -0.1, -0.1, - 0.0001, 0.0005, 0.0005, 0.001 }; - if( muonType == 1 ) { - double thisMaxi[] = { 1., 1., 1., 1., - 1., 1., 1., 1., - 1., 1., 1., 1., - 1., 1., 1., 1. }; - this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName ); - } else { - double thisMaxi[] = { -0.8, -0.8, -0.001, -0.001, - 0.005, 0.05, 0.05, 0.05, - -0.8, -0.8, -0.001, -0.001, - 0.005, 0.05, 0.05, 0.05 }; - this->setPar( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, thisStep, thisMini, thisMaxi, thisParName ); - } - } -}; - // ------------ ATTENTION ----------- // // Other functions are not in for now // // ---------------------------------- // diff --git a/trunk/MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitBase.h b/trunk/MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitBase.h index 7d98196731e7b..2f236dbd3fba8 100644 --- a/trunk/MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitBase.h +++ b/trunk/MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitBase.h @@ -19,6 +19,8 @@ class MuScleFitBase { public: MuScleFitBase(const edm::ParameterSet& iConfig) : + probabilitiesFileInPath_( iConfig.getUntrackedParameter( "ProbabilitiesFileInPath" , "MuonAnalysis/MomentumScaleCalibration/test/Probs_new_Horace_CTEQ_1000.root" ) ), + probabilitiesFile_( iConfig.getUntrackedParameter( "ProbabilitiesFile" , "" ) ), theMuonType_( iConfig.getParameter( "MuonType" ) ), theMuonLabel_( iConfig.getParameter( "MuonLabel" ) ), theRootFileName_( iConfig.getUntrackedParameter("OutputFileName") ), @@ -46,6 +48,9 @@ class MuScleFitBase /// Raed probability distributions from a local root file. void readProbabilityDistributionsFromFile(); + string probabilitiesFileInPath_; + string probabilitiesFile_; + int theMuonType_; edm::InputTag theMuonLabel_; string theRootFileName_; diff --git a/trunk/MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitUtils.h b/trunk/MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitUtils.h index 559f3eba5c6d6..59871915ec332 100644 --- a/trunk/MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitUtils.h +++ b/trunk/MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitUtils.h @@ -5,8 +5,8 @@ * * Provide basic functionalities useful for MuScleFit * - * $Date: 2009/08/07 15:57:06 $ - * $Revision: 1.9 $ + * $Date: 2009/08/14 12:10:16 $ + * $Revision: 1.10 $ * \author S. Bolognesi - INFN Torino / T. Dorigo - INFN Padova */ @@ -205,6 +205,21 @@ class MuScleFitUtils { static int iev_; + // Cuts on the muons to use in the fit + static double minMuonPt_; + static double maxMuonEta_; + + static bool debugMassResol_; + static struct massResolComponentsStruct + { + double dmdpt1; + double dmdpt2; + double dmdphi1; + double dmdphi2; + double dmdcotgth1; + double dmdcotgth2; + } massResolComponents; + /// Method to check if the mass value is within the mass window of the i-th resonance. static bool checkMassWindow( const double & mass, const int ires, const double & resMass, const double & leftFactor = 1., const double & rightFactor = 1. ); diff --git a/trunk/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFit.cc b/trunk/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFit.cc index faab881881c51..df4240e575b7f 100644 --- a/trunk/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFit.cc +++ b/trunk/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFit.cc @@ -1,8 +1,8 @@ // \class MuScleFit // Fitter of momentum scale and resolution from resonance decays to muon track pairs // -// $Date: 2009/09/08 09:52:17 $ -// $Revision: 1.56 $ +// $Date: 2009/09/24 07:47:28 $ +// $Revision: 1.57 $ // \author R. Bellan, C.Mariotti, S.Bolognesi - INFN Torino / T.Dorigo, M.De Mattia - INFN Padova // // Recent additions: @@ -91,10 +91,13 @@ // system. All functions in file MuScleFitUtils.cc have been suitably changed. // // ---------------------------------------------------------------------------------- -// Modifications by M.De Mattia 13/3/2009 -// -------------------------------------- +// Modifications by M. De Mattia 13/3/2009 +// --------------------------------------- // - The histograms map was moved to a base class (MuScleFitBase) from which this one inherits. // +// Modifications by M. De Mattia 20/7/2009 +// --------------------------------------- +// - Reworked background fit based on ranges. See comments in the code for more details. // --------------------------------------------------------------------------------------------- #include "MuScleFit.h" @@ -237,6 +240,13 @@ MuScleFit::MuScleFit( const ParameterSet& pset ) : MuScleFitBase( pset ), totalE // This must be set to true if using events generated with Sherpa MuScleFitUtils::sherpa_ = pset.getUntrackedParameter("Sherpa", false); + // Set the cuts on muons to be used in the fit + MuScleFitUtils::minMuonPt_ = pset.getUntrackedParameter("MinMuonPt", 0.); + MuScleFitUtils::maxMuonEta_ = pset.getUntrackedParameter("MaxMuonEta", 6.); + + MuScleFitUtils::debugMassResol_ = pset.getUntrackedParameter("DebugMassResol", false); + // MuScleFitUtils::massResolComponentsStruct MuScleFitUtils::massResolComponents; + // Read the Probs file from database. If false it searches the root file in // MuonAnalysis/MomentumScaleCalibration/test of the active release. // readPdfFromDB = pset.getParameter("readPdfFromDB"); @@ -344,11 +354,6 @@ void MuScleFit::beginOfJob (const EventSetup& eventSetup) { // ----------------- void MuScleFit::endOfJob () { if (debug_>0) cout << "[MuScleFit]: endOfJob" << endl; - - if( loopCounter == 0 ) { - plotter->writeHistoMap(); - delete plotter; - } } // New loop @@ -382,6 +387,12 @@ void MuScleFit::startingNewLoop (unsigned int iLoop) { // ------------------- edm::EDLooper::Status MuScleFit::endOfLoop (const edm::EventSetup& eventSetup, unsigned int iLoop) { + if( loopCounter == 0 ) { + // plotter->writeHistoMap(); + // The destructor will call the writeHistoMap after the cd to the output file + delete plotter; + } + cout << "Ending loop # " << iLoop << endl; // Write the histos to file @@ -693,18 +704,27 @@ edm::EDLooper::Status MuScleFit::duringLoop (const Event & event, const EventSet } if (prob>0) { deltalike = log(prob)*weight; // NB maximum likelihood --> deltalike is maximized - mapHisto_["hLikeVSMu"]->Fill (recMu1, deltalike); - mapHisto_["hLikeVSMu"]->Fill (recMu2, deltalike); - mapHisto_["hLikeVSMuMinus"]->Fill (recMu1, deltalike); - mapHisto_["hLikeVSMuPlus"]->Fill (recMu2, deltalike); + mapHisto_["hLikeVSMu"]->Fill(recMu1, deltalike); + mapHisto_["hLikeVSMu"]->Fill(recMu2, deltalike); + mapHisto_["hLikeVSMuMinus"]->Fill(recMu1, deltalike); + mapHisto_["hLikeVSMuPlus"]->Fill(recMu2, deltalike); double recoMass = (recMu1+recMu2).mass(); if( recoMass != 0 ) { - // IMPORTANT: massResol is a relative resolution - mapHisto_["hResolMassVSMu"]->Fill (recMu1, massResol, -1); - mapHisto_["hResolMassVSMu"]->Fill (recMu2, massResol, +1); - mapHisto_["hFunctionResolMassVSMu"]->Fill (recMu1, massResol/recoMass, -1); - mapHisto_["hFunctionResolMassVSMu"]->Fill (recMu2, massResol/recoMass, +1); + // IMPORTANT: massResol is not a relative resolution + mapHisto_["hResolMassVSMu"]->Fill(recMu1, massResol, -1); + mapHisto_["hResolMassVSMu"]->Fill(recMu2, massResol, +1); + mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu1, massResol/recoMass, -1); + mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu2, massResol/recoMass, +1); + } + + if( MuScleFitUtils::debugMassResol_ ) { + mapHisto_["hdMdPt1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdpt1, -1); + mapHisto_["hdMdPt2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdpt2, +1); + mapHisto_["hdMdPhi1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdphi1, -1); + mapHisto_["hdMdPhi2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdphi2, +1); + mapHisto_["hdMdCotgTh1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdcotgth1, -1); + mapHisto_["hdMdCotgTh2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdcotgth2, +1); } if( !MuScleFitUtils::speedup ) { diff --git a/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitBase.cc b/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitBase.cc index b8579ae7b1689..2e6c95730aa62 100644 --- a/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitBase.cc +++ b/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitBase.cc @@ -42,14 +42,23 @@ void MuScleFitBase::fillHistoMap(TFile* outputFile, unsigned int iLoop) { mapHisto_["hFunctionResolMassVSMu"] = new HResolutionVSPart( outputFile, "hFunctionResolMassVSMu", maxPt, 0, 0.1, 0, 0.1, true ); mapHisto_["hResolPtGenVSMu"] = new HResolutionVSPart( outputFile, "hResolPtGenVSMu", maxPt ); mapHisto_["hResolPtSimVSMu"] = new HResolutionVSPart( outputFile, "hResolPtSimVSMu", maxPt ); - mapHisto_["hResolEtaGenVSMu"] = new HResolutionVSPart( outputFile, "hResolEtaGenVSMu", maxPt ); - mapHisto_["hResolEtaSimVSMu"] = new HResolutionVSPart( outputFile, "hResolEtaSimVSMu", maxPt ); - mapHisto_["hResolThetaGenVSMu"] = new HResolutionVSPart( outputFile, "hResolThetaGenVSMu", maxPt ); - mapHisto_["hResolThetaSimVSMu"] = new HResolutionVSPart( outputFile, "hResolThetaSimVSMu", maxPt ); - mapHisto_["hResolCotgThetaGenVSMu"] = new HResolutionVSPart( outputFile, "hResolCotgThetaGenVSMu", maxPt ); - mapHisto_["hResolCotgThetaSimVSMu"] = new HResolutionVSPart( outputFile, "hResolCotgThetaSimVSMu", maxPt ); - mapHisto_["hResolPhiGenVSMu"] = new HResolutionVSPart( outputFile, "hResolPhiGenVSMu", maxPt ); - mapHisto_["hResolPhiSimVSMu"] = new HResolutionVSPart( outputFile, "hResolPhiSimVSMu", maxPt ); + mapHisto_["hResolEtaGenVSMu"] = new HResolutionVSPart( outputFile, "hResolEtaGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 ); + mapHisto_["hResolEtaSimVSMu"] = new HResolutionVSPart( outputFile, "hResolEtaSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 ); + mapHisto_["hResolThetaGenVSMu"] = new HResolutionVSPart( outputFile, "hResolThetaGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 ); + mapHisto_["hResolThetaSimVSMu"] = new HResolutionVSPart( outputFile, "hResolThetaSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 ); + mapHisto_["hResolCotgThetaGenVSMu"] = new HResolutionVSPart( outputFile, "hResolCotgThetaGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 ); + mapHisto_["hResolCotgThetaSimVSMu"] = new HResolutionVSPart( outputFile, "hResolCotgThetaSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 ); + mapHisto_["hResolPhiGenVSMu"] = new HResolutionVSPart( outputFile, "hResolPhiGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 ); + mapHisto_["hResolPhiSimVSMu"] = new HResolutionVSPart( outputFile, "hResolPhiSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 ); + + if( MuScleFitUtils::debugMassResol_ ) { + mapHisto_["hdMdPt1"] = new HResolutionVSPart( outputFile, "hdMdPt1", maxPt, 0, 100, -3.2, 3.2, true ); + mapHisto_["hdMdPt2"] = new HResolutionVSPart( outputFile, "hdMdPt2", maxPt, 0, 100, -3.2, 3.2, true ); + mapHisto_["hdMdPhi1"] = new HResolutionVSPart( outputFile, "hdMdPhi1", maxPt, 0, 100, -3.2, 3.2, true ); + mapHisto_["hdMdPhi2"] = new HResolutionVSPart( outputFile, "hdMdPhi2", maxPt, 0, 100, -3.2, 3.2, true ); + mapHisto_["hdMdCotgTh1"] = new HResolutionVSPart( outputFile, "hdMdCotgTh1", maxPt, 0, 100, -3.2, 3.2, true ); + mapHisto_["hdMdCotgTh2"] = new HResolutionVSPart( outputFile, "hdMdCotgTh2", maxPt, 0, 100, -3.2, 3.2, true ); + } HTH2D * recoGenHisto = new HTH2D(outputFile, "hPtRecoVsPtGen", "Pt reco vs Pt gen", 120, 0., 120., 120, 0, 120.); (*recoGenHisto)->SetXTitle("Pt gen (GeV)"); @@ -106,11 +115,23 @@ void MuScleFitBase::readProbabilityDistributionsFromFile() TFile * ProbsFile; if ( theMuonType_!=2 ) { //edm::FileInPath file("MuonAnalysis/MomentumScaleCalibration/test/Probs_new_1000_CTEQ.root"); - edm::FileInPath file("MuonAnalysis/MomentumScaleCalibration/test/Probs_new_Horace_CTEQ_1000.root"); - ProbsFile = new TFile (file.fullPath().c_str()); // NNBB need to reset this if MuScleFitUtils::nbins changes + // edm::FileInPath file("MuonAnalysis/MomentumScaleCalibration/test/Probs_new_Horace_CTEQ_1000.root"); + edm::FileInPath file(probabilitiesFileInPath_.c_str()); + // edm::FileInPath file("MuonAnalysis/MomentumScaleCalibration/test/Probs_merge.root"); + + if( probabilitiesFile_ != "" ) { + ProbsFile = new TFile (probabilitiesFile_.c_str()); + cout << "[MuScleFit-Constructor]: Reading TH2D probabilities from " << probabilitiesFile_ << endl; + } + else { + ProbsFile = new TFile (file.fullPath().c_str()); + cout << "[MuScleFit-Constructor]: Reading TH2D probabilities from " << probabilitiesFileInPath_ << endl; + } + // ProbsFile = new TFile ("Probs_new_1000_CTEQ.root"); // NNBB need to reset this if MuScleFitUtils::nbins changes //cout << "[MuScleFit-Constructor]: Reading TH2D probabilities from Probs_new_1000_CTEQ.root file" << endl; - cout << "[MuScleFit-Constructor]: Reading TH2D probabilities from Probs_new_Horace_CTEQ_1000.root file" << endl; + // cout << "[MuScleFit-Constructor]: Reading TH2D probabilities from Probs_new_Horace_CTEQ_1000.root file" << endl; + // cout << "[MuScleFit-Constructor]: Reading TH2D probabilities from Probs_merge.root file" << endl; } else { edm::FileInPath fileSM("MuonAnalysis/MomentumScaleCalibration/test/Probs_SM_1000.root"); ProbsFile = new TFile (fileSM.fullPath().c_str()); // NNBB need to reset this if MuScleFitUtils::nbins changes diff --git a/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitPlotter.cc b/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitPlotter.cc index f2801b40d9064..27049f2670bb0 100644 --- a/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitPlotter.cc +++ b/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitPlotter.cc @@ -1,8 +1,8 @@ // \class MuScleFitPlotter // Plotter for simulated,generated and reco info of muons // -// $Date: 2009/06/04 16:05:03 $ -// $Revision: 1.6 $ +// $Date: 2009/09/08 09:56:33 $ +// $Revision: 1.7 $ // \author C.Mariotti, S.Bolognesi - INFN Torino / T.Dorigo, M.De Mattia - INFN Padova // // ---------------------------------------------------------------------------------- @@ -64,7 +64,9 @@ void MuScleFitPlotter::fillGen1(Handle genParticles) ( pdgId==23 || pdgId==443 || pdgId==100443 || pdgId==553 || pdgId==100553 || pdgId==200553 ) ) { genRes = mcIter->p4(); - mapHisto["hGenRes"]->Fill(genRes); + if( pdgId == 23 ) mapHisto["hGenResZ"]->Fill(genRes); + else if( pdgId == 443 ) mapHisto["hGenResJpsi"]->Fill(genRes); + else if( pdgId == 553 ) mapHisto["hGenResUpsilon1S"]->Fill(genRes); } //Check if it's a muon from a resonance if( status==1 && pdgId==13 ) { @@ -85,9 +87,17 @@ void MuScleFitPlotter::fillGen1(Handle genParticles) cout<<"hgenmumu not found"<Fill(muFromRes.first+muFromRes.second); - mapHisto["hGenResVSMu"]->Fill( muFromRes.first, genRes, 1 ); - mapHisto["hGenResVSMu"]->Fill( muFromRes.second,genRes, -1 ); + mapHisto["hGenMuMuZ"]->Fill(muFromRes.first+muFromRes.second); + mapHisto["hGenMuMuJPsi"]->Fill(muFromRes.first+muFromRes.second); + mapHisto["hGenMuMuUpsilon1S"]->Fill(muFromRes.first+muFromRes.second); + + mapHisto["hGenResVSMuZ"]->Fill( muFromRes.first, genRes, 1 ); + mapHisto["hGenResVSMuZ"]->Fill( muFromRes.second,genRes, -1 ); + mapHisto["hGenResVSMuJPsi"]->Fill( muFromRes.first, genRes, 1 ); + mapHisto["hGenResVSMuJPsi"]->Fill( muFromRes.second,genRes, -1 ); + mapHisto["hGenResVSMuUpsilon1S"]->Fill( muFromRes.first, genRes, 1 ); + mapHisto["hGenResVSMuUpsilon1S"]->Fill( muFromRes.second,genRes, -1 ); + mapHisto["hGenResVsSelf"]->Fill( genRes, genRes, 1 ); } @@ -112,7 +122,12 @@ void MuScleFitPlotter::fillGen2(Handle evtMC) pdgId==553 || pdgId==100553 || pdgId==200553 ) ) { genRes = reco::Particle::LorentzVector((*part)->momentum().px(),(*part)->momentum().py(), (*part)->momentum().pz(),(*part)->momentum().e()); - mapHisto["hGenRes"]->Fill(genRes); + if( pdgId == 23 ) mapHisto["hGenResZ"]->Fill(genRes); + if( pdgId == 443 ) mapHisto["hGenResJPsi"]->Fill(genRes); + if( pdgId == 553 ) { + // cout << "genRes mass = " << CLHEP::HepLorentzVector(genRes.x(),genRes.y(),genRes.z(),genRes.t()).m() << endl; + mapHisto["hGenResUpsilon1S"]->Fill(genRes); + } } //Check if it's a muon from a resonance if (pdgId==13 && status==1) { @@ -140,9 +155,16 @@ void MuScleFitPlotter::fillGen2(Handle evtMC) } } } - mapHisto["hGenMuMu"]->Fill(muFromRes.first+muFromRes.second); - mapHisto["hGenResVSMu"]->Fill( muFromRes.first, genRes, 1 ); - mapHisto["hGenResVSMu"]->Fill( muFromRes.second,genRes, -1 ); + mapHisto["hGenMuMuZ"]->Fill(muFromRes.first+muFromRes.second); + mapHisto["hGenMuMuJPsi"]->Fill(muFromRes.first+muFromRes.second); + mapHisto["hGenMuMuUpsilon1S"]->Fill(muFromRes.first+muFromRes.second); + + mapHisto["hGenResVSMuZ"]->Fill( muFromRes.first, genRes, 1 ); + mapHisto["hGenResVSMuZ"]->Fill( muFromRes.second,genRes, -1 ); + mapHisto["hGenResVSMuJPsi"]->Fill( muFromRes.first, genRes, 1 ); + mapHisto["hGenResVSMuJPsi"]->Fill( muFromRes.second,genRes, -1 ); + mapHisto["hGenResVSMuUpsilon1S"]->Fill( muFromRes.first, genRes, 1 ); + mapHisto["hGenResVSMuUpsilon1S"]->Fill( muFromRes.second,genRes, -1 ); mapHisto["hGenResVsSelf"]->Fill( genRes, genRes, 1 ); } @@ -224,11 +246,19 @@ void MuScleFitPlotter::fillHistoMap() { // Generated Z and muons // --------------------- - mapHisto["hGenRes"] = new HParticle ("hGenRes",3.09685, 3.09695); - mapHisto["hGenMu"] = new HParticle ("hGenMu"); - mapHisto["hGenMuVSEta"] = new HPartVSEta ("hGenMuVSEta"); - mapHisto["hGenMuMu"] = new HParticle ("hGenMuMu",3.09685, 3.09695 ); - mapHisto["hGenResVSMu"] = new HMassVSPart ("hGenResVSMu",3.09685, 3.09695); + mapHisto["hGenResJPsi"] = new HParticle ("hGenResJPsi", 3.09685, 3.09695); + mapHisto["hGenResUpsilon1S"] = new HParticle ("hGenResUpsilon1S", 9., 11.); + mapHisto["hGenResZ"] = new HParticle ("hGenResZ", 60., 120.); + mapHisto["hGenMu"] = new HParticle ("hGenMu"); + mapHisto["hGenMuVSEta"] = new HPartVSEta ("hGenMuVSEta"); + + mapHisto["hGenMuMuJPsi"] = new HParticle ("hGenMuMuJPsi",3.09685, 3.09695 ); + mapHisto["hGenResVSMuJPsi"] = new HMassVSPart ("hGenResVSMuJPsi",3.09685, 3.09695); + mapHisto["hGenMuMuUpsilon1S"] = new HParticle ("hGenMuMuUpsilon1S", 9., 11.); + mapHisto["hGenResVSMuUpsilon1S"] = new HMassVSPart ("hGenResVSMuUpsilon1S", 9., 11.); + mapHisto["hGenMuMuZ"] = new HParticle ("hGenMuMuZ", 60., 120.); + mapHisto["hGenResVSMuZ"] = new HMassVSPart ("hGenResVSMuZ", 60., 120.); + mapHisto["hGenResVsSelf"] = new HMassVSPart ("hGenResVsSelf"); // Simulated resonance and muons diff --git a/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitUtils.cc b/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitUtils.cc index 5c44549656882..98ad42d471880 100644 --- a/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitUtils.cc +++ b/trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitUtils.cc @@ -1,7 +1,7 @@ /** See header file for a class description * - * $Date: 2009/08/14 12:05:22 $ - * $Revision: 1.17 $ + * $Date: 2009/09/08 09:56:33 $ + * $Revision: 1.18 $ * \author S. Bolognesi - INFN Torino / T. Dorigo, M. De Mattia - INFN Padova */ // Some notes: @@ -162,10 +162,10 @@ double MuScleFitUtils::x[][10000]; // Probability matrices and normalization values // --------------------------------------------- -int MuScleFitUtils::nbins = 1000; -double MuScleFitUtils::GLZValue[][1001][1001]; -double MuScleFitUtils::GLZNorm[][1001]; -double MuScleFitUtils::GLValue[][1001][1001]; +int MuScleFitUtils::nbins = 1000; +double MuScleFitUtils::GLZValue[][1001][1001]; +double MuScleFitUtils::GLZNorm[][1001]; +double MuScleFitUtils::GLValue[][1001][1001]; double MuScleFitUtils::GLNorm[][1001]; double MuScleFitUtils::ResMaxSigma[]; @@ -179,7 +179,7 @@ double MuScleFitUtils::massWindowHalfWidth[][3]; int MuScleFitUtils::MuonType; int MuScleFitUtils::MuonTypeForCheckMassWindow; -double MuScleFitUtils::ResGamma[] = {2.4952, 0.0000934, 0.000337, 0.000054, 0.000032, 0.000020}; +double MuScleFitUtils::ResGamma[] = {2.4952, 0.000020, 0.000032, 0.000054, 0.000317, 0.0000932 }; double MuScleFitUtils::ResMass[] = {91.1876, 10.3552, 10.0233, 9.4603, 3.68609, 3.0969}; double MuScleFitUtils::ResMassForBackground[] = { MuScleFitUtils::ResMass[0], (MuScleFitUtils::ResMass[1]+MuScleFitUtils::ResMass[2]+MuScleFitUtils::ResMass[3])/3, @@ -203,6 +203,12 @@ bool MuScleFitUtils::scaleFitNotDone_ = true; bool MuScleFitUtils::sherpa_ = false; +double MuScleFitUtils::minMuonPt_ = 0.; +double MuScleFitUtils::maxMuonEta_ = 6.; + +bool MuScleFitUtils::debugMassResol_; +MuScleFitUtils::massResolComponentsStruct MuScleFitUtils::massResolComponents; + int MuScleFitUtils::iev_ = 0; /////////////////////////////////////////////////////////////////////////////////////////////// @@ -287,40 +293,40 @@ pair MuScleFitUtils::findBestRecoRes( const vector< if (((*Muon1).charge()*(*Muon2).charge())>0) { continue; // This also gets rid of Muon1==Muon2... } - // Accept combinations only if both muons have |eta|<2.4 and pt>3; NOT FOR JPSI!! - // -------------------------------------------------------------- - // if ((*Muon1).p4().Pt()>3.0 && (*Muon2).p4().Pt()>3.0 && - // abs((*Muon1).p4().Eta())<2.4 && abs((*Muon2).p4().Eta())<2.4) { - double mcomb = ((*Muon1).p4()+(*Muon2).p4()).mass(); - double Y = ((*Muon1).p4()+(*Muon2).p4()).Eta(); - if (debug>1) { - cout<<"muon1 "<<(*Muon1).p4().Px()<<", "<<(*Muon1).p4().Py()<<", "<<(*Muon1).p4().Pz()<<", "<<(*Muon1).p4().E()<0 ) { - prob = massProb( mcomb, Y, ires, massResol ); - if( prob>maxprob ) { - if( (*Muon1).charge()<0 ) { // store first the mu minus and then the mu plus - recMuFromBestRes.first = (*Muon1).p4(); - recMuFromBestRes.second = (*Muon2).p4(); - } else { - recMuFromBestRes.first = (*Muon2).p4(); - recMuFromBestRes.second = (*Muon1).p4(); - } - ResFound = true; // NNBB we accept "resonances" even outside mass bounds - maxprob = prob; - } - double deltaMass = fabs(mcomb-ResMass[ires]); - if( deltaMassminMuonPt_ + // ------------------------------------------------------------------------------- + if( (*Muon1).p4().Pt() > minMuonPt_ && (*Muon2).p4().Pt() > minMuonPt_ && + fabs((*Muon1).p4().Eta()) < maxMuonEta_ && fabs((*Muon2).p4().Eta()) < maxMuonEta_ ) { + double mcomb = ((*Muon1).p4()+(*Muon2).p4()).mass(); + double Y = ((*Muon1).p4()+(*Muon2).p4()).Eta(); + if (debug>1) { + cout<<"muon1 "<<(*Muon1).p4().Px()<<", "<<(*Muon1).p4().Py()<<", "<<(*Muon1).p4().Pz()<<", "<<(*Muon1).p4().E()<0 ) { + prob = massProb( mcomb, Y, ires, massResol ); + if( prob>maxprob ) { + if( (*Muon1).charge()<0 ) { // store first the mu minus and then the mu plus + recMuFromBestRes.first = (*Muon1).p4(); + recMuFromBestRes.second = (*Muon2).p4(); + } else { + recMuFromBestRes.first = (*Muon2).p4(); + recMuFromBestRes.second = (*Muon1).p4(); + } + ResFound = true; // NNBB we accept "resonances" even outside mass bounds + maxprob = prob; + } + double deltaMass = fabs(mcomb-ResMass[ires]); + if( deltaMasssigmaPt( pt1,eta1,parval ); @@ -782,8 +797,10 @@ double MuScleFitUtils::massResolution( const lorentzVector& mu1, // Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to // multiply it by pt. double mass_res = sqrt(pow(dmdpt1*sigma_pt1*pt1,2)+pow(dmdpt2*sigma_pt2*pt2,2)+ - pow(dmdphi1*sigma_phi1,2)+pow(dmdphi2*sigma_phi2,2)+ - pow(dmdcotgth1*sigma_cotgth1,2)+pow(dmdcotgth2*sigma_cotgth2,2)); + pow(dmdphi1*sigma_phi1,2)+pow(dmdphi2*sigma_phi2,2)+ + pow(dmdcotgth1*sigma_cotgth1,2)+pow(dmdcotgth2*sigma_cotgth2,2)); + + // double mass_res = sqrt(pow(dmdpt1*sigma_pt1*pt1,2)+pow(dmdpt2*sigma_pt2*pt2,2)); if (debug>19) { cout << " Pt1=" << pt1 << " phi1=" << phi1 << " cotgth1=" << cos(theta1)/sin(theta1) << " - Pt2=" << pt2 @@ -873,6 +890,10 @@ double MuScleFitUtils::probability( const double & mass, const double & massReso const double GLvalue[][1001][1001], const double GLnorm[][1001], const int iRes, const int iY ) { + if( iRes == 0 && iY > 23 ) { + cout << "WARNING: rapidity bin selected = " << iY << " but there are only histograms for the first 24 bins" << endl; + } + double PS = 0.; bool insideProbMassWindow = true; // Interpolate the four values of GLZValue[] in the @@ -1070,6 +1091,11 @@ double MuScleFitUtils::massProb( const double & mass, const double & rapidity, c // NB max value of Z rapidity to be considered is 4. here // ------------------------------------------------------- + // ATTENTION: it should be: + // ------------------------ + // First check the Z, which is divided in 24 rapidity bins + // NB max value of Z rapidity to be considered is 2.4 here + // ------------------------------------------------------- @@ -1083,10 +1109,11 @@ double MuScleFitUtils::massProb( const double & mass, const double & rapidity, c + // ATTENTION: cut on Z rapidity at 2.4 since we only have histograms up to that value pair windowFactors = backgroundHandler->windowFactors( doBackgroundFit[loopCounter], 0 ); if( resfind[0]>0 && checkMassWindow( mass, 0, backgroundHandler->resMass( doBackgroundFit[loopCounter], 0 ), - windowFactors.first, windowFactors.second ) && fabs(rapidity)<4 ) { + windowFactors.first, windowFactors.second ) && fabs(rapidity)<2.4 ) { int iY = (int)(fabs(rapidity)*10.); resConsidered[0] = true; nres += 1;