#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, const double &eta1, const double &eta2) |
static double | massProb (const double &mass, const double &resEta, const double &rapidity, const double &massResol, double *parval, const bool doUseBkgrWindow, const double &eta1, const double &eta2) |
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, const ResolutionFunction &resolFunc) |
static double | massResolution (const lorentzVector &mu1, const lorentzVector &mu2, double *parval) |
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 double | deltaPhiMaxCut_ = 100. |
static double | deltaPhiMinCut_ = -100. |
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< double > | parResolMax |
static std::vector< double > | parResolMin |
static std::vector< int > | parResolOrder |
static std::vector< double > | parResolStep |
static std::vector< double > | parScale |
static std::vector< int > | parScaleFix |
static std::vector< double > | parScaleMax |
static std::vector< double > | parScaleMin |
static std::vector< int > | parScaleOrder |
static std::vector< double > | parScaleStep |
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 | separateRanges_ = true |
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 47 of file MuScleFitUtils.h.
MuScleFitUtils::MuScleFitUtils | ( | ) | [inline] |
Definition at line 52 of file MuScleFitUtils.h.
{};
virtual MuScleFitUtils::~MuScleFitUtils | ( | ) | [inline, virtual] |
Definition at line 56 of file MuScleFitUtils.h.
{};
lorentzVector MuScleFitUtils::applyBias | ( | const lorentzVector & | muon, |
const int | charge | ||
) | [static] |
Definition at line 430 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 450 of file MuScleFitUtils.cc.
References errorMatrix2Lands_multiChannel::id, and 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 470 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 400 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 1087 of file MuScleFitUtils.cc.
Referenced by computeWeight(), massProb(), and minimizeLikelihood().
double MuScleFitUtils::computeWeight | ( | const double & | mass, |
const int | iev, | ||
const bool | doUseBkgrWindow = false |
||
) | [static] |
Definition at line 1094 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 92 of file MuScleFitUtils.h.
References Pi.
Referenced by ResolutionAnalyzer::checkDeltaR(), MuScleFit::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 100 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 107 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 306 of file MuScleFitUtils.cc.
References gather_cfg::cout, debug, ires, scaleCards::mass, 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) { //rc2010 if (debug>0) std::cout << "muon_1_charge:"<<(*Muon1).charge() << std::endl; for (std::vector<reco::LeafCandidate>::const_iterator Muon2=Muon1+1; Muon2!=muons.end(); ++Muon2) { //rc2010 if (debug>0) std::cout << "after_2" << std::endl; 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 269 of file MuScleFitUtils.cc.
References ires, scaleCards::mass, 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 2217 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 2255 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 2179 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 1898 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 2011 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 489 of file MuScleFitUtils.cc.
References funct::cos(), create_public_lumi_plots::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 506 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 2121 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, | ||
const double & | eta1, | ||
const double & | eta2 | ||
) | [static] |
Definition at line 695 of file MuScleFitUtils.cc.
References errorMatrix2Lands_multiChannel::id, 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, eta1, eta2 ); 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, | ||
const double & | eta1, | ||
const double & | eta2 | ||
) | [static] |
Definition at line 860 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 ); std::pair<double, double> bgrResult = backgroundHandler->backgroundFunction( doBackgroundFit[loopCounter], &(parval[bgrParShift]), MuScleFitUtils::totalResNum, 0, resConsidered, ResMass, ResHalfWidth, MuonType, mass, eta1, eta2 ); 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 ); resConsidered, ResMass, ResHalfWidth, MuonType, mass, eta1, eta2 ); 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, | ||
const std::vector< double > & | parval | ||
) | [static] |
Definition at line 514 of file MuScleFitUtils.cc.
References errorMatrix2Lands_multiChannel::id, 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.: http://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] |
double MuScleFitUtils::massResolution | ( | const lorentzVector & | mu1, |
const lorentzVector & | mu2, | ||
const ResolutionFunction & | resolFunc | ||
) | [static] |
This method can be used outside MuScleFit. It gets the ResolutionFunction that must have been built with the parameters.
TO-DO: this method duplicates the code in the previous method. It should be changed to avoid the duplication.
Definition at line 647 of file MuScleFitUtils.cc.
References funct::cos(), create_public_lumi_plots::exp, scaleCards::mass, mMu2, funct::pow(), ResolutionFunction::sigmaCotgTh(), ResolutionFunction::sigmaPhi(), ResolutionFunction::sigmaPt(), funct::sin(), and mathSSE::sqrt().
{ double mass = (mu1+mu2).mass(); double pt1 = mu1.Pt(); double phi1 = mu1.Phi(); double eta1 = mu1.Eta(); double theta1 = 2*atan(exp(-eta1)); double pt2 = mu2.Pt(); double phi2 = mu2.Phi(); double eta2 = mu2.Eta(); double theta2 = 2*atan(exp(-eta2)); 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 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; // Resolution parameters: // ---------------------- double sigma_pt1 = resolFunc.sigmaPt( mu1 ); double sigma_pt2 = resolFunc.sigmaPt( mu2 ); double sigma_phi1 = resolFunc.sigmaPhi( mu1 ); double sigma_phi2 = resolFunc.sigmaPhi( mu2 ); double sigma_cotgth1 = resolFunc.sigmaCotgTh( mu1 ); double sigma_cotgth2 = resolFunc.sigmaCotgTh( mu2 ); // 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)); return mass_res; }
double MuScleFitUtils::massResolution | ( | const lorentzVector & | mu1, |
const lorentzVector & | mu2, | ||
double * | parval | ||
) | [static] |
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.
Definition at line 558 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, create_public_lumi_plots::exp, ires, LogDebug, scaleCards::mass, massResolComponents, mMu2, funct::pow(), resfind, ResHalfWidth, ResMass, ResMaxSigma, resolutionFunction, resolutionFunctionBase< T >::sigmaCotgTh(), resolutionFunctionBase< T >::sigmaPhi(), resolutionFunctionBase< T >::sigmaPt(), funct::sin(), and mathSSE::sqrt().
{ double mass = (mu1+mu2).mass(); double pt1 = mu1.Pt(); double phi1 = mu1.Phi(); double eta1 = mu1.Eta(); double theta1 = 2*atan(exp(-eta1)); double pt2 = mu2.Pt(); double phi2 = mu2.Phi(); double eta2 = mu2.Eta(); double theta2 = 2*atan(exp(-eta2)); 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 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*sigma_pt1*sigma_pt2); 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; didit = true; } } } return mass_res; }
static double MuScleFitUtils::massResolution | ( | const lorentzVector & | mu1, |
const lorentzVector & | mu2, | ||
std::auto_ptr< double > | parval | ||
) | [static] |
void MuScleFitUtils::minimizeLikelihood | ( | ) | [static] |
Definition at line 1127 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, scaleCards::mass, massWindowHalfWidth, minimumShapePlots_, minuitLoop_, MuonType, mergeVDriftHistosByStation::name, normalizationChanged_, parBgr, parBgrFix, parBgrOrder, parCrossSection, parCrossSectionOrder, parfix, scaleFunctionBase< T >::parNum(), CrossSectionHandler::parNum(), resolutionFunctionBase< T >::parNum(), parorder, parResol, parResolFix, parResolMax, parResolMin, parResolOrder, parResolStep, parScale, parScaleFix, parScaleMax, parScaleMin, parScaleOrder, parScaleStep, parvalue, ReducedSavedPair, BackgroundHandler::regionsParNum(), CrossSectionHandler::relativeCrossSections(), CrossSectionHandler::releaseParameters(), BackgroundHandler::rescale(), scaleFunctionBase< T >::resetParameters(), resfind, 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]; if( !parResolStep.empty() && !parResolMin.empty() && !parResolMax.empty() ) { MuScleFitUtils::resolutionFunctionForVec->setParameters( Start, Step, Mini, Maxi, ind, parname, parResol, parResolOrder, parResolStep, parResolMin, parResolMax, MuonType ); } else { 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(); if( !parScaleStep.empty() && !parScaleMin.empty() && !parScaleMax.empty() ) { MuScleFitUtils::scaleFunctionForVec->setParameters( &(Start[resParNum]), &(Step[resParNum]), &(Mini[resParNum]), &(Maxi[resParNum]), &(ind[resParNum]), &(parname[resParNum]), parScale, parScaleOrder, parScaleStep, parScaleMin, parScaleMax, MuonType ); } else { 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); double windowBorderShift = (windowBorder.second - windowBorder.first)*(1-protectionFactor)/2.; double windowBorderLeft = windowBorder.first + windowBorderShift; double windowBorderRight = windowBorder.second - windowBorderShift; 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); //Print some informations std::cout<<"MINUIT is starting the minimization for the iteration number "<<loopCounter<<std::endl; //Try to set iterations // rmin.SetMaxIterations(100000); std::cout<<"maxNumberOfIterations (just set) = "<<rmin.GetMaxIterations()<<std::endl; MuScleFitUtils::normalizationChanged_ = 0; // Maximum number of iterations arglis[0] = 100000; // tolerance arglis[1] = 0.1; // Run simplex first to get an initial estimate of the minimum if( startWithSimplex_ ) { rmin.mnexcm( "SIMPLEX", arglis, 0, ierror ); } rmin.mnexcm( "MIGRAD", arglis, 2, 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; std::stringstream iparStringName; iparStringName << ipar; rmin.mncomd( ("scan "+iparString.str()).c_str(), ierror ); if( ierror == 0 ) { TCanvas * canvas = new TCanvas(("likelihoodCanvas_loop_"+iLoopString.str()+"_oder_"+iorderString.str()+"_par_"+iparStringName.str()).c_str(), ("likelihood_"+iparStringName.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["+iparStringName.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 729 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 145 of file MuScleFitUtils.h.
Definition at line 157 of file MuScleFitUtils.h.
Referenced by computeWeight(), massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
TH1D * MuScleFitUtils::backgroundProb_ = 0 [static] |
Definition at line 168 of file MuScleFitUtils.h.
Referenced by massProb(), and minimizeLikelihood().
int MuScleFitUtils::BgrFitType = 0 [static] |
Definition at line 140 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood().
scaleFunctionBase< std::vector< double > > * MuScleFitUtils::biasFunction = 0 [static] |
Definition at line 133 of file MuScleFitUtils.h.
Referenced by applyBias(), and MuScleFit::MuScleFit().
int MuScleFitUtils::BiasType = 0 [static] |
Definition at line 131 of file MuScleFitUtils.h.
Referenced by MuScleFit::applyBias(), MuScleFit::checkParameters(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
bool MuScleFitUtils::computeMinosErrors_ [static] |
Definition at line 261 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
int MuScleFitUtils::counter_resprob = 0 [static] |
Definition at line 197 of file MuScleFitUtils.h.
Referenced by massResolution(), probability(), and MuScleFit::startingNewLoop().
double MuScleFitUtils::crossSection[6] [static] |
Definition at line 120 of file MuScleFitUtils.h.
Definition at line 153 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
int MuScleFitUtils::debug = 0 [static] |
Definition at line 112 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 248 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), ErrorsAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::fillHistograms(), MuScleFitBase::fillHistoMap(), massResolution(), and MuScleFit::MuScleFit().
double MuScleFitUtils::deltaPhiMaxCut_ = 100. [static] |
Definition at line 246 of file MuScleFitUtils.h.
Referenced by MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
double MuScleFitUtils::deltaPhiMinCut_ = -100. [static] |
Definition at line 245 of file MuScleFitUtils.h.
Referenced by MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
std::vector< int > MuScleFitUtils::doBackgroundFit [static] |
Definition at line 163 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), computeWeight(), massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::doCrossSectionFit [static] |
Definition at line 162 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::doResolFit [static] |
Definition at line 160 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::doScaleFit [static] |
Definition at line 161 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), MuScleFit::duringFastLoop(), likelihood(), minimizeLikelihood(), and MuScleFit::MuScleFit().
bool MuScleFitUtils::duringMinos_ = false [static] |
Definition at line 170 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood().
int MuScleFitUtils::FitStrategy = 1 [static] |
Definition at line 193 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::genPair [static] |
Definition at line 215 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), and MuScleFit::selectMuons().
double MuScleFitUtils::GLNorm [static] |
Definition at line 201 of file MuScleFitUtils.h.
Referenced by massProb(), MuScleFitBase::ProbForIntegral::operator()(), and MuScleFitBase::readProbabilityDistributionsFromFile().
double MuScleFitUtils::GLValue [static] |
Definition at line 200 of file MuScleFitUtils.h.
Referenced by massProb(), MuScleFitBase::ProbForIntegral::operator()(), and MuScleFitBase::readProbabilityDistributionsFromFile().
double MuScleFitUtils::GLZNorm [static] |
Definition at line 199 of file MuScleFitUtils.h.
Referenced by massProb(), MuScleFitBase::ProbForIntegral::operator()(), and MuScleFitBase::readProbabilityDistributionsFromFile().
double MuScleFitUtils::GLZValue [static] |
Definition at line 198 of file MuScleFitUtils.h.
Referenced by massProb(), MuScleFitBase::ProbForIntegral::operator()(), and MuScleFitBase::readProbabilityDistributionsFromFile().
int MuScleFitUtils::goodmuon = 0 [static] |
Definition at line 196 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 233 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), likelihood(), and MuScleFit::startingNewLoop().
TH1D * MuScleFitUtils::likelihoodInLoop_ = 0 [static] |
Definition at line 166 of file MuScleFitUtils.h.
Referenced by likelihood(), and minimizeLikelihood().
unsigned int MuScleFitUtils::loopCounter = 5 [static] |
Definition at line 127 of file MuScleFitUtils.h.
Referenced by computeWeight(), likelihood(), massProb(), minimizeLikelihood(), and MuScleFit::startingNewLoop().
Definition at line 252 of file MuScleFitUtils.cc.
Referenced by MuScleFit::duringFastLoop(), and massResolution().
double MuScleFitUtils::massWindowHalfWidth [static] |
Definition at line 116 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
double MuScleFitUtils::maxMuonEtaFirstRange_ = 6. [static] |
Definition at line 242 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
double MuScleFitUtils::maxMuonEtaSecondRange_ = 100. [static] |
Definition at line 244 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
double MuScleFitUtils::maxMuonPt_ = 100000000. [static] |
Definition at line 240 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
bool MuScleFitUtils::minimumShapePlots_ [static] |
Definition at line 262 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
double MuScleFitUtils::minMuonEtaFirstRange_ = -6. [static] |
Definition at line 241 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
double MuScleFitUtils::minMuonEtaSecondRange_ = -100. [static] |
Definition at line 243 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
double MuScleFitUtils::minMuonPt_ = 0. [static] |
Definition at line 239 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
int MuScleFitUtils::minuitLoop_ = 0 [static] |
Definition at line 165 of file MuScleFitUtils.h.
Referenced by likelihood(), massProb(), and minimizeLikelihood().
const double MuScleFitUtils::mMu2 = 0.011163612 [static] |
Definition at line 121 of file MuScleFitUtils.h.
Referenced by ResolutionAnalyzer::analyze(), ResolutionAnalyzer::fillMuonCollection(), MuScleFit::fillMuonCollection(), massResolution(), and ErrorsPropagationAnalyzer::massResolution().
const unsigned int MuScleFitUtils::motherPdgIdArray = {23, 100553, 100553, 553, 100443, 443} [static] |
Definition at line 125 of file MuScleFitUtils.h.
Referenced by findGenMuFromRes(), and findSimMuFromRes().
const double MuScleFitUtils::muMass = 0.105658 [static] |
Definition at line 122 of file MuScleFitUtils.h.
Referenced by fromPtEtaPhiToPxPyPz().
int MuScleFitUtils::MuonType [static] |
Definition at line 205 of file MuScleFitUtils.h.
Referenced by massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
int MuScleFitUtils::MuonTypeForCheckMassWindow [static] |
Definition at line 206 of file MuScleFitUtils.h.
Referenced by MuScleFit::MuScleFit().
int MuScleFitUtils::nbins = 1000 [static] |
Definition at line 204 of file MuScleFitUtils.h.
Referenced by probability(), and MuScleFitBase::readProbabilityDistributionsFromFile().
unsigned int MuScleFitUtils::normalizationChanged_ = 0 [static] |
Definition at line 225 of file MuScleFitUtils.h.
Referenced by likelihood(), and minimizeLikelihood().
bool MuScleFitUtils::normalizeLikelihoodByEventNumber_ = true [static] |
Definition at line 220 of file MuScleFitUtils.h.
Referenced by likelihood(), and MuScleFit::MuScleFit().
double MuScleFitUtils::oldNormalization_ = 0. [static] |
Definition at line 224 of file MuScleFitUtils.h.
Referenced by likelihood(), and MuScleFit::startingNewLoop().
std::vector< double > MuScleFitUtils::parBgr [static] |
Definition at line 183 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), MuScleFit::duringFastLoop(), massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parBgrFix [static] |
Definition at line 187 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parBgrOrder [static] |
Definition at line 191 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parBias [static] |
Definition at line 173 of file MuScleFitUtils.h.
Referenced by applyBias(), MuScleFit::checkParameters(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parCrossSection [static] |
Definition at line 182 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), MuScleFit::duringFastLoop(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parCrossSectionFix [static] |
Definition at line 186 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parCrossSectionOrder [static] |
Definition at line 190 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector<int> MuScleFitUtils::parfix [static] |
Definition at line 210 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood().
std::vector<int> MuScleFitUtils::parorder [static] |
Definition at line 211 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood().
std::vector< double > MuScleFitUtils::parResol [static] |
Definition at line 174 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 184 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parResolMax [static] |
Definition at line 177 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parResolMin [static] |
Definition at line 176 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parResolOrder [static] |
Definition at line 188 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parResolStep [static] |
Definition at line 175 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parScale [static] |
Definition at line 178 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), MuScleFit::duringFastLoop(), massProb(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parScaleFix [static] |
Definition at line 185 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parScaleMax [static] |
Definition at line 181 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parScaleMin [static] |
Definition at line 180 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< int > MuScleFitUtils::parScaleOrder [static] |
Definition at line 189 of file MuScleFitUtils.h.
Referenced by MuScleFit::checkParameters(), minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parScaleStep [static] |
Definition at line 179 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
std::vector< double > MuScleFitUtils::parSmear [static] |
Definition at line 172 of file MuScleFitUtils.h.
Referenced by applySmearing(), MuScleFit::checkParameters(), and MuScleFit::MuScleFit().
std::vector< std::vector< double > > MuScleFitUtils::parvalue [static] |
Definition at line 208 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), and minimizeLikelihood().
bool MuScleFitUtils::rapidityBinsForZ_ = true [static] |
Definition at line 231 of file MuScleFitUtils.h.
Referenced by massProb(), and MuScleFit::MuScleFit().
std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::ReducedSavedPair [static] |
Definition at line 214 of file MuScleFitUtils.h.
Referenced by likelihood(), and minimizeLikelihood().
std::vector< int > MuScleFitUtils::resfind [static] |
Definition at line 192 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 113 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 117 of file MuScleFitUtils.h.
Referenced by massProb().
double MuScleFitUtils::ResHalfWidth [static] |
Definition at line 203 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 118 of file MuScleFitUtils.h.
Referenced by findBestRecoRes(), massProb(), massResolution(), minimizeLikelihood(), MuScleFit::MuScleFit(), and probability().
double MuScleFitUtils::ResMaxSigma [static] |
Definition at line 202 of file MuScleFitUtils.h.
Referenced by massResolution(), probability(), and MuScleFitBase::readProbabilityDistributionsFromFile().
double MuScleFitUtils::ResMinMass = {-99, -99, -99, -99, -99, -99} [static] |
Definition at line 119 of file MuScleFitUtils.h.
Referenced by probability(), and MuScleFitBase::readProbabilityDistributionsFromFile().
int MuScleFitUtils::ResolFitType = 0 [static] |
Definition at line 134 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), MuScleFit::MuScleFit(), and ResolutionAnalyzer::ResolutionAnalyzer().
resolutionFunctionBase< double * > * MuScleFitUtils::resolutionFunction = 0 [static] |
Definition at line 135 of file MuScleFitUtils.h.
Referenced by ErrorsAnalyzer::fillHistograms(), ErrorsPropagationAnalyzer::fillHistograms(), massResolution(), ErrorsPropagationAnalyzer::massResolution(), MuScleFit::MuScleFit(), ResolutionAnalyzer::ResolutionAnalyzer(), and TestCorrection::TestCorrection().
resolutionFunctionBase< std::vector< double > > * MuScleFitUtils::resolutionFunctionForVec = 0 [static] |
Definition at line 136 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 222 of file MuScleFitUtils.h.
Referenced by likelihood(), and minimizeLikelihood().
std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::SavedPair [static] |
Definition at line 213 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 218 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood().
int MuScleFitUtils::ScaleFitType = 0 [static] |
Definition at line 137 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
scaleFunctionBase< double * > * MuScleFitUtils::scaleFunction = 0 [static] |
Definition at line 138 of file MuScleFitUtils.h.
Referenced by applyScale(), minimizeLikelihood(), and MuScleFit::MuScleFit().
scaleFunctionBase< std::vector< double > > * MuScleFitUtils::scaleFunctionForVec = 0 [static] |
Definition at line 139 of file MuScleFitUtils.h.
Referenced by MuScleFit::MuScleFit().
bool MuScleFitUtils::separateRanges_ = true [static] |
Definition at line 238 of file MuScleFitUtils.h.
Referenced by MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
bool MuScleFitUtils::sherpa_ = false [static] |
Definition at line 228 of file MuScleFitUtils.h.
Referenced by findGenMuFromRes(), and MuScleFit::MuScleFit().
TH1D * MuScleFitUtils::signalProb_ = 0 [static] |
Definition at line 167 of file MuScleFitUtils.h.
Referenced by massProb(), and minimizeLikelihood().
std::vector< std::pair< lorentzVector, lorentzVector > > MuScleFitUtils::simPair [static] |
Definition at line 216 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), and MuScleFit::selectMuons().
smearFunctionBase * MuScleFitUtils::smearFunction = 0 [static] |
Definition at line 130 of file MuScleFitUtils.h.
Referenced by applySmearing(), and MuScleFit::MuScleFit().
int MuScleFitUtils::SmearType = 0 [static] |
Definition at line 129 of file MuScleFitUtils.h.
Referenced by MuScleFit::applySmearing(), applySmearing(), MuScleFit::checkParameters(), MuScleFit::MuScleFit(), and MuScleFit::selectMuons().
bool MuScleFitUtils::speedup = false [static] |
Definition at line 194 of file MuScleFitUtils.h.
Referenced by MuScleFit::duringFastLoop(), MuScleFit::MuScleFit(), MuScleFit::selectMuons(), and MuScleFit::~MuScleFit().
bool MuScleFitUtils::startWithSimplex_ [static] |
Definition at line 260 of file MuScleFitUtils.h.
Referenced by minimizeLikelihood(), and MuScleFit::MuScleFit().
const int MuScleFitUtils::totalResNum = 6 [static] |
Definition at line 115 of file MuScleFitUtils.h.
Referenced by massProb().
bool MuScleFitUtils::useProbsFile_ = true [static] |
Definition at line 235 of file MuScleFitUtils.h.
Referenced by MuScleFit::beginOfJobInConstructor(), findBestRecoRes(), and MuScleFit::MuScleFit().
double MuScleFitUtils::x [static] |
Definition at line 195 of file MuScleFitUtils.h.
Referenced by applySmearing(), fitMass(), fitReso(), massProb(), and MuScleFit::MuScleFit().