|
|
|
#include <MuScleFitUtils.h>
Classes | |
| struct | byPt |
| struct | massResolComponentsStruct |
Public Member Functions | |
| MuScleFitUtils () | |
| virtual | ~MuScleFitUtils () |
Static Public Member Functions | |
| static lorentzVector | applyBias (const lorentzVector &muon, const int charge) |
| static lorentzVector | applyScale (const lorentzVector &muon, const std::vector< double > &parval, const int charge) |
| static lorentzVector | applyScale (const lorentzVector &muon, double *parval, const int charge) |
| static lorentzVector | applySmearing (const lorentzVector &muon) |
| static bool | checkMassWindow (const double &mass, const double &leftBorder, const double &rightBorder) |
| Method to check if the mass value is within the mass window of the i-th resonance. | |
| static double | computeWeight (const double &mass, const int iev, const bool doUseBkgrWindow=false) |
| static double | deltaPhi (const double &phi1, const double &phi2) |
| static double | deltaPhiNoFabs (const double &phi1, const double &phi2) |
| Without fabs at the end, used to have a symmetric distribution for the resolution fits and variance computations. | |
| static double | deltaR (const double &eta1, const double &eta2, const double &phi1, const double &phi2) |
| static std::pair < lorentzVector, lorentzVector > | findBestRecoRes (const std::vector< reco::LeafCandidate > &muons) |
| static std::pair< SimTrack, SimTrack > | findBestSimuRes (const std::vector< SimTrack > &simMuons) |
| static std::pair < lorentzVector, lorentzVector > | findGenMuFromRes (const edm::HepMCProduct *evtMC) |
| static std::pair < lorentzVector, lorentzVector > | findGenMuFromRes (const reco::GenParticleCollection *genParticles) |
| static std::pair < lorentzVector, lorentzVector > | findSimMuFromRes (const edm::Handle< edm::HepMCProduct > &evtMC, const edm::Handle< edm::SimTrackContainer > &simTracks) |
| static std::vector < TGraphErrors * > | fitMass (TH2F *histo) |
| static std::vector < TGraphErrors * > | fitReso (TH2F *histo) |
| static lorentzVector | fromPtEtaPhiToPxPyPz (const double *ptEtaPhiE) |
| static double | invDimuonMass (const lorentzVector &mu1, const lorentzVector &mu2) |
| static double | massProb (const double &mass, const double &rapidity, const int ires, const double &massResol) |
| static double | massProb (const double &mass, const double &resEta, const double &rapidity, const double &massResol, const std::vector< double > &parval, const bool doUseBkgrWindow=false) |
| static double | massProb (const double &mass, const double &resEta, const double &rapidity, const double &massResol, double *parval, const bool doUseBkgrWindow=false) |
| static double | massResolution (const lorentzVector &mu1, const lorentzVector &mu2, double *parval) |
| static double | massResolution (const lorentzVector &mu1, const lorentzVector &mu2, const std::vector< double > &parval) |
| static double | massResolution (const lorentzVector &mu1, const lorentzVector &mu2) |
| static double | massResolution (const lorentzVector &mu1, const lorentzVector &mu2, std::auto_ptr< double > parval) |
| static void | minimizeLikelihood () |
| static double | probability (const double &mass, const double &massResol, const double GLvalue[][1001][1001], const double GLnorm[][1001], const int iRes, const int iY) |
| Computes the probability given the mass, mass resolution and the arrays with the probabilities and the normalizations. | |
Static Public Attributes | |
| static const int | backgroundFunctionsRegions |
| static BackgroundHandler * | backgroundHandler |
| static TH1D * | backgroundProb_ = 0 |
| static int | BgrFitType = 0 |
| static scaleFunctionBase < std::vector< double > > * | biasFunction = 0 |
| static int | BiasType = 0 |
| static bool | computeMinosErrors_ |
| static int | counter_resprob = 0 |
| static double | crossSection [6] |
| static CrossSectionHandler * | crossSectionHandler |
| static int | debug = 0 |
| static bool | debugMassResol_ |
| static std::vector< int > | doBackgroundFit |
| static std::vector< int > | doCrossSectionFit |
| static std::vector< int > | doResolFit |
| static std::vector< int > | doScaleFit |
| static bool | duringMinos_ = false |
| static int | FitStrategy = 1 |
| static std::vector< std::pair < lorentzVector, lorentzVector > > | genPair |
| static double | GLNorm [6][1001] |
| static double | GLValue [6][1001][1001] |
| static double | GLZNorm [40][1001] |
| static double | GLZValue [40][1001][1001] |
| static int | goodmuon = 0 |
| static int | iev_ = 0 |
| static TH1D * | likelihoodInLoop_ = 0 |
| static unsigned int | loopCounter = 5 |
| static struct MuScleFitUtils::massResolComponentsStruct | massResolComponents |
| static double | massWindowHalfWidth [3][6] |
| static double | maxMuonEtaFirstRange_ = 6. |
| static double | maxMuonEtaSecondRange_ = 100. |
| static double | maxMuonPt_ = 100000000. |
| static bool | minimumShapePlots_ |
| static double | minMuonEtaFirstRange_ = -6. |
| static double | minMuonEtaSecondRange_ = -100. |
| static double | minMuonPt_ = 0. |
| static int | minuitLoop_ = 0 |
| static const double | mMu2 = 0.011163612 |
| static const unsigned int | motherPdgIdArray [6] = {23, 100553, 100553, 553, 100443, 443} |
| static const double | muMass = 0.105658 |
| static int | MuonType |
| static int | MuonTypeForCheckMassWindow |
| static int | nbins = 1000 |
| static unsigned int | normalizationChanged_ = 0 |
| static bool | normalizeLikelihoodByEventNumber_ = true |
| static double | oldNormalization_ = 0. |
| static std::vector< double > | parBgr |
| static std::vector< int > | parBgrFix |
| static std::vector< int > | parBgrOrder |
| static std::vector< double > | parBias |
| static std::vector< double > | parCrossSection |
| static std::vector< int > | parCrossSectionFix |
| static std::vector< int > | parCrossSectionOrder |
| static std::vector< int > | parfix |
| static std::vector< int > | parorder |
| static std::vector< double > | parResol |
| static std::vector< int > | parResolFix |
| static std::vector< int > | parResolOrder |
| static std::vector< double > | parScale |
| static std::vector< int > | parScaleFix |
| static std::vector< int > | parScaleOrder |
| static std::vector< double > | parSmear |
| static std::vector < std::vector< double > > | parvalue |
| static bool | rapidityBinsForZ_ = true |
| static std::vector< std::pair < lorentzVector, lorentzVector > > | ReducedSavedPair |
| static std::vector< int > | resfind |
| static bool | ResFound = false |
| static double | ResGamma [6] = {2.4952, 0.000020, 0.000032, 0.000054, 0.000317, 0.0000932 } |
| static double | ResHalfWidth [6] |
| static double | ResMass [6] = {91.1876, 10.3552, 10.0233, 9.4603, 3.68609, 3.0969} |
| static double | ResMaxSigma [6] |
| static double | ResMinMass [6] = {-99, -99, -99, -99, -99, -99} |
| static int | ResolFitType = 0 |
| static resolutionFunctionBase < double * > * | resolutionFunction = 0 |
| static resolutionFunctionBase < std::vector< double > > * | resolutionFunctionForVec = 0 |
| static TMinuit * | rminPtr_ = 0 |
| static std::vector< std::pair < lorentzVector, lorentzVector > > | SavedPair |
| static bool | scaleFitNotDone_ = true |
| static int | ScaleFitType = 0 |
| static scaleFunctionBase < double * > * | scaleFunction = 0 |
| static scaleFunctionBase < std::vector< double > > * | scaleFunctionForVec = 0 |
| static bool | sherpa_ = false |
| static TH1D * | signalProb_ = 0 |
| static std::vector< std::pair < lorentzVector, lorentzVector > > | simPair |
| static smearFunctionBase * | smearFunction = 0 |
| static int | SmearType = 0 |
| static bool | speedup = false |
| static bool | startWithSimplex_ |
| static const int | totalResNum = 6 |
| static bool | useProbsFile_ = true |
| static double | x [7][10000] |
Definition at line 46 of file MuScleFitUtils.h.
| MuScleFitUtils::MuScleFitUtils | ( | ) | [inline] |
Definition at line 51 of file MuScleFitUtils.h.
{};
| virtual MuScleFitUtils::~MuScleFitUtils | ( | ) | [inline, virtual] |
Definition at line 55 of file MuScleFitUtils.h.
{};
| lorentzVector MuScleFitUtils::applyBias | ( | const lorentzVector & | muon, |
| const int | charge | ||
| ) | [static] |
Definition at line 417 of file MuScleFitUtils.cc.
References biasFunction, gather_cfg::cout, debug, fromPtEtaPhiToPxPyPz(), parBias, and scaleFunctionBase< T >::scale().
{
double ptEtaPhiE[4] = {muon.Pt(),muon.Eta(),muon.Phi(),muon.E()};
if (MuScleFitUtils::debug>1) std::cout << "pt before bias = " << ptEtaPhiE[0] << std::endl;
// Use functors (although not with the () operator)
// Note that we always pass pt, eta and phi, but internally only the needed
// values are used.
// The functors used are takend from the same group used for the scaling
// thus the name of the method used is "scale".
ptEtaPhiE[0] = biasFunction->scale(ptEtaPhiE[0], ptEtaPhiE[1], ptEtaPhiE[2], chg, MuScleFitUtils::parBias);
if (MuScleFitUtils::debug>1) std::cout << "pt after bias = " << ptEtaPhiE[0] << std::endl;
return( fromPtEtaPhiToPxPyPz(ptEtaPhiE) );
}
| lorentzVector MuScleFitUtils::applyScale | ( | const lorentzVector & | muon, |
| const std::vector< double > & | parval, | ||
| const int | charge | ||
| ) | [static] |
Definition at line 437 of file MuScleFitUtils.cc.
References AlCaHLTBitMon_ParallelJobs::p.
Referenced by MuScleFit::duringFastLoop(), and likelihood().
{
double * p = new double[(int)(parval.size())];
// Replaced by auto_ptr, which handles delete at the end
// std::auto_ptr<double> p(new double[(int)(parval.size())]);
// Removed auto_ptr, check massResolution for an explanation.
int id = 0;
for (std::vector<double>::const_iterator it=parval.begin(); it!=parval.end(); ++it, ++id) {
//(&*p)[id] = *it;
// Also ok would be (p.get())[id] = *it;
p[id] = *it;
}
lorentzVector tempScaleVec( applyScale (muon, p, chg) );
delete[] p;
return tempScaleVec;
}
| lorentzVector MuScleFitUtils::applyScale | ( | const lorentzVector & | muon, |
| double * | parval, | ||
| const int | charge | ||
| ) | [static] |
Definition at line 457 of file MuScleFitUtils.cc.
References gather_cfg::cout, debug, fromPtEtaPhiToPxPyPz(), parResol, scaleFunctionBase< T >::scale(), scaleFunction, and edm::shift.
{
double ptEtaPhiE[4] = {muon.Pt(),muon.Eta(),muon.Phi(),muon.E()};
int shift = parResol.size();
if (MuScleFitUtils::debug>1) std::cout << "pt before scale = " << ptEtaPhiE[0] << std::endl;
// the address of parval[shift] is passed as pointer to double. Internally it is used as a normal array, thus:
// array[0] = parval[shift], array[1] = parval[shift+1], ...
ptEtaPhiE[0] = scaleFunction->scale(ptEtaPhiE[0], ptEtaPhiE[1], ptEtaPhiE[2], chg, &(parval[shift]));
if (MuScleFitUtils::debug>1) std::cout << "pt after scale = " << ptEtaPhiE[0] << std::endl;
return( fromPtEtaPhiToPxPyPz(ptEtaPhiE) );
}
| lorentzVector MuScleFitUtils::applySmearing | ( | const lorentzVector & | muon | ) | [static] |
Definition at line 387 of file MuScleFitUtils.cc.
References gather_cfg::cout, debug, eta(), fromPtEtaPhiToPxPyPz(), goodmuon, i, parSmear, phi, smearFunctionBase::smear(), smearFunction, SmearType, x, and detailsBasic3DVector::y.
{
double pt = muon.Pt();
double eta = muon.Eta();
double phi = muon.Phi();
double E = muon.E();
double y[7];
for (int i=0; i<SmearType+3; i++) {
y[i] = x[i][goodmuon%10000];
}
// Use the smear function selected in the constructor
smearFunction->smear( pt, eta, phi, y, parSmear );
if (debug>9) {
std::cout << "Smearing Pt,eta,phi = " << pt << " " << eta << " "
<< phi << "; x = ";
for (int i=0; i<SmearType+3; i++) {
std::cout << y[i];
}
std::cout << std::endl;
}
double ptEtaPhiE[4] = {pt, eta, phi, E};
return( fromPtEtaPhiToPxPyPz(ptEtaPhiE) );
}
| bool MuScleFitUtils::checkMassWindow | ( | const double & | mass, |
| const double & | leftBorder, | ||
| const double & | rightBorder | ||
| ) | [inline, static] |
Method to check if the mass value is within the mass window of the i-th resonance.
Definition at line 1072 of file MuScleFitUtils.cc.
Referenced by computeWeight(), massProb(), and minimizeLikelihood().
{
return( (mass > leftBorder) && (mass < rightBorder) );
}
| double MuScleFitUtils::computeWeight | ( | const double & | mass, |
| const int | iev, | ||
| const bool | doUseBkgrWindow = false |
||
| ) | [static] |
Definition at line 1079 of file MuScleFitUtils.cc.
References backgroundHandler, checkMassWindow(), gather_cfg::cout, debug, doBackgroundFit, ires, loopCounter, resfind, CommonMethods::weight(), and BackgroundHandler::windowBorders().
Referenced by MuScleFit::duringFastLoop(), and likelihood().
{
// Compute weight for this event
// -----------------------------
double weight = 0.;
// Take the highest-mass resonance within bounds
// NB this must be revised once credible estimates of the relative xs of Y(1S), (2S), and (3S)
// are made. Those are priors in the decision of which resonance to assign to an in-between event.
// -----------------------------------------------------------------------------------------------
if( doUseBkgrWindow && (debug > 0) ) std::cout << "using backgrond window for mass = " << mass << std::endl;
// bool useBackgroundWindow = (doBackgroundFit[loopCounter] || doUseBkgrWindow);
bool useBackgroundWindow = (doBackgroundFit[loopCounter]);
for (int ires=0; ires<6; ires++) {
if (resfind[ires]>0 && weight==0.) {
// std::pair<double, double> windowFactor = backgroundHandler->windowFactors( useBackgroundWindow, ires );
std::pair<double, double> windowBorder = backgroundHandler->windowBorders( useBackgroundWindow, ires );
// if( checkMassWindow(mass, ires, backgroundHandler->resMass( useBackgroundWindow, ires ),
// windowFactor.first, windowFactor.second) ) {
if( checkMassWindow(mass, windowBorder.first, windowBorder.second) ) {
weight = 1.0;
if( doUseBkgrWindow && (debug > 0) ) std::cout << "setting weight to = " << weight << std::endl;
}
}
}
return weight;
}
| static double MuScleFitUtils::deltaPhi | ( | const double & | phi1, |
| const double & | phi2 | ||
| ) | [inline, static] |
Definition at line 88 of file MuScleFitUtils.h.
References Pi.
Referenced by MuScleFit::checkDeltaR(), ResolutionAnalyzer::checkDeltaR(), deltaPhiNoFabs(), deltaR(), and HDelta::Fill().
| static double MuScleFitUtils::deltaPhiNoFabs | ( | const double & | phi1, |
| const double & | phi2 | ||
| ) | [inline, static] |
Without fabs at the end, used to have a symmetric distribution for the resolution fits and variance computations.
Definition at line 96 of file MuScleFitUtils.h.
References deltaPhi(), and Pi.
Referenced by ResolutionAnalyzer::analyze(), HCovarianceVSParts::Fill(), and MuScleFit::fillComparisonHistograms().
| static double MuScleFitUtils::deltaR | ( | const double & | eta1, |
| const double & | eta2, | ||
| const double & | phi1, | ||
| const double & | phi2 | ||
| ) | [inline, static] |
Definition at line 103 of file MuScleFitUtils.h.
References deltaPhi(), funct::pow(), and mathSSE::sqrt().
| std::pair< lorentzVector, lorentzVector > MuScleFitUtils::findBestRecoRes | ( | const std::vector< reco::LeafCandidate > & | muons | ) | [static] |
Definition at line 297 of file MuScleFitUtils.cc.
References gather_cfg::cout, debug, ires, massProb(), massResolution(), maxMuonEtaFirstRange_, maxMuonEtaSecondRange_, maxMuonPt_, minMuonEtaFirstRange_, minMuonEtaSecondRange_, minMuonPt_, parResol, resfind, ResFound, ResMass, and useProbsFile_.
Referenced by TestCorrection::analyze(), and MuScleFit::selectMuons().
{
// NB this routine returns the resonance, but it also sets the ResFound flag, which
// is used in MuScleFit to decide whether to use the event or not.
// --------------------------------------------------------------------------------
if (debug>0) std::cout << "In findBestRecoRes" << std::endl;
ResFound = false;
std::pair<lorentzVector, lorentzVector> recMuFromBestRes;
// Choose the best resonance using its mass probability
// ----------------------------------------------------
double maxprob = -0.1;
double minDeltaMass = 999999;
std::pair<reco::LeafCandidate,reco::LeafCandidate> bestMassMuons;
for (std::vector<reco::LeafCandidate>::const_iterator Muon1=muons.begin(); Muon1!=muons.end(); ++Muon1) {
for (std::vector<reco::LeafCandidate>::const_iterator Muon2=Muon1+1; Muon2!=muons.end(); ++Muon2) {
if (((*Muon1).charge()*(*Muon2).charge())>0) {
continue; // This also gets rid of Muon1==Muon2...
}
// To allow the selection of ranges at negative and positive eta independently we define two
// ranges of eta: (minMuonEtaFirstRange_, maxMuonEtaFirstRange_) and (minMuonEtaSecondRange_, maxMuonEtaSecondRange_).
// If the interval selected is simmetric, one only needs to specify the first range. The second has
// default values that accept all muons (minMuonEtaSecondRange_ = -100., maxMuonEtaSecondRange_ = 100.).
double pt1 = (*Muon1).p4().Pt();
double pt2 = (*Muon2).p4().Pt();
double eta1 = (*Muon1).p4().Eta();
double eta2 = (*Muon2).p4().Eta();
if( pt1 > minMuonPt_ && pt1 < maxMuonPt_ &&
pt2 > minMuonPt_ && pt2 < maxMuonPt_ &&
eta1 > minMuonEtaFirstRange_ && eta1 < maxMuonEtaFirstRange_ &&
eta2 > minMuonEtaFirstRange_ && eta2 < maxMuonEtaFirstRange_ &&
eta1 > minMuonEtaSecondRange_ && eta1 < maxMuonEtaSecondRange_ &&
eta2 > minMuonEtaSecondRange_ && eta2 < maxMuonEtaSecondRange_ ) {
double mcomb = ((*Muon1).p4()+(*Muon2).p4()).mass();
double Y = ((*Muon1).p4()+(*Muon2).p4()).Rapidity();
if (debug>1) {
std::cout<<"muon1 "<<(*Muon1).p4().Px()<<", "<<(*Muon1).p4().Py()<<", "<<(*Muon1).p4().Pz()<<", "<<(*Muon1).p4().E()<<std::endl;
std::cout<<"muon2 "<<(*Muon2).p4().Px()<<", "<<(*Muon2).p4().Py()<<", "<<(*Muon2).p4().Pz()<<", "<<(*Muon2).p4().E()<<std::endl;
std::cout<<"mcomb "<<mcomb<<std::endl;}
double massResol = 0.;
if( useProbsFile_ ) {
massResol = massResolution ((*Muon1).p4(), (*Muon2).p4(), parResol);
}
double prob = 0;
for( int ires=0; ires<6; ires++ ) {
if( resfind[ires]>0 ) {
if( useProbsFile_ ) {
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;
}
// if( ResMass[ires] == 0 ) {
// std::cout << "Error: ResMass["<<ires<<"] = " << ResMass[ires] << std::endl;
// exit(1);
// }
double deltaMass = fabs(mcomb-ResMass[ires])/ResMass[ires];
if( deltaMass<minDeltaMass ){
bestMassMuons = std::make_pair((*Muon1),(*Muon2));
minDeltaMass = deltaMass;
}
}
}
}
}
}
//If outside mass window (maxprob==0) then take the two muons with best invariant mass
//(anyway they will not be used in the likelihood calculation, only to fill plots)
if(!maxprob){
if(bestMassMuons.first.charge()<0){
recMuFromBestRes.first = bestMassMuons.first.p4();
recMuFromBestRes.second = bestMassMuons.second.p4();
}
else{
recMuFromBestRes.second = bestMassMuons.first.p4();
recMuFromBestRes.first = bestMassMuons.second.p4();
}
}
return recMuFromBestRes;
}
| std::pair< SimTrack, SimTrack > MuScleFitUtils::findBestSimuRes | ( | const std::vector< SimTrack > & | simMuons | ) | [static] |
Definition at line 260 of file MuScleFitUtils.cc.
References ires, massProb(), and resfind.
Referenced by MuScleFitPlotter::fillSim().
{
std::pair<SimTrack, SimTrack> simMuFromBestRes;
double maxprob = -0.1;
// Double loop on muons
// --------------------
for (std::vector<SimTrack>::const_iterator simMu1=simMuons.begin(); simMu1!=simMuons.end(); simMu1++) {
for (std::vector<SimTrack>::const_iterator simMu2=simMu1+1; simMu2!=simMuons.end(); simMu2++) {
if (((*simMu1).charge()*(*simMu2).charge())>0) {
continue; // this also gets rid of simMu1==simMu2...
}
// Choose the best resonance using its mass. Check Z, Y(3S,2S,1S), Psi(2S), J/Psi in order
// ---------------------------------------------------------------------------------------
double mcomb = ((*simMu1).momentum()+(*simMu2).momentum()).mass();
double Y = ((*simMu1).momentum()+(*simMu2).momentum()).Rapidity();
for (int ires=0; ires<6; ires++) {
if (resfind[ires]>0) {
double prob = massProb( mcomb, Y, ires, 0. );
if (prob>maxprob) {
simMuFromBestRes.first = (*simMu1);
simMuFromBestRes.second = (*simMu2);
maxprob = prob;
}
}
}
}
}
// Return most likely combination of muons making a resonance
// ----------------------------------------------------------
return simMuFromBestRes;
}
| std::pair< lorentzVector, lorentzVector > MuScleFitUtils::findGenMuFromRes | ( | const edm::HepMCProduct * | evtMC | ) | [static] |
Definition at line 2164 of file MuScleFitUtils.cc.
References edm::HepMCProduct::GetEvent(), ires, motherPdgIdArray, resfind, and sherpa_.
{
const HepMC::GenEvent* Evt = evtMC->GetEvent();
std::pair<lorentzVector,lorentzVector> muFromRes;
//Loop on generated particles
for (HepMC::GenEvent::particle_const_iterator part=Evt->particles_begin();
part!=Evt->particles_end(); part++) {
if (fabs((*part)->pdg_id())==13 && (*part)->status()==1) {
bool fromRes = false;
for (HepMC::GenVertex::particle_iterator mother = (*part)->production_vertex()->particles_begin(HepMC::ancestors);
mother != (*part)->production_vertex()->particles_end(HepMC::ancestors); ++mother) {
unsigned int motherPdgId = (*mother)->pdg_id();
// For sherpa the resonance is not saved. The muons from the resonance can be identified
// by having as mother a muon of status 3.
if( sherpa_ ) {
if( motherPdgId == 13 && (*mother)->status() == 3 ) fromRes = true;
}
else {
for( int ires = 0; ires < 6; ++ires ) {
if( motherPdgId == motherPdgIdArray[ires] && resfind[ires] ) fromRes = true;
}
}
}
if(fromRes){
if((*part)->pdg_id()==13)
// muFromRes.first = (*part)->momentum();
muFromRes.first = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
(*part)->momentum().pz(),(*part)->momentum().e()));
else
muFromRes.second = (lorentzVector((*part)->momentum().px(),(*part)->momentum().py(),
(*part)->momentum().pz(),(*part)->momentum().e()));
}
}
}
return muFromRes;
}
| std::pair< lorentzVector, lorentzVector > MuScleFitUtils::findGenMuFromRes | ( | const reco::GenParticleCollection * | genParticles | ) | [static] |
Definition at line 2202 of file MuScleFitUtils.cc.
References gather_cfg::cout, debug, ires, motherPdgIdArray, and resfind.
Referenced by MuScleFitGenFilter::filter().
{
std::pair<lorentzVector,lorentzVector> muFromRes;
//Loop on generated particles
if( debug>0 ) std::cout << "Starting loop on " << genParticles->size() << " genParticles" << std::endl;
for( reco::GenParticleCollection::const_iterator part=genParticles->begin(); part!=genParticles->end(); ++part ) {
if (fabs(part->pdgId())==13 && part->status()==1) {
bool fromRes = false;
unsigned int motherPdgId = part->mother()->pdgId();
if( debug>0 ) {
std::cout << "Found a muon with mother: " << motherPdgId << std::endl;
}
for( int ires = 0; ires < 6; ++ires ) {
if( motherPdgId == motherPdgIdArray[ires] && resfind[ires] ) fromRes = true;
}
if(fromRes){
if(part->pdgId()==13) {
muFromRes.first = part->p4();
if( debug>0 ) std::cout << "Found a genMuon + : " << muFromRes.first << std::endl;
// muFromRes.first = (lorentzVector(part->p4().px(),part->p4().py(),
// part->p4().pz(),part->p4().e()));
}
else {
muFromRes.second = part->p4();
if( debug>0 ) std::cout << "Found a genMuon - : " << muFromRes.second << std::endl;
// muFromRes.second = (lorentzVector(part->p4().px(),part->p4().py(),
// part->p4().pz(),part->p4().e()));
}
}
}
}
return muFromRes;
}
| std::pair< lorentzVector, lorentzVector > MuScleFitUtils::findSimMuFromRes | ( | const edm::Handle< edm::HepMCProduct > & | evtMC, |
| const edm::Handle< edm::SimTrackContainer > & | simTracks | ||
| ) | [static] |
Definition at line 2126 of file MuScleFitUtils.cc.
References configurableAnalysis::GenParticle, ires, motherPdgIdArray, and resfind.
Referenced by MuScleFitPlotter::fillGenSim().
{
//Loop on simulated tracks
std::pair<lorentzVector, lorentzVector> simMuFromRes;
for( edm::SimTrackContainer::const_iterator simTrack=simTracks->begin(); simTrack!=simTracks->end(); ++simTrack ) {
//Chose muons
if (fabs((*simTrack).type())==13) {
//If tracks from IP than find mother
if ((*simTrack).genpartIndex()>0) {
HepMC::GenParticle* gp = evtMC->GetEvent()->barcode_to_particle ((*simTrack).genpartIndex());
if( gp != 0 ) {
for (HepMC::GenVertex::particle_iterator mother = gp->production_vertex()->particles_begin(HepMC::ancestors);
mother!=gp->production_vertex()->particles_end(HepMC::ancestors); ++mother) {
bool fromRes = false;
unsigned int motherPdgId = (*mother)->pdg_id();
for( int ires = 0; ires < 6; ++ires ) {
if( motherPdgId == motherPdgIdArray[ires] && resfind[ires] ) fromRes = true;
}
if( fromRes ) {
if(gp->pdg_id() == 13)
simMuFromRes.first = lorentzVector(simTrack->momentum().px(),simTrack->momentum().py(),
simTrack->momentum().pz(),simTrack->momentum().e());
else
simMuFromRes.second = lorentzVector(simTrack->momentum().px(),simTrack->momentum().py(),
simTrack->momentum().pz(),simTrack->momentum().e());
}
}
}
// else LogDebug("MuScleFitUtils") << "WARNING: no matching genParticle found for simTrack" << std::endl;
}
}
}
return simMuFromRes;
}
| std::vector< TGraphErrors * > MuScleFitUtils::fitMass | ( | TH2F * | histo | ) | [static] |
Definition at line 1845 of file MuScleFitUtils.cc.
References cont, gather_cfg::cout, debug, alignCSCRings::e, i, j, lorentzianPeak(), mergeVDriftHistosByStation::name, python::entryComment::results, and x.
{
if (MuScleFitUtils::debug>0) std::cout << "Fitting " << histo->GetName() << std::endl;
std::vector<TGraphErrors *> results;
// Results of the fit
// ------------------
std::vector<double> Ftop;
std::vector<double> Fwidth;
std::vector<double> Fmass;
std::vector<double> Etop;
std::vector<double> Ewidth;
std::vector<double> Emass;
std::vector<double> Fchi2;
// X bin center and width
// ----------------------
std::vector<double> Xcenter;
std::vector<double> Ex;
// Fit with lorentzian peak
// ------------------------
TF1 *fitFcn = new TF1 ("fitFcn", lorentzianPeak, 70, 110, 3);
fitFcn->SetParameters (100, 3, 91);
fitFcn->SetParNames ("Ftop", "Fwidth", "Fmass");
fitFcn->SetLineWidth (2);
// Fit slices projected along Y from bins in X
// -------------------------------------------
double cont_min = 20; // Minimum number of entries
Int_t binx = histo->GetXaxis()->GetNbins();
// TFile *f= new TFile("prova.root", "recreate");
// histo->Write();
for (int i=1; i<=binx; i++) {
TH1 * histoY = histo->ProjectionY ("", i, i);
// histoY->Write();
double cont = histoY->GetEntries();
if (cont>cont_min) {
histoY->Fit ("fitFcn", "0", "", 70, 110);
double *par = fitFcn->GetParameters();
double *err = fitFcn->GetParErrors();
Ftop.push_back(par[0]);
Fwidth.push_back(par[1]);
Fmass.push_back(par[2]);
Etop.push_back(err[0]);
Ewidth.push_back(err[1]);
Emass.push_back(err[2]);
double chi2 = fitFcn->GetChisquare();
Fchi2.push_back(chi2);
double xx = histo->GetXaxis()->GetBinCenter(i);
Xcenter.push_back(xx);
double ex = 0; // FIXME: you can use the bin width
Ex.push_back(ex);
}
}
// f->Close();
// Put the fit results in arrays for TGraphErrors
// ----------------------------------------------
const int nn = Fmass.size();
double *x = new double[nn];
double *ym = new double[nn];
double *e = new double[nn];
double *eym = new double[nn];
double *yw = new double[nn];
double *eyw = new double[nn];
double *yc = new double[nn];
for (int j=0; j<nn; j++) {
x[j] = Xcenter[j];
ym[j] = Fmass[j];
eym[j] = Emass[j];
yw[j] = Fwidth[j];
eyw[j] = Ewidth[j];
yc[j] = Fchi2[j];
e[j] = Ex[j];
}
// Create TGraphErrors
// -------------------
TString name = histo->GetName();
TGraphErrors *grM = new TGraphErrors (nn, x, ym, e, eym);
grM->SetTitle (name+"_M");
grM->SetName (name+"_M");
TGraphErrors *grW = new TGraphErrors (nn, x, yw, e, eyw);
grW->SetTitle (name+"_W");
grW->SetName (name+"_W");
TGraphErrors *grC = new TGraphErrors (nn, x, yc, e, e);
grC->SetTitle (name+"_chi2");
grC->SetName (name+"_chi2");
// Cleanup
// -------
delete x;
delete ym;
delete eym;
delete yw;
delete eyw;
delete yc;
delete e;
delete fitFcn;
results.push_back(grM);
results.push_back(grW);
results.push_back(grC);
return results;
}
| std::vector< TGraphErrors * > MuScleFitUtils::fitReso | ( | TH2F * | histo | ) | [static] |
Definition at line 1958 of file MuScleFitUtils.cc.
References cont, gather_cfg::cout, alignCSCRings::e, Gaussian(), i, j, mergeVDriftHistosByStation::name, python::entryComment::results, and x.
{
std::cout << "Fitting " << histo->GetName() << std::endl;
std::vector<TGraphErrors *> results;
// Results from fit
// ----------------
std::vector<double> maxs;
std::vector<double> means;
std::vector<double> sigmas;
std::vector<double> chi2s;
std::vector<double> Emaxs;
std::vector<double> Emeans;
std::vector<double> Esigmas;
// X bin center and width
// ----------------------
std::vector<double> Xcenter;
std::vector<double> Ex;
// Fit with a gaussian
// -------------------
TF1 *fitFcn = new TF1 ("fitFunc", Gaussian, -0.2, 0.2, 3);
fitFcn->SetParameters (100, 0, 0.02);
fitFcn->SetParNames ("max", "mean", "sigma");
fitFcn->SetLineWidth (2);
// Fit slices projected along Y from bins in X
// -------------------------------------------
double cont_min = 20; // Minimum number of entries
Int_t binx = histo->GetXaxis()->GetNbins();
for (int i=1; i<=binx; i++) {
TH1 * histoY = histo->ProjectionY ("", i, i);
double cont = histoY->GetEntries();
if (cont>cont_min) {
histoY->Fit ("fitFunc", "0", "", -0.2, 0.2);
double *par = fitFcn->GetParameters();
double *err = fitFcn->GetParErrors();
maxs.push_back (par[0]);
means.push_back (par[1]);
sigmas.push_back (par[2]);
Emaxs.push_back (err[0]);
Emeans.push_back (err[1]);
Esigmas.push_back (err[2]);
double chi2 = fitFcn->GetChisquare();
chi2s.push_back (chi2);
double xx = histo->GetXaxis()->GetBinCenter(i);
Xcenter.push_back (xx);
double ex = 0; // FIXME: you can use the bin width
Ex.push_back (ex);
}
}
// Put the fit results in arrays for TGraphErrors
// ----------------------------------------------
const int nn = means.size();
double *x = new double[nn];
double *ym = new double[nn];
double *e = new double[nn];
double *eym = new double[nn];
double *yw = new double[nn];
double *eyw = new double[nn];
double *yc = new double[nn];
for (int j=0; j<nn; j++) {
x[j] = Xcenter[j];
ym[j] = means[j];
eym[j] = Emeans[j];
// yw[j] = maxs[j];
// eyw[j] = Emaxs[j];
yw[j] = sigmas[j];
eyw[j] = Esigmas[j];
yc[j] = chi2s[j];
e[j] = Ex[j];
}
// Create TGraphErrors
// -------------------
TString name = histo->GetName();
TGraphErrors *grM = new TGraphErrors (nn, x, ym, e, eym);
grM->SetTitle (name+"_mean");
grM->SetName (name+"_mean");
TGraphErrors *grW = new TGraphErrors (nn, x, yw, e, eyw);
grW->SetTitle (name+"_sigma");
grW->SetName (name+"_sigma");
TGraphErrors *grC = new TGraphErrors (nn, x, yc, e, e);
grC->SetTitle (name+"_chi2");
grC->SetName (name+"_chi2");
// Cleanup
// -------
delete x;
delete ym;
delete eym;
delete yw;
delete eyw;
delete yc;
delete e;
delete fitFcn;
results.push_back (grM);
results.push_back (grW);
results.push_back (grC);
return results;
}
| lorentzVector MuScleFitUtils::fromPtEtaPhiToPxPyPz | ( | const double * | ptEtaPhiE | ) | [static] |
Definition at line 476 of file MuScleFitUtils.cc.
References funct::cos(), funct::exp(), muMass, funct::sin(), mathSSE::sqrt(), and tmp.
Referenced by applyBias(), applyScale(), and applySmearing().
{
double px = ptEtaPhiE[0]*cos(ptEtaPhiE[2]);
double py = ptEtaPhiE[0]*sin(ptEtaPhiE[2]);
double tmp = 2*atan(exp(-ptEtaPhiE[1]));
double pz = ptEtaPhiE[0]*cos(tmp)/sin(tmp);
double E = sqrt(px*px+py*py+pz*pz+muMass*muMass);
// lorentzVector corrMu(px,py,pz,E);
// To fix memory leaks, this is to be substituted with
// std::auto_ptr<lorentzVector> corrMu(new lorentzVector(px, py, pz, E));
return lorentzVector(px,py,pz,E);
}
| double MuScleFitUtils::invDimuonMass | ( | const lorentzVector & | mu1, |
| const lorentzVector & | mu2 | ||
| ) | [inline, static] |
Definition at line 493 of file MuScleFitUtils.cc.
Referenced by likelihood(), and minimizeLikelihood().
{
return (mu1+mu2).mass();
}
| double MuScleFitUtils::massProb | ( | const double & | mass, |
| const double & | rapidity, | ||
| const int | ires, | ||
| const double & | massResol | ||
| ) | [static] |
Definition at line 2068 of file MuScleFitUtils.cc.
References gather_cfg::cout, debug, alignCSCRings::e, GL, ires, np, P, Pi, ResGamma, ResMass, w(), and x.
Referenced by MuScleFit::duringFastLoop(), findBestRecoRes(), findBestSimuRes(), likelihood(), and massProb().
{
// This routine computes the likelihood that a given measured mass "measMass" is
// the result of resonance #ires if the resolution expected for the two muons is massResol
// ---------------------------------------------------------------------------------------
double P = 0.;
// Return Lorentz value for zero resolution cases (like MC)
// --------------------------------------------------------
if (massResol==0.) {
if (debug>9) std::cout << "Mass, gamma , mref, width, P: " << mass
<< " " << ResGamma[ires] << " " << ResMass[ires]<< " " << massResol
<< " : used Lorentz P-value" << std::endl;
return (0.5*ResGamma[ires]/TMath::Pi())/((mass-ResMass[ires])*(mass-ResMass[ires])+
.25*ResGamma[ires]*ResGamma[ires]);
}
// NB defined as below, P is not a "probability" but a likelihood that we observe
// a dimuon mass "mass", given mRef, gamma, and massResol. It is what we need for the
// fit which finds the best resolution parameters, though. A definition which is
// more properly a probability is given below (in massProb2()), where the result
// cannot be used to fit resolution parameters because the fit would always prefer
// to set the res parameters to the minimum possible value (best resolution),
// to have a probability close to one of observing any mass.
// -------------------------------------------------------------------------------
// NNBB the following two lines have been replaced with the six following them,
// which provide an improvement of a factor 9 in speed of execution at a
// negligible price in precision.
// ----------------------------------------------------------------------------
// GL->SetParameters(gamma,mRef,mass,massResol);
// P = GL->Integral(mass-5*massResol, mass+5*massResol);
Int_t np = 100;
double * x = new double[np];
double * w = new double[np];
GL->SetParameters (ResGamma[ires], ResMass[ires], mass, massResol);
GL->CalcGaussLegendreSamplingPoints (np, x, w, 0.1e-15);
P = GL->IntegralFast (np, x, w, ResMass[ires]-10*ResGamma[ires], ResMass[ires]+10*ResGamma[ires]);
delete[] x;
delete[] w;
// If we are too far away we set P to epsilon and forget about this event
// ----------------------------------------------------------------------
if (P<1.0e-12) {
P = 1.0e-12;
if (debug>9) std::cout << "Mass, gamma , mref, width, P: " << mass
<< " " << ResGamma[ires] << " " << ResMass[ires] << " " << massResol
<< ": used epsilon" << std::endl;
return P;
}
if (debug>9) std::cout << "Mass, gamma , mref, width, P: " << mass
<< " " << ResGamma[ires] << " " << ResMass[ires] << " " << massResol
<< " " << P << std::endl;
return P;
}
| double MuScleFitUtils::massProb | ( | const double & | mass, |
| const double & | resEta, | ||
| const double & | rapidity, | ||
| const double & | massResol, | ||
| const std::vector< double > & | parval, | ||
| const bool | doUseBkgrWindow = false |
||
| ) | [static] |
Definition at line 686 of file MuScleFitUtils.cc.
References massProb(), and AlCaHLTBitMon_ParallelJobs::p.
{
#ifdef USE_CALLGRIND
CALLGRIND_START_INSTRUMENTATION;
#endif
double * p = new double[(int)(parval.size())];
// Replaced by auto_ptr, which handles delete at the end
// Removed auto_ptr, check massResolution for an explanation.
// std::auto_ptr<double> p(new double[(int)(parval.size())]);
std::vector<double>::const_iterator it = parval.begin();
int id = 0;
for ( ; it!=parval.end(); ++it, ++id) {
// (&*p)[id] = *it;
p[id] = *it;
}
// p must be passed by value as below:
double massProbability = massProb( mass, resEta, rapidity, massResol, p, doUseBkgrWindow );
delete[] p;
#ifdef USE_CALLGRIND
CALLGRIND_STOP_INSTRUMENTATION;
CALLGRIND_DUMP_STATS;
#endif
return massProbability;
}
| double MuScleFitUtils::massProb | ( | const double & | mass, |
| const double & | resEta, | ||
| const double & | rapidity, | ||
| const double & | massResol, | ||
| double * | parval, | ||
| const bool | doUseBkgrWindow = false |
||
| ) | [static] |
Definition at line 851 of file MuScleFitUtils.cc.
References BackgroundHandler::backgroundFunction(), backgroundHandler, backgroundProb_, checkMassWindow(), gather_cfg::cout, crossSectionHandler, debug, doBackgroundFit, GLNorm, GLValue, GLZNorm, GLZValue, i, ires, j, loopCounter, minuitLoop_, MuonType, P, parBgr, CrossSectionHandler::parNum(), parResol, parScale, probability(), rapidityBinsForZ_, CrossSectionHandler::relativeCrossSections(), resfind, ResHalfWidth, BackgroundHandler::resMass(), ResMass, signalProb_, totalResNum, and BackgroundHandler::windowBorders().
{
// This routine computes the likelihood that a given measured mass "measMass" is
// the result of a reference mass ResMass[] if the resolution
// expected for the two muons is massResol.
// This version includes two parameters (the last two in parval, by default)
// to size up the background fraction and its relative normalization with respect
// to the signal shape.
//
// We model the signal probability with a Lorentz L(M,H) of resonance mass M and natural width H
// convoluted with a gaussian G(m,s) of measured mass m and expected mass resolution s,
// by integrating over the intersection of the supports of L and G (which can be made to coincide with
// the region where L is non-zero, given that H<<s in most cases) the product L(M,H)*G(m,s) dx as follows:
//
// GL(m,s) = Int(M-10H,M+10H) [ L(x-M,H) * G(x-m,s) ] dx
//
// The above convolution is computed numerically by an independent root macro, Probs.C, which outputs
// the values in six 1001x1001 grids, one per resonance.
//
// NB THe following block of explanations for background models is outdated, see detailed
// explanations where the code computes PB.
// +++++++++++++++++++++++
// For the background, instead, we have two choices: a linear and an exponential model.
// * For the linear model, we choose a one-parameter form whereby the line is automatically normalized
// in the support [x1,x2] where we defined our "signal region", as follows:
//
// B(x;b) = 1/(x2-x1) + {x - (x2+x1)/2} * b
//
// Defined as above, B(x) is a line passing for the point of coordinates (x1+x2)/2, 1/(x2-x1),
// whose slope b has as support the interval ( -2/[(x1-x2)*(x1+x2)], 2/[(x1-x2)*(x1+x2)] )
// so that B(x) is always positive-definite in [x1,x2].
//
// * For the exponential model, we define B(x;b) as
//
// B(x;b) = b * { exp(-b*x) / [exp(-b*x1)-exp(-b*x2)] }
//
// This one-parameter definition is automatically normalized to unity in [x1,x2], with a parameter
// b which has to be positive in order for the slope to be negative.
// Please note that this model is not useful in most circumstances; a more useful form would be one
// which included a linear component.
// ++++++++++++++++++++++
//
// Once GL(m,s) and B(x;b) are defined, we introduce a further parameter a, such that we can have the
// likelihood control the relative fraction of signal and background. We first normalize GL(m,s) for
// any given s by taking the integral
//
// Int(x1,x2) GL(m,s) dm = K_s
//
// We then define the probability as
//
// P(m,s,a,b) = GL(m,s)/K_s * a + B(x,b) * (1-a)
//
// with a taking on values in the interval [0,1].
// Defined as above, the probability is well-behaved, in the sense that it has a value between 0 and 1,
// and the four parameters m,s,a,b fully control its shape.
//
// It is to be noted that the formulation above requires the computation of two rather time-consuming
// integrals. The one defining GL(m,s) can be stored in a TH2D and loaded by the constructor from a
// file suitably prepared, and this will save loads of computing time.
// ----------------------------------------------------------------------------------------------------
double P = 0.;
int crossSectionParShift = parResol.size() + parScale.size();
// Take the relative cross sections
std::vector<double> relativeCrossSections = crossSectionHandler->relativeCrossSections(&(parval[crossSectionParShift]), resfind);
// for( unsigned int i=0; i<relativeCrossSections.size(); ++i ) {
// std::cout << "relativeCrossSections["<<i<<"] = " << relativeCrossSections[i] << std::endl;
// std::cout << "parval["<<crossSectionParShift+i<<"] = " << parval[crossSectionParShift+i] << std::endl;
// }
// int bgrParShift = crossSectionParShift + parCrossSection.size();
int bgrParShift = crossSectionParShift + crossSectionHandler->parNum();
double Bgrp1 = 0.;
// double Bgrp2 = 0.;
// double Bgrp3 = 0.;
// NB defined as below, P is a non-rigorous "probability" that we observe
// a dimuon mass "mass", given ResMass[], gamma, and massResol. It is what we need for the
// fit which finds the best resolution parameters, though. A definition which is
// more properly a probability is given below (in massProb2()), where the result
// cannot be used to fit resolution parameters because the fit would always prefer
// to set the res parameters to the minimum possible value (best resolution),
// to have a probability close to one of observing any mass.
// -------------------------------------------------------------------------------
// Determine what resonance(s) we have to deal with
// NB for now we assume equal xs for each resonance
// so we do not assign them different weights
// ------------------------------------------------
double PS[6] = {0.};
double PB = 0.;
double PStot[6] = {0.};
// Should be removed because it is not used
bool resConsidered[6] = {false};
bool useBackgroundWindow = (doBackgroundFit[loopCounter] || doUseBkgrWindow);
// bool useBackgroundWindow = (doBackgroundFit[loopCounter]);
// First check the Z, which is divided in 24 rapidity bins
// NB max value of Z rapidity to be considered is 2.4 here
// -------------------------------------------------------
// Do this only if we want to use the rapidity bins for the Z
if( MuScleFitUtils::rapidityBinsForZ_ ) {
// ATTENTION: cut on Z rapidity at 2.4 since we only have histograms up to that value
// std::pair<double, double> windowFactors = backgroundHandler->windowFactors( useBackgroundWindow, 0 );
std::pair<double, double> windowBorders = backgroundHandler->windowBorders( useBackgroundWindow, 0 );
if( resfind[0]>0
// && checkMassWindow( mass, 0,
// backgroundHandler->resMass( useBackgroundWindow, 0 ),
// windowFactors.first, windowFactors.second )
&& checkMassWindow( mass, windowBorders.first, windowBorders.second )
// && fabs(rapidity)<2.4
) {
int iY = (int)(fabs(rapidity)*10.);
if( iY > 23 ) iY = 23;
if (MuScleFitUtils::debug>1) std::cout << "massProb:resFound = 0, rapidity bin =" << iY << std::endl;
// In this case the last value is the rapidity bin
PS[0] = probability(mass, massResol, GLZValue, GLZNorm, 0, iY);
if( PS[0] != PS[0] ) {
std::cout << "ERROR: PS[0] = nan, setting it to 0" << std::endl;
PS[0] = 0;
}
std::pair<double, double> bgrResult = backgroundHandler->backgroundFunction( doBackgroundFit[loopCounter],
&(parval[bgrParShift]), MuScleFitUtils::totalResNum, 0,
resConsidered, ResMass, ResHalfWidth, MuonType, mass, resEta );
Bgrp1 = bgrResult.first;
// When fitting the background we have only one Bgrp1
// When not fitting the background we have many only in a superposition region and this case is treated
// separately after this loop
PB = bgrResult.second;
if( PB != PB ) PB = 0;
PStot[0] = (1-Bgrp1)*PS[0] + Bgrp1*PB;
// PStot[0] *= crossSectionHandler->crossSection(0);
// PStot[0] *= parval[crossSectionParShift];
PStot[0] *= relativeCrossSections[0];
// std::cout << "PStot["<<0<<"] = " << "(1-"<<Bgrp1<<")*"<<PS[0]<<" + "<<Bgrp1<<"*"<<PB<<" = " << PStot[0] << std::endl;
}
else {
if( debug > 0 ) {
std::cout << "Mass = " << mass << " outside range with rapidity = " << rapidity << std::endl;
std::cout << "with resMass = " << backgroundHandler->resMass( useBackgroundWindow, 0 ) << " and left border = " << windowBorders.first << " right border = " << windowBorders.second << std::endl;
}
}
}
// Next check the other resonances
// -------------------------------
int firstRes = 1;
if( !MuScleFitUtils::rapidityBinsForZ_ ) firstRes = 0;
for( int ires=firstRes; ires<6; ++ires ) {
if( resfind[ires] > 0 ) {
// First is left, second is right (returns (1,1) in the case of resonances, it could be improved avoiding the call in this case)
// std::pair<double, double> windowFactor = backgroundHandler->windowFactors( useBackgroundWindow, ires );
std::pair<double, double> windowBorder = backgroundHandler->windowBorders( useBackgroundWindow, ires );
if( checkMassWindow(mass, windowBorder.first, windowBorder.second) ) {
if (MuScleFitUtils::debug>1) std::cout << "massProb:resFound = " << ires << std::endl;
// In this case the rapidity value is instead the resonance index again.
PS[ires] = probability(mass, massResol, GLValue, GLNorm, ires, ires);
std::pair<double, double> bgrResult = backgroundHandler->backgroundFunction( doBackgroundFit[loopCounter],
&(parval[bgrParShift]), MuScleFitUtils::totalResNum, ires,
resConsidered, ResMass, ResHalfWidth, MuonType, mass, resEta );
Bgrp1 = bgrResult.first;
PB = bgrResult.second;
if( PB != PB ) PB = 0;
PStot[ires] = (1-Bgrp1)*PS[ires] + Bgrp1*PB;
if( MuScleFitUtils::debug>0 ) std::cout << "PStot["<<ires<<"] = " << "(1-"<<Bgrp1<<")*"<<PS[ires]<<" + "<<Bgrp1<<"*"<<PB<<" = " << PStot[ires] << std::endl;
PStot[ires] *= relativeCrossSections[ires];
}
}
}
for( int i=0; i<6; ++i ) {
P += PStot[i];
}
if( MuScleFitUtils::signalProb_ != 0 && MuScleFitUtils::backgroundProb_ != 0 ) {
double PStotTemp = 0.;
for( int i=0; i<6; ++i ) {
PStotTemp += PS[i]*relativeCrossSections[i];
}
if( PStotTemp != PStotTemp ) {
std::cout << "ERROR: PStotTemp = nan!!!!!!!!!" << std::endl;
int parnumber = (int)(parResol.size()+parScale.size()+crossSectionHandler->parNum()+parBgr.size());
for( int i=0; i<6; ++i ) {
std::cout << "PS[i] = " << PS[i] << std::endl;
if( PS[i] != PS[i] ) {
std::cout << "mass = " << mass << std::endl;
std::cout << "massResol = " << massResol << std::endl;
for( int j=0; j<parnumber; ++j ) {
std::cout << "parval["<<j<<"] = " << parval[j] << std::endl;
}
}
}
}
if( PStotTemp == PStotTemp ) {
MuScleFitUtils::signalProb_->SetBinContent(MuScleFitUtils::minuitLoop_, MuScleFitUtils::signalProb_->GetBinContent(MuScleFitUtils::minuitLoop_) + PStotTemp);
}
if (debug>0) std::cout << "mass = " << mass << ", P = " << P << ", PStot = " << PStotTemp << ", PB = " << PB << ", bgrp1 = " << Bgrp1 << std::endl;
MuScleFitUtils::backgroundProb_->SetBinContent(MuScleFitUtils::minuitLoop_, MuScleFitUtils::backgroundProb_->GetBinContent(MuScleFitUtils::minuitLoop_) + PB);
}
return P;
}
| double MuScleFitUtils::massResolution | ( | const lorentzVector & | mu1, |
| const lorentzVector & | mu2, | ||
| double * | parval | ||
| ) | [static] |
Definition at line 528 of file MuScleFitUtils.cc.
References funct::cos(), counter_resprob, gather_cfg::cout, resolutionFunctionBase< T >::covPt1Pt2(), debug, debugMassResol_, MuScleFitUtils::massResolComponentsStruct::dmdcotgth1, MuScleFitUtils::massResolComponentsStruct::dmdcotgth2, MuScleFitUtils::massResolComponentsStruct::dmdphi1, MuScleFitUtils::massResolComponentsStruct::dmdphi2, MuScleFitUtils::massResolComponentsStruct::dmdpt1, MuScleFitUtils::massResolComponentsStruct::dmdpt2, funct::exp(), ires, LogDebug, massResolComponents, mMu2, funct::pow(), resfind, ResHalfWidth, ResMass, ResMaxSigma, resolutionFunction, resolutionFunctionBase< T >::sigmaCotgTh(), resolutionFunctionBase< T >::sigmaPhi(), resolutionFunctionBase< T >::sigmaPt(), funct::sin(), and mathSSE::sqrt().
{
// We use the following formula:
//
// M = sqrt ( (E1+E2)^2 - (P1+P2)^2 )
//
// where we express E and P as a function of Pt, phi, and theta:
//
// E = sqrt ( Pt^2*(1+cotg(theta)^2) + M_mu^2 )
// Px = Pt*cos(phi), Py = Pt*sin(phi), Pz = Pt*cotg(theta)
//
// from which we find
//
// M = sqrt( 2*M_mu^2 + 2*sqrt(Pt1^2/sin(theta1)^2 + M_mu^2)*sqrt(Pt2^2/sin(theta2)^2 + M_mu^2) -
// 2*Pt1*Pt2* ( cos(phi1-phi2) + cotg(theta1)*cotg(theta2) ) )
//
// and derive WRT Pt1, Pt2, phi1, phi2, theta1, theta2 to get the resolution.
// --------------------------------------------------------------------------
double mass = (mu1+mu2).mass();
double pt1 = mu1.Pt();
double phi1 = mu1.Phi();
double eta1 = mu1.Eta();
double theta1 = 2*atan(exp(-eta1));
// double cotgTheta1 = cos(theta1)/sin(theta1);
double pt2 = mu2.Pt();
double phi2 = mu2.Phi();
double eta2 = mu2.Eta();
double theta2 = 2*atan(exp(-eta2));
// double cotgTheta2 = cos(theta2)/sin(theta2);
// double mass_check = sqrt(2*mMu2+2*sqrt(std::pow(pt1/sin(theta1),2)+mMu2)*sqrt(std::pow(pt2/sin(theta2),2)+mMu2)-
// 2*pt1*pt2*(cos(phi1-phi2)+1/(tan(theta1)*tan(theta2))));
// ATTENTION: need to compute 1/tan(theta) as cos(theta)/sin(theta) because the latter diverges for theta=pi/2
// -----------------------------------------------------------------------------------------------------------
double dmdpt1 = (pt1/std::pow(sin(theta1),2)*sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2))-
pt2*(cos(phi1-phi2)+cos(theta1)*cos(theta2)/(sin(theta1)*sin(theta2))))/mass;
double dmdpt2 = (pt2/std::pow(sin(theta2),2)*sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2))-
pt1*(cos(phi2-phi1)+cos(theta2)*cos(theta1)/(sin(theta2)*sin(theta1))))/mass;
double dmdphi1 = pt1*pt2/mass*sin(phi1-phi2);
double dmdphi2 = pt2*pt1/mass*sin(phi2-phi1);
// double dmdtheta1 = (-std::pow(pt1/sin(theta1),2)/tan(theta1)*
// sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2))+
// 2*pt1*pt2/(tan(theta2)*std::pow(sin(theta1),2)))/mass;
// double dmdtheta2 = (-std::pow(pt2/sin(theta2),2)/tan(theta2)*
// sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2))+
// 2*pt2*pt1/(tan(theta1)*std::pow(sin(theta2),2)))/mass;
double dmdcotgth1 = (pt1*pt1*cos(theta1)/sin(theta1)*
sqrt((std::pow(pt2/sin(theta2),2)+mMu2)/(std::pow(pt1/sin(theta1),2)+mMu2)) -
pt1*pt2*cos(theta2)/sin(theta2))/mass;
double dmdcotgth2 = (pt2*pt2*cos(theta2)/sin(theta2)*
sqrt((std::pow(pt1/sin(theta1),2)+mMu2)/(std::pow(pt2/sin(theta2),2)+mMu2)) -
pt2*pt1*cos(theta1)/sin(theta1))/mass;
if( debugMassResol_ ) {
massResolComponents.dmdpt1 = dmdpt1;
massResolComponents.dmdpt2 = dmdpt2;
massResolComponents.dmdphi1 = dmdphi1;
massResolComponents.dmdphi2 = dmdphi2;
massResolComponents.dmdcotgth1 = dmdcotgth1;
massResolComponents.dmdcotgth2 = dmdcotgth2;
}
// Resolution parameters:
// ----------------------
double sigma_pt1 = resolutionFunction->sigmaPt( pt1,eta1,parval );
double sigma_pt2 = resolutionFunction->sigmaPt( pt2,eta2,parval );
double sigma_phi1 = resolutionFunction->sigmaPhi( pt1,eta1,parval );
double sigma_phi2 = resolutionFunction->sigmaPhi( pt2,eta2,parval );
double sigma_cotgth1 = resolutionFunction->sigmaCotgTh( pt1,eta1,parval );
double sigma_cotgth2 = resolutionFunction->sigmaCotgTh( pt2,eta2,parval );
double cov_pt1pt2 = resolutionFunction->covPt1Pt2( pt1, eta1, pt2, eta2, parval );
// Sigma_Pt is defined as a relative sigmaPt/Pt for this reason we need to
// multiply it by pt.
double mass_res = sqrt(std::pow(dmdpt1*sigma_pt1*pt1,2)+std::pow(dmdpt2*sigma_pt2*pt2,2)+
std::pow(dmdphi1*sigma_phi1,2)+std::pow(dmdphi2*sigma_phi2,2)+
std::pow(dmdcotgth1*sigma_cotgth1,2)+std::pow(dmdcotgth2*sigma_cotgth2,2)+
2*dmdpt1*dmdpt2*cov_pt1pt2);
// if( sigma_cotgth1 < 0 || sigma_cotgth2 < 0 ) {
// std::cout << "WARNING: sigma_cotgth1 = " << sigma_cotgth1 << std::endl;
// std::cout << "WARNING: sigma_cotgth2 = " << sigma_cotgth2 << std::endl;
// std::cout << "mass_res = " << mass_res << std::endl;
// }
// double mass_res = sqrt(std::pow(dmdpt1*sigma_pt1*pt1,2)+std::pow(dmdpt2*sigma_pt2*pt2,2));
if (debug>19) {
std::cout << " Pt1=" << pt1 << " phi1=" << phi1 << " cotgth1=" << cos(theta1)/sin(theta1) << " - Pt2=" << pt2
<< " phi2=" << phi2 << " cotgth2=" << cos(theta2)/sin(theta2) << std::endl;
std::cout << " P[0]="
<< parval[0] << " P[1]=" << parval[1] << "P[2]=" << parval[2] << " P[3]=" << parval[3] << std::endl;
std::cout << " Dmdpt1= " << dmdpt1 << " dmdpt2= " << dmdpt2 << " sigma_pt1="
<< sigma_pt1 << " sigma_pt2=" << sigma_pt2 << std::endl;
std::cout << " Dmdphi1= " << dmdphi1 << " dmdphi2= " << dmdphi2 << " sigma_phi1="
<< sigma_phi1 << " sigma_phi2=" << sigma_phi2 << std::endl;
std::cout << " Dmdcotgth1= " << dmdcotgth1 << " dmdcotgth2= " << dmdcotgth2
<< " sigma_cotgth1="
<< sigma_cotgth1 << " sigma_cotgth2=" << sigma_cotgth2 << std::endl;
std::cout << " Mass resolution (pval) for muons of Pt = " << pt1 << " " << pt2
<< " : " << mass << " +- " << mass_res << std::endl;
}
// Debug std::cout
// ----------
bool didit = false;
for (int ires=0; ires<6; ires++) {
if (!didit && resfind[ires]>0 && fabs(mass-ResMass[ires])<ResHalfWidth[ires]) {
if (mass_res>ResMaxSigma[ires] && counter_resprob<100) {
counter_resprob++;
LogDebug("MuScleFitUtils") << "RESOLUTION PROBLEM: ires=" << ires << std::endl;
// std::cout << "RESOLUTION PROBLEM: ires=" << ires << std::endl;
// std::cout << "---------------------------" << std::endl;
// std::cout << " Pt1=" << pt1 << " phi1=" << phi1 << " cotgth1=" << cos(theta1)/sin(theta1) << " - Pt2=" << pt2
// << " phi2=" << phi2 << " cotgth2=" << cos(theta2)/sin(theta2) << std::endl;
// if (ResolFitType==1)
// std::cout << " P[0]="
// << parval[0] << " P[1]=" << parval[1] << "P[2]=" << parval[2] << " P[3]=" << parval[3]
// << " P[4]=" << parval[4] << " P[5]=" << parval[5] << std::endl;
// if (ResolFitType==2)
// std::cout << " P[0]="
// << parval[0] << " P[1]=" << parval[1] << "P[2]=" << parval[2] << " P[3]=" << parval[3] << std::endl;
// std::cout << " Dmdpt1= " << dmdpt1 << " dmdpt2= " << dmdpt2 << " sigma_pt1="
// << sigma_pt1 << " sigma_pt2=" << sigma_pt2 << std::endl;
// std::cout << " Dmdphi1= " << dmdphi1 << " dmdphi2= " << dmdphi2 << " sigma_phi1="
// << sigma_phi1 << " sigma_phi2=" << sigma_phi2 << std::endl;
// std::cout << " Dmdcotgth1= " << dmdcotgth1 << " dmdcotgth2= " << dmdcotgth2
// << " sigma_cotgth1="
// << sigma_cotgth1 << " sigma_cotgth2=" << sigma_cotgth2 << std::endl;
// std::cout << " Mass resolution (pval) for muons of Pt = " << pt1 << " " << pt2
// << " : " << mass << " +- " << mass_res << std::endl;
// std::cout << "---------------------------" << std::endl;
didit = true;
}
}
}
// if( mass_res != mass_res ) {
// std::cout << "MASS_RESOL_NAN PROBLEM:" << std::endl;
// std::cout << "---------------------------" << std::endl;
// std::cout << " Pt1=" << pt1 << " phi1=" << phi1 << " cotgth1=" << cos(theta1)/sin(theta1) << " - Pt2=" << pt2
// << " phi2=" << phi2 << " cotgth2=" << cos(theta2)/sin(theta2) << std::endl;
// std::cout << " P[0]=" << parval[0] << " P[1]=" << parval[1] << "P[2]=" << parval[2] << " P[3]=" << parval[3] << " P[4]=" << parval[4] << std::endl;
// std::cout << " Dmdpt1= " << dmdpt1 << " dmdpt2= " << dmdpt2 << " sigma_pt1=" << sigma_pt1 << " sigma_pt2=" << sigma_pt2 << std::endl;
// std::cout << " Dmdphi1= " << dmdphi1 << " dmdphi2= " << dmdphi2 << " sigma_phi1=" << sigma_phi1 << " sigma_phi2=" << sigma_phi2 << std::endl;
// std::cout << " Dmdcotgth1= " << dmdcotgth1 << " dmdcotgth2= " << dmdcotgth2 << " sigma_cotgth1=" << sigma_cotgth1 << " sigma_cotgth2=" << sigma_cotgth2 << std::endl;
// std::cout << " Mass resolution (pval) for muons of Pt = " << pt1 << " " << pt2 << " : " << mass << " +- " << mass_res << std::endl;
// std::cout << "---------------------------" << std::endl;
// }
return mass_res;
}
| double MuScleFitUtils::massResolution | ( | const lorentzVector & | mu1, |
| const lorentzVector & | mu2, | ||
| const std::vector< double > & | parval | ||
| ) | [static] |
Definition at line 501 of file MuScleFitUtils.cc.
References massResolution(), and AlCaHLTBitMon_ParallelJobs::p.
{
// double * p = new double[(int)(parval.size())];
// Replaced by auto_ptr, which handles delete at the end
// --------- //
// ATTENTION //
// --------- //
// auto_ptr calls delete, not delete[] and thus it must
// not be used with arrays. There are alternatives see
// e.g.: https://www.gotw.ca/gotw/042.htm. The best
// alternative seems to be to switch to vector though.
// std::auto_ptr<double> p(new double[(int)(parval.size())]);
double * p = new double[(int)(parval.size())];
std::vector<double>::const_iterator it = parval.begin();
int id = 0;
for ( ; it!=parval.end(); ++it, ++id) {
// (&*p)[id] = *it;
p[id] = *it;
}
double massRes = massResolution (mu1, mu2, p);
delete[] p;
return massRes;
}
| static double MuScleFitUtils::massResolution | ( | const lorentzVector & | mu1, |
| const lorentzVector & | mu2 | ||
| ) | [static] |
| static double MuScleFitUtils::massResolution | ( | const lorentzVector & | mu1, |
| const lorentzVector & | mu2, | ||
| std::auto_ptr< double > | parval | ||
| ) | [static] |
| void MuScleFitUtils::minimizeLikelihood | ( | ) | [static] |
Definition at line 1112 of file MuScleFitUtils.cc.
References backgroundHandler, backgroundProb_, BgrFitType, svgfig::canvas(), CastorDataFrameFilter_impl::check(), checkMassWindow(), computeMinosErrors_, gather_cfg::cout, crossSectionHandler, debug, doBackgroundFit, doCrossSectionFit, doResolFit, doScaleFit, duringMinos_, FitStrategy, i, invDimuonMass(), iorder, ires, likelihood(), likelihoodInLoop_, loopCounter, massWindowHalfWidth, minimumShapePlots_, minuitLoop_, MuonType, mergeVDriftHistosByStation::name, normalizationChanged_, parBgr, parBgrFix, parBgrOrder, parCrossSection, parCrossSectionOrder, parfix, scaleFunctionBase< T >::parNum(), CrossSectionHandler::parNum(), resolutionFunctionBase< T >::parNum(), parorder, parResol, parResolFix, parResolOrder, parScale, parScaleFix, parScaleOrder, parvalue, ReducedSavedPair, BackgroundHandler::regionsParNum(), CrossSectionHandler::relativeCrossSections(), CrossSectionHandler::releaseParameters(), BackgroundHandler::rescale(), scaleFunctionBase< T >::resetParameters(), resfind, BackgroundHandler::resMass(), ResMass, ResolFitType, rminPtr_, SavedPair, scaleFitNotDone_, ScaleFitType, scaleFunction, CrossSectionHandler::setParameters(), BackgroundHandler::setParameters(), resolutionFunctionBase< T >::setParameters(), scaleFunctionBase< T >::setParameters(), signalProb_, startWithSimplex_, cmsPerfCommons::Step, BackgroundHandler::unlockParameter(), and BackgroundHandler::windowBorders().
Referenced by MuScleFit::endOfFastLoop().
{
// Output file with fit parameters resulting from minimization
// -----------------------------------------------------------
ofstream FitParametersFile;
FitParametersFile.open ("FitParameters.txt", std::ios::app);
FitParametersFile << "Fitting with resolution, scale, bgr function # "
<< ResolFitType << " " << ScaleFitType << " " << BgrFitType
<< " - Iteration " << loopCounter << std::endl;
// Fill parvalue and other vectors needed for the fitting
// ------------------------------------------------------
// ----- //
// FIXME //
// ----- //
// this was changed to verify the possibility that fixed parameters influence the errors.
// It must be 0 otherwise the parameters for resonances will not be passed by minuit (will be always 0).
// Should be removed.
int parForResonanceWindows = 0;
// int parnumber = (int)(parResol.size()+parScale.size()+parCrossSection.size()+parBgr.size() - parForResonanceWindows);
int parnumber = (int)(parResol.size()+parScale.size()+crossSectionHandler->parNum()+parBgr.size() - parForResonanceWindows);
int parnumberAll = (int)(parResol.size()+parScale.size()+crossSectionHandler->parNum()+parBgr.size());
// parvalue is a std::vector<std::vector<double> > storing all the parameters from all the loops
parvalue.push_back(parResol);
std::vector<double> *tmpVec = &(parvalue.back());
// If this is not the first loop we want to start from neutral values
// Otherwise the scale will start with values correcting again a bias
// that is already corrected.
if( scaleFitNotDone_ ) {
tmpVec->insert( tmpVec->end(), parScale.begin(), parScale.end() );
std::cout << "scaleFitNotDone: tmpVec->size() = " << tmpVec->size() << std::endl;
}
else {
scaleFunction->resetParameters(tmpVec);
std::cout << "scaleFitDone: tmpVec->size() = " << tmpVec->size() << std::endl;
}
tmpVec->insert( tmpVec->end(), parCrossSection.begin(), parCrossSection.end() );
tmpVec->insert( tmpVec->end(), parBgr.begin(), parBgr.end() );
int i = 0;
std::vector<double>::const_iterator it = tmpVec->begin();
for( ; it != tmpVec->end(); ++it, ++i ) {
std::cout << "tmpVec["<<i<<"] = " << *it << std::endl;
}
// Empty vector of size = number of cross section fitted parameters. Note that the cross section
// fit works in a different way than the others and it uses ratios of the paramters passed via cfg.
// We use this empty vector for compatibility with the rest of the structure.
std::vector<int> crossSectionParNumSizeVec( MuScleFitUtils::crossSectionHandler->parNum(), 0 );
std::vector<int> parfix(parResolFix);
parfix.insert( parfix.end(), parScaleFix.begin(), parScaleFix.end() );
parfix.insert( parfix.end(), crossSectionParNumSizeVec.begin(), crossSectionParNumSizeVec.end() );
parfix.insert( parfix.end(), parBgrFix.begin(), parBgrFix.end() );
std::vector<int> parorder(parResolOrder);
parorder.insert( parorder.end(), parScaleOrder.begin(), parScaleOrder.end() );
parorder.insert( parorder.end(), crossSectionParNumSizeVec.begin(), crossSectionParNumSizeVec.end() );
parorder.insert( parorder.end(), parBgrOrder.begin(), parBgrOrder.end() );
// This is filled later
std::vector<double> parerr(3*parnumberAll,0.);
if (debug>19) {
std::cout << "[MuScleFitUtils-minimizeLikelihood]: Parameters before likelihood " << std::endl;
for (unsigned int i=0; i<(unsigned int)parnumberAll; i++) {
std::cout << " Par # " << i << " = " << parvalue[loopCounter][i] << " : free = " << parfix[i] << "; order = "
<< parorder[i] << std::endl;
}
}
// // Background rescaling from regions to resonances
// // -----------------------------------------------
// // If we are in a loop > 0 and we are not fitting the background, but we have fitted it in the previous iteration
// if( loopCounter > 0 && !(doBackgroundFit[loopCounter]) && doBackgroundFit[loopCounter-1] ) {
// // This rescales from regions to resonances
// int localMuonType = MuonType;
// if( MuonType > 2 ) localMuonType = 2;
// backgroundHandler->rescale( parBgr, ResMass, massWindowHalfWidth[localMuonType],
// MuScleFitUtils::SavedPair);
// }
// Init Minuit
// -----------
TMinuit rmin (parnumber);
rminPtr_ = &rmin;
rmin.SetFCN (likelihood); // Unbinned likelihood
// Standard initialization of minuit parameters:
// sets input to be $stdin, output to be $stdout
// and saving to a file.
rmin.mninit (5, 6, 7);
int ierror = 0;
int istat;
double arglis[4];
arglis[0] = FitStrategy; // Strategy 1 or 2
// 1 standard
// 2 try to improve minimum (slower)
rmin.mnexcm ("SET STR", arglis, 1, ierror);
arglis[0] = 10001;
// Set the random seed for the generator used in SEEk to a fixed value for reproducibility
rmin.mnexcm("SET RAN", arglis, 1, ierror);
// Set fit parameters
// ------------------
double * Start = new double[parnumberAll];
double * Step = new double[parnumberAll];
double * Mini = new double[parnumberAll];
double * Maxi = new double[parnumberAll];
int * ind = new int[parnumberAll]; // Order of release of parameters
TString * parname = new TString[parnumberAll];
MuScleFitUtils::resolutionFunctionForVec->setParameters( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, MuonType );
// Take the number of parameters in the resolutionFunction and displace the arrays passed to the scaleFunction
int resParNum = MuScleFitUtils::resolutionFunctionForVec->parNum();
MuScleFitUtils::scaleFunctionForVec->setParameters( &(Start[resParNum]), &(Step[resParNum]), &(Mini[resParNum]), &(Maxi[resParNum]),
&(ind[resParNum]), &(parname[resParNum]), parScale, parScaleOrder, MuonType );
// Initialize cross section parameters
int crossSectionParShift = resParNum + MuScleFitUtils::scaleFunctionForVec->parNum();
MuScleFitUtils::crossSectionHandler->setParameters( &(Start[crossSectionParShift]), &(Step[crossSectionParShift]), &(Mini[crossSectionParShift]),
&(Maxi[crossSectionParShift]), &(ind[crossSectionParShift]), &(parname[crossSectionParShift]),
parCrossSection, parCrossSectionOrder, resfind );
// Initialize background parameters
int bgrParShift = crossSectionParShift + crossSectionHandler->parNum();
MuScleFitUtils::backgroundHandler->setParameters( &(Start[bgrParShift]), &(Step[bgrParShift]), &(Mini[bgrParShift]), &(Maxi[bgrParShift]),
&(ind[bgrParShift]), &(parname[bgrParShift]), parBgr, parBgrOrder, MuonType );
for( int ipar=0; ipar<parnumber; ++ipar ) {
std::cout << "parname["<<ipar<<"] = " << parname[ipar] << std::endl;
std::cout << "Start["<<ipar<<"] = " << Start[ipar] << std::endl;
std::cout << "Step["<<ipar<<"] = " << Step[ipar] << std::endl;
std::cout << "Mini["<<ipar<<"] = " << Mini[ipar] << std::endl;
std::cout << "Maxi["<<ipar<<"] = " << Maxi[ipar] << std::endl;
rmin.mnparm( ipar, parname[ipar], Start[ipar], Step[ipar], Mini[ipar], Maxi[ipar], ierror );
// Testing without limits
// rmin.mnparm( ipar, parname[ipar], Start[ipar], Step[ipar], 0, 0, ierror );
}
// Do minimization
// ---------------
if (debug>19)
std::cout << "[MuScleFitUtils-minimizeLikelihood]: Starting minimization" << std::endl;
double fmin;
double fdem;
double errdef;
int npari;
int nparx;
rmin.mnexcm ("call fcn", arglis, 1, ierror);
// First, fix all parameters
// -------------------------
if (debug>19)
std::cout << "[MuScleFitUtils-minimizeLikelihood]: First fix all parameters ...";
for (int ipar=0; ipar<parnumber; ipar++) {
rmin.FixParameter (ipar);
}
// Then release them in the specified order and refit
// --------------------------------------------------
if (debug>19) std::cout << " Then release them in order..." << std::endl;
TString name;
double pval;
double pmin;
double pmax;
double errp;
double errl;
double errh;
int ivar;
double erro;
double cglo;
int n_times = 0;
// n_times = number of loops required to unlock all parameters.
if (debug>19) std::cout << "Before scale parNum" << std::endl;
int scaleParNum = scaleFunction->parNum();
if (debug>19) std::cout << "After scale parNum" << std::endl;
// edm::LogInfo("minimizeLikelihood") << "number of parameters for scaleFunction = " << scaleParNum << std::endl;
// edm::LogInfo("minimizeLikelihood") << "number of parameters for resolutionFunction = " << resParNum << std::endl;
// edm::LogInfo("minimizeLikelihood") << "number of parameters for cross section = " << crossSectionHandler->parNum() << std::endl;
// edm::LogInfo("minimizeLikelihood") << "number of parameters for backgroundFunction = " << parBgr.size() << std::endl;
std::cout << "number of parameters for scaleFunction = " << scaleParNum << std::endl;
std::cout << "number of parameters for resolutionFunction = " << resParNum << std::endl;
std::cout << "number of parameters for cross section = " << crossSectionHandler->parNum() << std::endl;
std::cout << "number of parameters for backgroundFunction = " << parBgr.size() << std::endl;
// edm::LogInfo("minimizeLikelihood") << "number of parameters for backgroundFunction = " << backgroundFunction->parNum() << std::endl;
for (int i=0; i<parnumber; i++) {
// NB ind[] has been set as parorder[] previously
if (n_times<ind[i]) {
edm::LogInfo("minimizeLikelihood") << "n_times = " << n_times << ", ind["<<i<<"] = " << ind[i] << ", scaleParNum = " << scaleParNum << ", doScaleFit["<<loopCounter<<"] = " << doScaleFit[loopCounter] << std::endl;
// Set the n_times only if we will do the fit
if ( i<resParNum ) {
if( doResolFit[loopCounter] ) n_times = ind[i];
}
else if( i<resParNum+scaleParNum ) {
if( doScaleFit[loopCounter] ) n_times = ind[i];
}
else if( doBackgroundFit[loopCounter] ) n_times = ind[i];
}
}
for (int iorder=0; iorder<n_times+1; iorder++) { // Repeat fit n_times times
std::cout << "Starting minimization " << iorder << " of " << n_times << std::endl;
bool somethingtodo = false;
// Use parameters from cfg to select which fit to do
// -------------------------------------------------
if( doResolFit[loopCounter] ) {
// Release resolution parameters and fit them
// ------------------------------------------
for( unsigned int ipar=0; ipar<parResol.size(); ++ipar ) {
if( parfix[ipar]==0 && ind[ipar]==iorder ) {
rmin.Release( ipar );
somethingtodo = true;
}
}
}
if( doScaleFit[loopCounter] ) {
// Release scale parameters and fit them
// -------------------------------------
for( unsigned int ipar=parResol.size(); ipar<parResol.size()+parScale.size(); ++ipar ) {
if( parfix[ipar]==0 && ind[ipar]==iorder ) { // parfix=0 means parameter is free
rmin.Release( ipar );
somethingtodo = true;
}
}
scaleFitNotDone_ = false;
}
unsigned int crossSectionParShift = parResol.size()+parScale.size();
if( doCrossSectionFit[loopCounter] ) {
// Release cross section parameters and fit them
// ---------------------------------------------
// Note that only cross sections of resonances that are being fitted are released
bool doCrossSection = crossSectionHandler->releaseParameters( rmin, resfind, parfix, ind, iorder, crossSectionParShift );
if( doCrossSection ) somethingtodo = true;
}
if( doBackgroundFit[loopCounter] ) {
// Release background parameters and fit them
// ------------------------------------------
// for( int ipar=parResol.size()+parScale.size(); ipar<parnumber; ++ipar ) {
// Free only the parameters for the regions, as the resonances intervals are never used to fit the background
unsigned int bgrParShift = crossSectionParShift+crossSectionHandler->parNum();
for( unsigned int ipar = bgrParShift; ipar < bgrParShift+backgroundHandler->regionsParNum(); ++ipar ) {
// Release only those parameters for the resonances we are fitting
if( parfix[ipar]==0 && ind[ipar]==iorder && backgroundHandler->unlockParameter(resfind, ipar - bgrParShift) ) {
rmin.Release( ipar );
somethingtodo = true;
}
}
}
// OK, now do minimization if some parameter has been released
// -----------------------------------------------------------
if( somethingtodo ) {
// #ifdef DEBUG
std::stringstream fileNum;
fileNum << loopCounter;
minuitLoop_ = 0;
char name[50];
sprintf(name, "likelihoodInLoop_%d_%d", loopCounter, iorder);
TH1D * tempLikelihoodInLoop = new TH1D(name, "likelihood value in minuit loop", 10000, 0, 10000);
likelihoodInLoop_ = tempLikelihoodInLoop;
char signalProbName[50];
sprintf(signalProbName, "signalProb_%d_%d", loopCounter, iorder);
TH1D * tempSignalProb = new TH1D(signalProbName, "signal probability", 10000, 0, 10000);
signalProb_ = tempSignalProb;
char backgroundProbName[50];
sprintf(backgroundProbName, "backgroundProb_%d_%d", loopCounter, iorder);
TH1D * tempBackgroundProb = new TH1D(backgroundProbName, "background probability", 10000, 0, 10000);
backgroundProb_ = tempBackgroundProb;
// #endif
// Before we start the minimization we create a list of events with only the events inside a smaller
// window than the one in which the probability is != 0. We will compute the probability for all those
// events and hopefully the margin will avoid them to get a probability = 0 (which causes a discontinuity
// in the likelihood function). The width of this smaller window must be optimized, but we can start using
// an 90% of the normalization window.
double protectionFactor = 0.9;
MuScleFitUtils::ReducedSavedPair.clear();
for( unsigned int nev=0; nev<MuScleFitUtils::SavedPair.size(); ++nev ) {
const lorentzVector * recMu1 = &(MuScleFitUtils::SavedPair[nev].first);
const lorentzVector * recMu2 = &(MuScleFitUtils::SavedPair[nev].second);
double mass = MuScleFitUtils::invDimuonMass( *recMu1, *recMu2 );
// Test all resonances to see if the mass is inside at least one of the windows
bool check = false;
for( int ires = 0; ires < 6; ++ires ) {
// std::pair<double, double> windowFactor = backgroundHandler->windowFactors( doBackgroundFit[loopCounter], ires );
std::pair<double, double> windowBorder = backgroundHandler->windowBorders( doBackgroundFit[loopCounter], ires );
// if( resfind[ires] && checkMassWindow( mass, ires, backgroundHandler->resMass( doBackgroundFit[loopCounter], ires ),
// 0.9*windowFactor.first, 0.9*windowFactor.second ) ) {
double resMassValue = backgroundHandler->resMass( doBackgroundFit[loopCounter], ires );
double windowBorderLeft = resMassValue - protectionFactor*(resMassValue - windowBorder.first);
double windowBorderRight = resMassValue + protectionFactor*(windowBorder.second - resMassValue);
if( resfind[ires] && checkMassWindow( mass, windowBorderLeft, windowBorderRight ) ) {
check = true;
}
}
if( check ) {
MuScleFitUtils::ReducedSavedPair.push_back(std::make_pair(*recMu1, *recMu2));
}
}
std::cout << "Fitting with " << MuScleFitUtils::ReducedSavedPair.size() << " events" << std::endl;
// rmin.SetMaxIterations(500*parnumber);
MuScleFitUtils::normalizationChanged_ = 0;
// Maximum number of iterations
arglis[0] = 5000;
// Run simplex first to get an initial estimate of the minimum
if( startWithSimplex_ ) {
rmin.mnexcm( "simplex", arglis, 0, ierror );
}
rmin.mnexcm( "mini", arglis, 0, ierror );
// #ifdef DEBUG
likelihoodInLoop_->Write();
signalProb_->Write();
backgroundProb_->Write();
delete tempLikelihoodInLoop;
delete tempSignalProb;
delete tempBackgroundProb;
likelihoodInLoop_ = 0;
signalProb_ = 0;
backgroundProb_ = 0;
// #endif
// Compute again the error matrix
rmin.mnexcm( "hesse", arglis, 0, ierror );
// Peform minos error analysis.
if( computeMinosErrors_ ) {
duringMinos_ = true;
rmin.mnexcm( "minos", arglis, 0, ierror );
duringMinos_ = false;
}
if( normalizationChanged_ > 1 ) {
std::cout << "WARNING: normalization changed during fit meaning that events exited from the mass window. This causes a discontinuity in the likelihood function. Please check the scan of the likelihood as a function of the parameters to see if there are discontinuities around the minimum." << std::endl;
}
}
// bool notWritten = true;
for (int ipar=0; ipar<parnumber; ipar++) {
rmin.mnpout (ipar, name, pval, erro, pmin, pmax, ivar);
// Save parameters in parvalue[] vector
// ------------------------------------
if (ierror!=0 && debug>0) {
std::cout << "[MuScleFitUtils-minimizeLikelihood]: ierror!=0, bogus pars" << std::endl;
}
// for (int ipar=0; ipar<parnumber; ipar++) {
// rmin.mnpout (ipar, name, pval, erro, pmin, pmax, ivar);
parvalue[loopCounter][ipar] = pval;
// }
// int ilax2 = 0;
// Double_t val2pl, val2mi;
// rmin.mnmnot (ipar+1, ilax2, val2pl, val2mi);
rmin.mnerrs (ipar, errh, errl, errp, cglo);
// Set error on params
// -------------------
if (errp!=0) {
parerr[3*ipar] = errp;
} else {
parerr[3*ipar] = (((errh)>(fabs(errl)))?(errh):(fabs(errl)));
}
parerr[3*ipar+1] = errl;
parerr[3*ipar+2] = errh;
if( ipar == 0 ) {
FitParametersFile << " Resolution fit parameters:" << std::endl;
}
if( ipar == int(parResol.size()) ) {
FitParametersFile << " Scale fit parameters:" << std::endl;
}
if( ipar == int(parResol.size()+parScale.size()) ) {
FitParametersFile << " Cross section fit parameters:" << std::endl;
}
if( ipar == int(parResol.size()+parScale.size()+crossSectionHandler->parNum()) ) {
FitParametersFile << " Background fit parameters:" << std::endl;
}
// if( ipar >= int(parResol.size()+parScale.size()) && ipar < int(parResol.size()+parScale.size()+crossSectionHandler->parNum()) && notWritted ) {
// std::vector<double> relativeCrossSections = crossSectionHandler->relativeCrossSections(&(parvalue[loopCounter][parResol.size()+parScale.size()]));
// std::vector<double>::const_iterator it = relativeCrossSections.begin();
// for( ; it != relativeCrossSections.end(); ++it ) {
// FitParametersFile << " Results of the fit: parameter " << ipar << " has value "
// << *it << "+-" << 0
// << " + " << 0 << " - " << 0
// << " /t/t (" << 0 << ")" << std::endl;
// }
// notWritten = false;
// }
// else {
FitParametersFile << " Results of the fit: parameter " << ipar << " has value "
<< pval << "+-" << parerr[3*ipar]
<< " + " << parerr[3*ipar+1] << " - " << parerr[3*ipar+2]
// << " \t\t (" << parname[ipar] << ")"
<< std::endl;
// }
}
rmin.mnstat (fmin, fdem, errdef, npari, nparx, istat); // NNBB Commented for a check!
FitParametersFile << std::endl;
if( minimumShapePlots_ ) {
// Create plots of the minimum vs parameters
// -----------------------------------------
// Keep this after the parameters filling because it recomputes the values and it can compromise the fit results.
if( somethingtodo ) {
std::stringstream iorderString;
iorderString << iorder;
std::stringstream iLoopString;
iLoopString << loopCounter;
for (int ipar=0; ipar<parnumber; ipar++) {
if( parfix[ipar] == 1 ) continue;
std::cout << "plotting parameter = " << ipar+1 << std::endl;
std::stringstream iparString;
iparString << ipar+1;
rmin.mncomd( ("scan "+iparString.str()).c_str(), ierror );
if( ierror == 0 ) {
TCanvas * canvas = new TCanvas(("likelihoodCanvas_loop_"+iLoopString.str()+"_oder_"+iorderString.str()+"_par_"+iparString.str()).c_str(), ("likelihood_"+iparString.str()).c_str(), 1000, 800);
canvas->cd();
// arglis[0] = ipar;
// rmin.mnexcm( "sca", arglis, 0, ierror );
TGraph * graph = (TGraph*)rmin.GetPlot();
graph->Draw("AP");
// graph->SetTitle(("parvalue["+iparString.str()+"]").c_str());
graph->SetTitle(parname[ipar]);
// graph->Write();
canvas->Write();
}
}
// // Draw contours of the fit
// TCanvas * canvas = new TCanvas(("contourCanvas_oder_"+iorderString.str()).c_str(), "contour", 1000, 800);
// canvas->cd();
// TGraph * contourGraph = (TGraph*)rmin.Contour(4, 2, 4);
// if( (rmin.GetStatus() == 0) || (rmin.GetStatus() >= 3) ) {
// contourGraph->Draw("AP");
// }
// else {
// std::cout << "Contour graph error: status = " << rmin.GetStatus() << std::endl;
// }
// canvas->Write();
}
}
} // end loop on iorder
FitParametersFile.close();
std::cout << "[MuScleFitUtils-minimizeLikelihood]: Parameters after likelihood " << std::endl;
for (unsigned int ipar=0; ipar<(unsigned int)parnumber; ipar++) {
std::cout << ipar << " " << parvalue[loopCounter][ipar] << " : free = "
<< parfix[ipar] << "; order = " << parorder[ipar] << std::endl;
}
// Put back parvalue into parResol, parScale, parCrossSection, parBgr
// ------------------------------------------------------------------
for( int i=0; i<(int)(parResol.size()); ++i ) {
parResol[i] = parvalue[loopCounter][i];
}
for( int i=0; i<(int)(parScale.size()); ++i ) {
parScale[i] = parvalue[loopCounter][i+parResol.size()];
}
parCrossSection = crossSectionHandler->relativeCrossSections(&(parvalue[loopCounter][parResol.size()+parScale.size()]), resfind);
for( unsigned int i=0; i<parCrossSection.size(); ++i ) {
// parCrossSection[i] = parvalue[loopCounter][i+parResol.size()+parScale.size()];
std::cout << "relative cross section["<<i<<"] = " << parCrossSection[i] << std::endl;
}
// Save only the fitted background parameters
for( unsigned int i = 0; i<(parBgr.size() - parForResonanceWindows); ++i ) {
parBgr[i] = parvalue[loopCounter][i+parResol.size()+parScale.size()+crossSectionHandler->parNum()];
}
// Background rescaling from regions to resonances
// -----------------------------------------------
// Only if we fitted the background
if( doBackgroundFit[loopCounter] ) {
// This rescales from regions to resonances
int localMuonType = MuonType;
if( MuonType > 2 ) localMuonType = 2;
backgroundHandler->rescale( parBgr, ResMass, massWindowHalfWidth[localMuonType],
MuScleFitUtils::ReducedSavedPair );
}
// Delete the arrays used to set some parameters
delete[] Start;
delete[] Step;
delete[] Mini;
delete[] Maxi;
delete[] ind;
delete[] parname;
}
| double MuScleFitUtils::probability | ( | const double & | mass, |
| const double & | massResol, | ||
| const double | GLvalue[][1001][1001], | ||
| const double | GLnorm[][1001], | ||
| const int | iRes, | ||
| const int | iY | ||
| ) | [static] |
Computes the probability given the mass, mass resolution and the arrays with the probabilities and the normalizations.
After the introduction of the rapidity bins for the Z the probability method works in the following way:
Definition at line 720 of file MuScleFitUtils.cc.
References counter_resprob, gather_cfg::cout, debug, LogDebug, nbins, ResHalfWidth, ResMass, ResMaxSigma, and ResMinMass.
Referenced by massProb(), and MuScleFitBase::ProbForIntegral::operator()().
{
if( iRes == 0 && iY > 23 ) {
std::cout << "WARNING: rapidity bin selected = " << iY << " but there are only histograms for the first 24 bins" << std::endl;
}
double PS = 0.;
bool insideProbMassWindow = true;
// Interpolate the four values of GLZValue[] in the
// grid square within which the (mass,sigma) values lay
// ----------------------------------------------------
// This must be done with respect to the width used in the computation of the probability distribution,
// so that the bin 0 really matches the bin 0 of that distribution.
// double fracMass = (mass-(ResMass[iRes]-ResHalfWidth[iRes]))/(2*ResHalfWidth[iRes]);
double fracMass = (mass - ResMinMass[iRes])/(2*ResHalfWidth[iRes]);
if (debug>1) std::cout << std::setprecision(9)<<"mass ResMinMass[iRes] ResHalfWidth[iRes] ResHalfWidth[iRes]"
<< mass << " "<<ResMinMass[iRes]<<" "<<ResHalfWidth[iRes]<<" "<<ResHalfWidth[iRes]<<std::endl;
int iMassLeft = (int)(fracMass*(double)nbins);
int iMassRight = iMassLeft+1;
double fracMassStep = (double)nbins*(fracMass - (double)iMassLeft/(double)nbins);
if (debug>1) std::cout<<"nbins iMassLeft fracMass "<<nbins<<" "<<iMassLeft<<" "<<fracMass<<std::endl;
// Simple protections for the time being: the region where we fit should not include
// values outside the boundaries set by ResMass-ResHalfWidth : ResMass+ResHalfWidth
// ---------------------------------------------------------------------------------
if (iMassLeft<0) {
edm::LogInfo("probability") << "WARNING: fracMass=" << fracMass << ", iMassLeft="
<< iMassLeft << "; mass = " << mass << " and bounds are " << ResMinMass[iRes]
<< ":" << ResMinMass[iRes]+2*ResHalfWidth[iRes] << " - iMassLeft set to 0" << std::endl;
iMassLeft = 0;
iMassRight = 1;
insideProbMassWindow = false;
}
if (iMassRight>nbins) {
edm::LogInfo("probability") << "WARNING: fracMass=" << fracMass << ", iMassRight="
<< iMassRight << "; mass = " << mass << " and bounds are " << ResMinMass[iRes]
<< ":" << ResMass[iRes]+2*ResHalfWidth[iRes] << " - iMassRight set to " << nbins-1 << std::endl;
iMassLeft = nbins-1;
iMassRight = nbins;
insideProbMassWindow = false;
}
double fracSigma = (massResol/ResMaxSigma[iRes]);
int iSigmaLeft = (int)(fracSigma*(double)nbins);
int iSigmaRight = iSigmaLeft+1;
double fracSigmaStep = (double)nbins * (fracSigma - (double)iSigmaLeft/(double)nbins);
// std::cout << "massResol = " << massResol << std::endl;
// std::cout << "ResMaxSigma["<<iRes<<"] = " << ResMaxSigma[iRes] << std::endl;
// std::cout << "fracSigma = " << fracSigma << std::endl;
// std::cout << "nbins = " << nbins << std::endl;
// std::cout << "ISIGMALEFT = " << iSigmaLeft << std::endl;
// std::cout << "ISIGMARIGHT = " << iSigmaRight << std::endl;
// std::cout << "fracSigmaStep = " << fracSigmaStep << std::endl;
// Simple protections for the time being: they should not affect convergence, since
// ResMaxSigma is set to very large values, and if massResol exceeds them the fit
// should not get any prize for that (for large sigma, the prob. distr. becomes flat)
// ----------------------------------------------------------------------------------
if (iSigmaLeft<0) {
edm::LogInfo("probability") << "WARNING: fracSigma = " << fracSigma << ", iSigmaLeft="
<< iSigmaLeft << ", with massResol = " << massResol << " and ResMaxSigma[iRes] = "
<< ResMaxSigma[iRes] << " - iSigmaLeft set to 0" << std::endl;
iSigmaLeft = 0;
iSigmaRight = 1;
}
if (iSigmaRight>nbins ) {
if (counter_resprob<100)
edm::LogInfo("probability") << "WARNING: fracSigma = " << fracSigma << ", iSigmaRight="
<< iSigmaRight << ", with massResol = " << massResol << " and ResMaxSigma[iRes] = "
<< ResMaxSigma[iRes] << " - iSigmaRight set to " << nbins-1 << std::endl;
iSigmaLeft = nbins-1;
iSigmaRight = nbins;
}
// If f11,f12,f21,f22 are the values at the four corners, one finds by linear interpolation the
// formula below for PS
// --------------------------------------------------------------------------------------------
if( insideProbMassWindow ) {
double f11 = 0.;
if (GLnorm[iY][iSigmaLeft]>0)
f11 = GLvalue[iY][iMassLeft][iSigmaLeft] / GLnorm[iY][iSigmaLeft];
double f12 = 0.;
if (GLnorm[iY][iSigmaRight]>0)
f12 = GLvalue[iY][iMassLeft][iSigmaRight] / GLnorm[iY][iSigmaRight];
double f21 = 0.;
if (GLnorm[iY][iSigmaLeft]>0)
f21 = GLvalue[iY][iMassRight][iSigmaLeft] / GLnorm[iY][iSigmaLeft];
double f22 = 0.;
if (GLnorm[iY][iSigmaRight]>0)
f22 = GLvalue[iY][iMassRight][iSigmaRight] / GLnorm[iY][iSigmaRight];
PS = f11 + (f12-f11)*fracSigmaStep + (f21-f11)*fracMassStep +
(f22-f21-f12+f11)*fracMassStep*fracSigmaStep;
if (PS>0.1 || debug>1) LogDebug("MuScleFitUtils") << "iRes = " << iRes << " PS=" << PS << " f11,f12,f21,f22="
<< f11 << " " << f12 << " " << f21 << " " << f22 << " "
<< " fSS=" << fracSigmaStep << " fMS=" << fracMassStep << " iSL, iSR="
<< iSigmaLeft << " " << iSigmaRight
<< " GLvalue["<<iY<<"]["<<iMassLeft<<"] = " << GLvalue[iY][iMassLeft][iSigmaLeft]
<< " GLnorm["<<iY<<"]["<<iSigmaLeft<<"] = " << GLnorm[iY][iSigmaLeft] << std::endl;
// if (PS>0.1) std::cout << "iRes = " << iRes << " PS=" << PS << " f11,f12,f21,f22="
// << f11 << " " << f12 << " " << f21 << " " << f22 << " "
// << " fSS=" << fracSigmaStep << " fMS=" << fracMassStep << " iSL, iSR="
// << iSigmaLeft << " " << iSigmaRight << " GLV,GLN="
// << GLvalue[iY][iMassLeft][iSigmaLeft]
// << " " << GLnorm[iY][iSigmaLeft] << std::endl;
}
else {
edm::LogInfo("probability") << "outside mass probability window. Setting PS["<<iRes<<"] = 0" << std::endl;
}
// if( PS != PS ) {
// std::cout << "ERROR: PS = " << PS << " for iRes = " << iRes << std::endl;
// std::cout << "mass = " << mass << ", massResol = " << massResol << std::endl;
// std::cout << "fracMass = " << fracMass << ", iMassLeft = " << iMassLeft
// << ", iMassRight = " << iMassRight << ", fracMassStep = " << fracMassStep << std::endl;
// std::cout << "fracSigma = " << fracSigma << ", iSigmaLeft = " << iSigmaLeft
// << ", iSigmaRight = " << iSigmaRight << ", fracSigmaStep = " << fracSigmaStep << std::endl;
// std::cout << "ResMaxSigma["<<iRes<<"] = " << ResMaxSigma[iRes] << std::endl;
// std::cout << "GLvalue["<<iY<<"]["<<iMassLeft<<"] = " << GLvalue[iY][iMassLeft][iSigmaLeft]
// << " GLnorm["<<iY<<"]["<<iSigmaLeft<<"] = " << GLnorm[iY][iSigmaLeft] << std::endl;
// }
return PS;
}
const int MuScleFitUtils::backgroundFunctionsRegions [static] |
Definition at line 141 of file MuScleFitUtils.h.
Definition at line 153 of file MuScleFitUtils.h.
Referenced by computeWeight(), massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
TH1D * MuScleFitUtils::backgroundProb_ = 0 [static] |
Definition at line 164 of file MuScleFitUtils.h.
Referenced by massProb(), and minimizeLikelihood().
int MuScleFitUtils::BgrFitType = 0 [static] |
Definition at line 136 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood().
scaleFunctionBase< std::vector< double > > * MuScleFitUtils::biasFunction = 0 [static] |
Definition at line 129 of file MuScleFitUtils.h.
Referenced by applyBias(), and MuScleFit::MuScleFit().
int MuScleFitUtils::BiasType = 0 [static] |
Definition at line 127 of file MuScleFitUtils.h.
Referenced by MuScleFit::applyBias(), MuScleFit::checkParameters(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
bool MuScleFitUtils::computeMinosErrors_ [static] |
Definition at line 248 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
int MuScleFitUtils::counter_resprob = 0 [static] |
Definition at line 187 of file MuScleFitUtils.h.
Referenced by massResolution(), probability(), and MuScleFit::startingNewLoop().
double MuScleFitUtils::crossSection[6] [static] |
Definition at line 116 of file MuScleFitUtils.h.
Definition at line 149 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
int MuScleFitUtils::debug = 0 [static] |
Definition at line 108 of file MuScleFitUtils.h.
Referenced by applyBias(), applyScale(), applySmearing(), computeWeight(), ErrorsAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::fillHistograms(), findBestRecoRes(), findGenMuFromRes(), fitMass(), likelihood(), massProb(), massResolution(), minimizeLikelihood(), MuScleFit::MuScleFit(), and probability().
bool MuScleFitUtils::debugMassResol_ [static] |
Definition at line 235 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), ErrorsAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::fillHistograms(), MuScleFitBase::fillHistoMap(), massResolution(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::doBackgroundFit [static] |
Definition at line 159 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), computeWeight(), massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::doCrossSectionFit [static] |
Definition at line 158 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::doResolFit [static] |
Definition at line 156 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::doScaleFit [static] |
Definition at line 157 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), MuScleFit::duringFastLoop(), likelihood(), minimizeLikelihood(), and MuScleFit::MuScleFit().
bool MuScleFitUtils::duringMinos_ = false [static] |
Definition at line 166 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood().
int MuScleFitUtils::FitStrategy = 1 [static] |
Definition at line 183 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::genPair [static] |
Definition at line 205 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), and MuScleFit::selectMuons().
double MuScleFitUtils::GLNorm [static] |
Definition at line 191 of file MuScleFitUtils.h.
Referenced by massProb(), MuScleFitBase::ProbForIntegral::operator()(), and MuScleFitBase::readProbabilityDistributionsFromFile().
double MuScleFitUtils::GLValue [static] |
Definition at line 190 of file MuScleFitUtils.h.
Referenced by massProb(), MuScleFitBase::ProbForIntegral::operator()(), and MuScleFitBase::readProbabilityDistributionsFromFile().
double MuScleFitUtils::GLZNorm [static] |
Definition at line 189 of file MuScleFitUtils.h.
Referenced by massProb(), MuScleFitBase::ProbForIntegral::operator()(), and MuScleFitBase::readProbabilityDistributionsFromFile().
double MuScleFitUtils::GLZValue [static] |
Definition at line 188 of file MuScleFitUtils.h.
Referenced by massProb(), MuScleFitBase::ProbForIntegral::operator()(), and MuScleFitBase::readProbabilityDistributionsFromFile().
int MuScleFitUtils::goodmuon = 0 [static] |
Definition at line 186 of file MuScleFitUtils.h.
Referenced by MuScleFit::applyBias(), MuScleFit::applySmearing(), applySmearing(), MuScleFit::fillMuonCollection(), ResolutionAnalyzer::fillMuonCollection(), MuScleFit::MuScleFit(), and MuScleFit::startingNewLoop().
int MuScleFitUtils::iev_ = 0 [static] |
Definition at line 223 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), likelihood(), and MuScleFit::startingNewLoop().
TH1D * MuScleFitUtils::likelihoodInLoop_ = 0 [static] |
Definition at line 162 of file MuScleFitUtils.h.
Referenced by likelihood(), and minimizeLikelihood().
unsigned int MuScleFitUtils::loopCounter = 5 [static] |
Definition at line 123 of file MuScleFitUtils.h.
Referenced by computeWeight(), likelihood(), massProb(), minimizeLikelihood(), and MuScleFit::startingNewLoop().
Definition at line 243 of file MuScleFitUtils.cc.
Referenced by MuScleFit::duringFastLoop(), and massResolution().
double MuScleFitUtils::massWindowHalfWidth [static] |
Definition at line 112 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
double MuScleFitUtils::maxMuonEtaFirstRange_ = 6. [static] |
Definition at line 231 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
double MuScleFitUtils::maxMuonEtaSecondRange_ = 100. [static] |
Definition at line 233 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
double MuScleFitUtils::maxMuonPt_ = 100000000. [static] |
Definition at line 229 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
bool MuScleFitUtils::minimumShapePlots_ [static] |
Definition at line 249 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
double MuScleFitUtils::minMuonEtaFirstRange_ = -6. [static] |
Definition at line 230 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
double MuScleFitUtils::minMuonEtaSecondRange_ = -100. [static] |
Definition at line 232 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
double MuScleFitUtils::minMuonPt_ = 0. [static] |
Definition at line 228 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
int MuScleFitUtils::minuitLoop_ = 0 [static] |
Definition at line 161 of file MuScleFitUtils.h.
Referenced by likelihood(), massProb(), and minimizeLikelihood().
const double MuScleFitUtils::mMu2 = 0.011163612 [static] |
Definition at line 117 of file MuScleFitUtils.h.
Referenced by ResolutionAnalyzer::analyze(), MuScleFit::fillMuonCollection(), ResolutionAnalyzer::fillMuonCollection(), and massResolution().
const unsigned int MuScleFitUtils::motherPdgIdArray = {23, 100553, 100553, 553, 100443, 443} [static] |
Definition at line 121 of file MuScleFitUtils.h.
Referenced by findGenMuFromRes(), and findSimMuFromRes().
const double MuScleFitUtils::muMass = 0.105658 [static] |
Definition at line 118 of file MuScleFitUtils.h.
Referenced by fromPtEtaPhiToPxPyPz().
int MuScleFitUtils::MuonType [static] |
Definition at line 195 of file MuScleFitUtils.h.
Referenced by massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
int MuScleFitUtils::MuonTypeForCheckMassWindow [static] |
Definition at line 196 of file MuScleFitUtils.h.
Referenced by MuScleFit::MuScleFit().
int MuScleFitUtils::nbins = 1000 [static] |
Definition at line 194 of file MuScleFitUtils.h.
Referenced by probability(), and MuScleFitBase::readProbabilityDistributionsFromFile().
unsigned int MuScleFitUtils::normalizationChanged_ = 0 [static] |
Definition at line 215 of file MuScleFitUtils.h.
Referenced by likelihood(), and minimizeLikelihood().
bool MuScleFitUtils::normalizeLikelihoodByEventNumber_ = true [static] |
Definition at line 210 of file MuScleFitUtils.h.
Referenced by likelihood(), and MuScleFit::MuScleFit().
double MuScleFitUtils::oldNormalization_ = 0. [static] |
Definition at line 214 of file MuScleFitUtils.h.
Referenced by likelihood(), and MuScleFit::startingNewLoop().
std::vector< double > MuScleFitUtils::parBgr [static] |
Definition at line 173 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), MuScleFit::duringFastLoop(), massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parBgrFix [static] |
Definition at line 177 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parBgrOrder [static] |
Definition at line 181 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parBias [static] |
Definition at line 169 of file MuScleFitUtils.h.
Referenced by applyBias(), MuScleFit::checkParameters(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parCrossSection [static] |
Definition at line 172 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), MuScleFit::duringFastLoop(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parCrossSectionFix [static] |
Definition at line 176 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parCrossSectionOrder [static] |
Definition at line 180 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector<int> MuScleFitUtils::parfix [static] |
Definition at line 200 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood().
std::vector<int> MuScleFitUtils::parorder [static] |
Definition at line 201 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood().
std::vector< double > MuScleFitUtils::parResol [static] |
Definition at line 170 of file MuScleFitUtils.h.
Referenced by ResolutionAnalyzer::analyze(), applyScale(), MuScleFit::checkParameters(), MuScleFit::duringFastLoop(), findBestRecoRes(), massProb(), minimizeLikelihood(), MuScleFit::MuScleFit(), ResolutionAnalyzer::ResolutionAnalyzer(), and TestCorrection::TestCorrection().
std::vector< int > MuScleFitUtils::parResolFix [static] |
Definition at line 174 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parResolOrder [static] |
Definition at line 178 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parScale [static] |
Definition at line 171 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), MuScleFit::duringFastLoop(), massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parScaleFix [static] |
Definition at line 175 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parScaleOrder [static] |
Definition at line 179 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parSmear [static] |
Definition at line 168 of file MuScleFitUtils.h.
Referenced by applySmearing(), MuScleFit::checkParameters(), and MuScleFit::MuScleFit().
std::vector< std::vector< double > > MuScleFitUtils::parvalue [static] |
Definition at line 198 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), and minimizeLikelihood().
bool MuScleFitUtils::rapidityBinsForZ_ = true [static] |
Definition at line 221 of file MuScleFitUtils.h.
Referenced by massProb(), and MuScleFit::MuScleFit().
std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::ReducedSavedPair [static] |
Definition at line 204 of file MuScleFitUtils.h.
Referenced by likelihood(), and minimizeLikelihood().
std::vector< int > MuScleFitUtils::resfind [static] |
Definition at line 182 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), computeWeight(), ErrorsAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::fillHistograms(), ResolutionAnalyzer::fillHistoMap(), MuScleFitBase::fillHistoMap(), findBestRecoRes(), findBestSimuRes(), findGenMuFromRes(), findSimMuFromRes(), massProb(), massResolution(), minimizeLikelihood(), MuScleFit::MuScleFit(), MuScleFitGenFilter::MuScleFitGenFilter(), MuScleFitBase::readProbabilityDistributionsFromFile(), ResolutionAnalyzer::ResolutionAnalyzer(), and TestCorrection::TestCorrection().
bool MuScleFitUtils::ResFound = false [static] |
Definition at line 109 of file MuScleFitUtils.h.
Referenced by TestCorrection::analyze(), MuScleFit::duringFastLoop(), findBestRecoRes(), and MuScleFit::selectMuons().
double MuScleFitUtils::ResGamma = {2.4952, 0.000020, 0.000032, 0.000054, 0.000317, 0.0000932 } [static] |
Definition at line 113 of file MuScleFitUtils.h.
Referenced by massProb().
double MuScleFitUtils::ResHalfWidth [static] |
Definition at line 193 of file MuScleFitUtils.h.
Referenced by massProb(), massResolution(), probability(), and MuScleFitBase::readProbabilityDistributionsFromFile().
double MuScleFitUtils::ResMass = {91.1876, 10.3552, 10.0233, 9.4603, 3.68609, 3.0969} [static] |
Definition at line 114 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), massProb(), massResolution(), minimizeLikelihood(), MuScleFit::MuScleFit(), and probability().
double MuScleFitUtils::ResMaxSigma [static] |
Definition at line 192 of file MuScleFitUtils.h.
Referenced by massResolution(), probability(), and MuScleFitBase::readProbabilityDistributionsFromFile().
double MuScleFitUtils::ResMinMass = {-99, -99, -99, -99, -99, -99} [static] |
Definition at line 115 of file MuScleFitUtils.h.
Referenced by probability(), and MuScleFitBase::readProbabilityDistributionsFromFile().
int MuScleFitUtils::ResolFitType = 0 [static] |
Definition at line 130 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), MuScleFit::MuScleFit(), and ResolutionAnalyzer::ResolutionAnalyzer().
resolutionFunctionBase< double * > * MuScleFitUtils::resolutionFunction = 0 [static] |
Definition at line 131 of file MuScleFitUtils.h.
Referenced by ErrorsAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::fillHistograms(), massResolution(), MuScleFit::MuScleFit(), ResolutionAnalyzer::ResolutionAnalyzer(), and TestCorrection::TestCorrection().
resolutionFunctionBase< std::vector< double > > * MuScleFitUtils::resolutionFunctionForVec = 0 [static] |
Definition at line 132 of file MuScleFitUtils.h.
Referenced by ResolutionAnalyzer::analyze(), MuScleFit::duringFastLoop(), MuScleFit::MuScleFit(), ResolutionAnalyzer::ResolutionAnalyzer(), and TestCorrection::TestCorrection().
TMinuit * MuScleFitUtils::rminPtr_ = 0 [static] |
Definition at line 212 of file MuScleFitUtils.h.
Referenced by likelihood(), and minimizeLikelihood().
std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::SavedPair [static] |
Definition at line 203 of file MuScleFitUtils.h.
Referenced by TestCorrection::analyze(), MuScleFit::duringFastLoop(), MuScleFit::endOfLoop(), likelihood(), minimizeLikelihood(), MuScleFit::selectMuons(), and MuScleFit::~MuScleFit().
bool MuScleFitUtils::scaleFitNotDone_ = true [static] |
Definition at line 208 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood().
int MuScleFitUtils::ScaleFitType = 0 [static] |
Definition at line 133 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
scaleFunctionBase< double * > * MuScleFitUtils::scaleFunction = 0 [static] |
Definition at line 134 of file MuScleFitUtils.h.
Referenced by applyScale(), minimizeLikelihood(), and MuScleFit::MuScleFit().
scaleFunctionBase< std::vector< double > > * MuScleFitUtils::scaleFunctionForVec = 0 [static] |
Definition at line 135 of file MuScleFitUtils.h.
Referenced by MuScleFit::MuScleFit().
bool MuScleFitUtils::sherpa_ = false [static] |
Definition at line 218 of file MuScleFitUtils.h.
Referenced by findGenMuFromRes(), and MuScleFit::MuScleFit().
TH1D * MuScleFitUtils::signalProb_ = 0 [static] |
Definition at line 163 of file MuScleFitUtils.h.
Referenced by massProb(), and minimizeLikelihood().
std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::simPair [static] |
Definition at line 206 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), and MuScleFit::selectMuons().
smearFunctionBase * MuScleFitUtils::smearFunction = 0 [static] |
Definition at line 126 of file MuScleFitUtils.h.
Referenced by applySmearing(), and MuScleFit::MuScleFit().
int MuScleFitUtils::SmearType = 0 [static] |
Definition at line 125 of file MuScleFitUtils.h.
Referenced by MuScleFit::applySmearing(), applySmearing(), MuScleFit::checkParameters(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
bool MuScleFitUtils::speedup = false [static] |
Definition at line 184 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), MuScleFit::MuScleFit(), MuScleFit::selectMuons(), and MuScleFit::~MuScleFit().
bool MuScleFitUtils::startWithSimplex_ [static] |
Definition at line 247 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
const int MuScleFitUtils::totalResNum = 6 [static] |
Definition at line 111 of file MuScleFitUtils.h.
Referenced by massProb().
bool MuScleFitUtils::useProbsFile_ = true [static] |
Definition at line 225 of file MuScleFitUtils.h.
Referenced by MuScleFit::beginOfJobInConstructor(), findBestRecoRes(), and MuScleFit::MuScleFit().
double MuScleFitUtils::x [static] |
Definition at line 185 of file MuScleFitUtils.h.
Referenced by applySmearing(), fitMass(), fitReso(), massProb(), and MuScleFit::MuScleFit().
1.7.3