Skip to content

Commit

Permalink
---
Browse files Browse the repository at this point in the history 
yaml
---
r: 74719
b: "refs/heads/CMSSW_7_1_X"
c: 7437c0c
h: "refs/heads/CMSSW_7_1_X"
i:
  74717: 9b6c880
  74715: 44e759c
  74711: ddfe56e
  74703: 4012db2
  74687: 583ddbf
v: v3
  • Loading branch information
Marco De Mattia committed Oct 5, 2009
1 parent 18cb608 commit f550046
Show file tree
Hide file tree
Showing 8 changed files with 240 additions and 158 deletions.
2 changes: 1 addition & 1 deletion [refs]
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
---
refs/heads/gh-pages: 09c786f70121f131b3715aaf3464996502bbeb7e
"refs/heads/CMSSW_7_1_X": d4aa8aecaf0b3374aa37a33991d54d223e9df430
"refs/heads/CMSSW_7_1_X": 7437c0c5a85104af36f281189eb38fa410d5881b
98 changes: 31 additions & 67 deletions trunk/MuonAnalysis/MomentumScaleCalibration/interface/Functions.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1121,16 +1121,17 @@ class resolutionFunctionType11 : public resolutionFunctionBase<T> {
};

/**
* 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 T>
class resolutionFunctionType12 : public resolutionFunctionBase<T> {
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);
Expand All @@ -1148,32 +1149,51 @@ class resolutionFunctionType12 : public resolutionFunctionBase<T> {
}
// 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<int> & 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 );
}
}
Expand Down Expand Up @@ -1235,62 +1255,6 @@ class resolutionFunctionType12 : public resolutionFunctionBase<T> {
// }
// };

/// 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 T>
class resolutionFunctionType13 : public resolutionFunctionBase<T> {
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<int> & 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 //
// ---------------------------------- //
Expand Down
5 changes: 5 additions & 0 deletions trunk/MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@ class MuScleFitBase
{
public:
MuScleFitBase(const edm::ParameterSet& iConfig) :
probabilitiesFileInPath_( iConfig.getUntrackedParameter<string>( "ProbabilitiesFileInPath" , "MuonAnalysis/MomentumScaleCalibration/test/Probs_new_Horace_CTEQ_1000.root" ) ),
probabilitiesFile_( iConfig.getUntrackedParameter<string>( "ProbabilitiesFile" , "" ) ),
theMuonType_( iConfig.getParameter<int>( "MuonType" ) ),
theMuonLabel_( iConfig.getParameter<edm::InputTag>( "MuonLabel" ) ),
theRootFileName_( iConfig.getUntrackedParameter<string>("OutputFileName") ),
Expand Down Expand Up @@ -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_;
Expand Down
19 changes: 17 additions & 2 deletions trunk/MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitUtils.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
*/

Expand Down Expand Up @@ -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. );

Expand Down
56 changes: 38 additions & 18 deletions trunk/MuonAnalysis/MomentumScaleCalibration/plugins/MuScleFit.cc
View file
Original file line number Diff line number Diff line change
@@ -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:
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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<bool>("Sherpa", false);

// Set the cuts on muons to be used in the fit
MuScleFitUtils::minMuonPt_ = pset.getUntrackedParameter<double>("MinMuonPt", 0.);
MuScleFitUtils::maxMuonEta_ = pset.getUntrackedParameter<double>("MaxMuonEta", 6.);

MuScleFitUtils::debugMassResol_ = pset.getUntrackedParameter<bool>("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<bool>("readPdfFromDB");
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 ) {
Expand Down
43 changes: 32 additions & 11 deletions trunk/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitBase.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)");
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit f550046

Please sign in to comment.