Reference Manual

Go to the documentation of this file.
00001 #ifndef MUSCLEFITBASE_C
00002 #define MUSCLEFITBASE_C
00003 
00004 #include "MuonAnalysis/MomentumScaleCalibration/interface/MuScleFitBase.h"
00005 #include "FWCore/ParameterSet/interface/FileInPath.h"
00006 
00007 void MuScleFitBase::fillHistoMap(TFile* outputFile, unsigned int iLoop) {
00008   //Reconstructed muon kinematics
00009   //-----------------------------
00010   outputFile->cd();
00011   // If no Z is required, use a smaller mass range.
00012   double minMass = 0.;
00013   double maxMass = 200.;
00014   double maxPt = 100.;
00015   double yMaxEta = 4.;
00016   double yMaxPt = 2.;
00017   if( MuScleFitUtils::resfind[0] != 1 ) {
00018     maxMass = 20.;
00019     maxPt = 20.;
00020     yMaxEta = 0.2;
00021     yMaxPt = 0.2;
00022     // If running on standalone muons we need to expand the window range
00023     if( theMuonType_ == 2 ) {
00024       yMaxEta = 20.;
00025     }
00026   }
00027 
00028   LogDebug("MuScleFitBase") << "Creating new histograms" << std::endl;
00029 
00030   mapHisto_["hRecBestMu"]      = new HParticle ("hRecBestMu", minMass, maxMass, maxPt);
00031   mapHisto_["hRecBestMuVSEta"] = new HPartVSEta ("hRecBestMuVSEta", minMass, maxMass, maxPt);
00032   //mapHisto_["hRecBestMu_Acc"]  = new HParticle ("hRecBestMu_Acc", minMass, maxMass);
00033   mapHisto_["hDeltaRecBestMu"] = new HDelta ("hDeltaRecBestMu");
00034 
00035   mapHisto_["hRecBestRes"]          = new HParticle   ("hRecBestRes", minMass, maxMass, maxPt);
00036   mapHisto_["hRecBestResAllEvents"] = new HParticle   ("hRecBestResAllEvents", minMass, maxMass, maxPt);
00037   //mapHisto_["hRecBestRes_Acc"] = new HParticle   ("hRecBestRes_Acc", minMass, maxMass);
00038   // If not finding Z, use a smaller mass window
00039   mapHisto_["hRecBestResVSMu"] = new HMassVSPart ("hRecBestResVSMu", minMass, maxMass, maxPt);
00040   mapHisto_["hRecBestResVSRes"] = new HMassVSPart ("hRecBestResVSRes", minMass, maxMass, maxPt);
00041   //Generated Mass versus pt
00042   mapHisto_["hGenResVSMu"] = new HMassVSPart ("hGenResVSMu", minMass, maxMass, maxPt);
00043   // Likelihood values VS muon variables
00044   // -------------------------------------
00045   mapHisto_["hLikeVSMu"]       = new HLikelihoodVSPart ("hLikeVSMu");
00046   mapHisto_["hLikeVSMuMinus"]  = new HLikelihoodVSPart ("hLikeVSMuMinus");
00047   mapHisto_["hLikeVSMuPlus"]   = new HLikelihoodVSPart ("hLikeVSMuPlus");
00048 
00049   //Resolution VS muon kinematic
00050   //----------------------------
00051   mapHisto_["hResolMassVSMu"]         = new HResolutionVSPart( outputFile, "hResolMassVSMu", maxPt, 0., yMaxEta, 0., yMaxPt, true );
00052   mapHisto_["hFunctionResolMassVSMu"] = new HResolutionVSPart( outputFile, "hFunctionResolMassVSMu", maxPt, 0, 0.1, 0, 0.1, true );
00053   mapHisto_["hResolPtGenVSMu"]        = new HResolutionVSPart( outputFile, "hResolPtGenVSMu", maxPt, -0.1, 0.1, -0.1, 0.1 );
00054   mapHisto_["hResolPtSimVSMu"]        = new HResolutionVSPart( outputFile, "hResolPtSimVSMu", maxPt, -0.1, 0.1, -0.1, 0.1 );
00055   mapHisto_["hResolEtaGenVSMu"]       = new HResolutionVSPart( outputFile, "hResolEtaGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 );
00056   mapHisto_["hResolEtaSimVSMu"]       = new HResolutionVSPart( outputFile, "hResolEtaSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 );
00057   mapHisto_["hResolThetaGenVSMu"]     = new HResolutionVSPart( outputFile, "hResolThetaGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 );
00058   mapHisto_["hResolThetaSimVSMu"]     = new HResolutionVSPart( outputFile, "hResolThetaSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 );
00059   mapHisto_["hResolCotgThetaGenVSMu"] = new HResolutionVSPart( outputFile, "hResolCotgThetaGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 );
00060   mapHisto_["hResolCotgThetaSimVSMu"] = new HResolutionVSPart( outputFile, "hResolCotgThetaSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 );
00061   mapHisto_["hResolPhiGenVSMu"]       = new HResolutionVSPart( outputFile, "hResolPhiGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 );
00062   mapHisto_["hResolPhiSimVSMu"]       = new HResolutionVSPart( outputFile, "hResolPhiSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02 );
00063 
00064   if( MuScleFitUtils::debugMassResol_ ) {
00065     mapHisto_["hdMdPt1"] = new HResolutionVSPart( outputFile, "hdMdPt1", maxPt, 0, 100, -3.2, 3.2, true );
00066     mapHisto_["hdMdPt2"] = new HResolutionVSPart( outputFile, "hdMdPt2", maxPt, 0, 100, -3.2, 3.2, true );
00067     mapHisto_["hdMdPhi1"] = new HResolutionVSPart( outputFile, "hdMdPhi1", maxPt, 0, 100, -3.2, 3.2, true );
00068     mapHisto_["hdMdPhi2"] = new HResolutionVSPart( outputFile, "hdMdPhi2", maxPt, 0, 100, -3.2, 3.2, true );
00069     mapHisto_["hdMdCotgTh1"] = new HResolutionVSPart( outputFile, "hdMdCotgTh1", maxPt, 0, 100, -3.2, 3.2, true );
00070     mapHisto_["hdMdCotgTh2"] = new HResolutionVSPart( outputFile, "hdMdCotgTh2", maxPt, 0, 100, -3.2, 3.2, true );
00071   }
00072 
00073   HTH2D * recoGenHisto = new HTH2D(outputFile, "hPtRecoVsPtGen", "Pt reco vs Pt gen", "hPtRecoVsPtGen", 120, 0., 120., 120, 0, 120.);
00074   (*recoGenHisto)->SetXTitle("Pt gen (GeV)");
00075   (*recoGenHisto)->SetYTitle("Pt reco (GeV)");
00076   mapHisto_["hPtRecoVsPtGen"] = recoGenHisto;
00077   HTH2D * recoSimHisto = new HTH2D(outputFile, "hPtRecoVsPtSim", "Pt reco vs Pt sim", "hPtRecoVsPtSim", 120, 0., 120., 120, 0, 120.);
00078   (*recoSimHisto)->SetXTitle("Pt sim (GeV)");
00079   (*recoSimHisto)->SetYTitle("Pt reco (GeV)");
00080   mapHisto_["hPtRecoVsPtSim"] = recoSimHisto;
00081   // Resolutions from resolution functions
00082   // -------------------------------------
00083   mapHisto_["hFunctionResolPt"]        = new HFunctionResolution( outputFile, "hFunctionResolPt", maxPt );
00084   mapHisto_["hFunctionResolCotgTheta"] = new HFunctionResolution( outputFile, "hFunctionResolCotgTheta", maxPt );
00085   mapHisto_["hFunctionResolPhi"]       = new HFunctionResolution( outputFile, "hFunctionResolPhi", maxPt );
00086 
00087   // Mass probability histograms
00088   // ---------------------------
00089   // The word "profile" is added to the title automatically
00090   mapHisto_["hMass_P"]      = new HTProfile( outputFile, "Mass_P", "Mass probability", 4000, 0., 200., 0., 50. );
00091   mapHisto_["hMass_fine_P"] = new HTProfile( outputFile, "Mass_fine_P", "Mass probability", 4000, 0., 20., 0., 50. );
00092   mapHisto_["hMass_Probability"]      = new HTH1D( outputFile, "Mass_Probability", "Mass probability", 4000, 0., 200.);
00093   mapHisto_["hMass_fine_Probability"] = new HTH1D( outputFile, "Mass_fine_Probability", "Mass probability", 4000, 0., 20.);
00094   mapHisto_["hMassProbVsMu"] = new HMassVSPartProfile( "hMassProbVsMu", minMass, maxMass, maxPt );
00095   mapHisto_["hMassProbVsRes"] = new HMassVSPartProfile( "hMassProbVsRes", minMass, maxMass, maxPt );
00096   mapHisto_["hMassProbVsMu_fine"] = new HMassVSPartProfile( "hMassProbVsMu_fine", minMass, maxMass, maxPt );
00097   mapHisto_["hMassProbVsRes_fine"] = new HMassVSPartProfile( "hMassProbVsRes_fine", minMass, maxMass, maxPt );
00098 
00099   // (M_reco-M_gen)/M_gen vs (pt, eta) of the muons from MC
00100   mapHisto_["hDeltaMassOverGenMassVsPt"] = new HTH2D( outputFile, "DeltaMassOverGenMassVsPt", "DeltaMassOverGenMassVsPt", "DeltaMassOverGenMass", 200, 0, maxPt, 200, -0.2, 0.2 );
00101   mapHisto_["hDeltaMassOverGenMassVsEta"] = new HTH2D( outputFile, "DeltaMassOverGenMassVsEta", "DeltaMassOverGenMassVsEta", "DeltaMassOverGenMass", 200, -3., 3., 200, -0.2, 0.2 );
00102 
00103   // Square of mass resolution vs (pt, eta) of the muons from MC
00104   mapHisto_["hMassResolutionVsPtEta"] = new HCovarianceVSxy( "Mass", "Mass", 100, 0., maxPt, 60, -3, 3, outputFile->mkdir("MassCovariance") );
00105   // Mass resolution vs (pt, eta) from resolution function
00106   mapHisto_["hFunctionResolMass"] = new HFunctionResolution( outputFile, "hFunctionResolMass", maxPt );
00107 }
00108 
00109 void MuScleFitBase::clearHistoMap() {
00110   for (std::map<std::string, Histograms*>::const_iterator histo=mapHisto_.begin();
00111        histo!=mapHisto_.end(); histo++) {
00112     delete (*histo).second;
00113   }
00114 }
00115 
00116 void MuScleFitBase::writeHistoMap( const unsigned int iLoop ) {
00117   for (std::map<std::string, Histograms*>::const_iterator histo=mapHisto_.begin();
00118        histo!=mapHisto_.end(); histo++) {
00119     // This is to avoid writing into subdirs. Need a workaround.
00120     theFiles_[iLoop]->cd();
00121     (*histo).second->Write();
00122   }
00123 }
00124 
00125 void MuScleFitBase::readProbabilityDistributionsFromFile()
00126 {
00127   TH2D * GLZ[24];
00128   TH2D * GL[6];
00129   TFile * ProbsFile;
00130   if( probabilitiesFile_ != "" ) {
00131     ProbsFile = new TFile (probabilitiesFile_.c_str());
00132     std::cout << "[MuScleFit-Constructor]: Reading TH2D probabilities from " << probabilitiesFile_ << std::endl;
00133   }
00134   else {
00135     // edm::FileInPath file("MuonAnalysis/MomentumScaleCalibration/test/Probs_new_1000_CTEQ.root");
00136     // edm::FileInPath file("MuonAnalysis/MomentumScaleCalibration/test/Probs_new_Horace_CTEQ_1000.root");
00137     // edm::FileInPath file("MuonAnalysis/MomentumScaleCalibration/test/Probs_merge.root");
00138     edm::FileInPath file(probabilitiesFileInPath_.c_str());
00139     ProbsFile = new TFile (file.fullPath().c_str());
00140     std::cout << "[MuScleFit-Constructor]: Reading TH2D probabilities from " << probabilitiesFileInPath_ << std::endl;
00141   }
00142 
00143   ProbsFile->cd();
00144   if( theMuonType_!=2 && MuScleFitUtils::resfind[0]) {
00145     for ( int i=0; i<24; i++ ) {
00146       char nameh[6];
00147       sprintf (nameh,"GLZ%d",i);
00148       GLZ[i] = dynamic_cast<TH2D*>(ProbsFile->Get(nameh));
00149     }
00150   }
00151   if( theMuonType_==2 && MuScleFitUtils::resfind[0]) 
00152     GL[0] = dynamic_cast<TH2D*> (ProbsFile->Get("GL0"));
00153   if(MuScleFitUtils::resfind[1]) 
00154     GL[1] = dynamic_cast<TH2D*> (ProbsFile->Get("GL1"));
00155   if(MuScleFitUtils::resfind[2]) 
00156     GL[2] = dynamic_cast<TH2D*> (ProbsFile->Get("GL2"));
00157   if(MuScleFitUtils::resfind[3]) 
00158     GL[3] = dynamic_cast<TH2D*> (ProbsFile->Get("GL3"));
00159   if(MuScleFitUtils::resfind[4]) 
00160     GL[4] = dynamic_cast<TH2D*> (ProbsFile->Get("GL4"));
00161   if(MuScleFitUtils::resfind[5]) 
00162     GL[5] = dynamic_cast<TH2D*> (ProbsFile->Get("GL5"));
00163 
00164   // Read the limits for M and Sigma axis for each pdf
00165   // Note: we assume all the Z histograms to have the same limits
00166   // x is mass, y is sigma
00167   if(MuScleFitUtils::resfind[0] && theMuonType_!=2) {
00168     MuScleFitUtils::ResHalfWidth[0] = (GLZ[0]->GetXaxis()->GetXmax() - GLZ[0]->GetXaxis()->GetXmin())/2.;
00169     MuScleFitUtils::ResMaxSigma[0] = (GLZ[0]->GetYaxis()->GetXmax() - GLZ[0]->GetYaxis()->GetXmin());
00170     MuScleFitUtils::ResMinMass[0] = (GLZ[0]->GetXaxis()->GetXmin());
00171   }
00172   if(MuScleFitUtils::resfind[0] && theMuonType_==2) {
00173     MuScleFitUtils::ResHalfWidth[0] = (GL[0]->GetXaxis()->GetXmax() - GL[0]->GetXaxis()->GetXmin())/2.;
00174     MuScleFitUtils::ResMaxSigma[0] = (GL[0]->GetYaxis()->GetXmax() - GL[0]->GetYaxis()->GetXmin());
00175     MuScleFitUtils::ResMinMass[0] = (GL[0]->GetXaxis()->GetXmin());
00176   }
00177   for( int i=1; i<6; ++i ) {
00178     if(MuScleFitUtils::resfind[i]){
00179       MuScleFitUtils::ResHalfWidth[i] = (GL[i]->GetXaxis()->GetXmax() - GL[i]->GetXaxis()->GetXmin())/2.;
00180       MuScleFitUtils::ResMaxSigma[i] = (GL[i]->GetYaxis()->GetXmax() - GL[i]->GetYaxis()->GetXmin());
00181       MuScleFitUtils::ResMinMass[i] = (GL[i]->GetXaxis()->GetXmin());
00182      // if( debug_>2 ) {
00183       std::cout << "MuScleFitUtils::ResHalfWidth["<<i<<"] = " << MuScleFitUtils::ResHalfWidth[i] << std::endl;
00184       std::cout << "MuScleFitUtils::ResMaxSigma["<<i<<"] = " << MuScleFitUtils::ResMaxSigma[i] << std::endl;
00185       // }
00186     }
00187   }
00188 
00189   // Extract normalization for mass slice in Y bins of Z
00190   // ---------------------------------------------------
00191   if(MuScleFitUtils::resfind[0] && theMuonType_!=2) {
00192     for (int iY=0; iY<24; iY++) {
00193       int nBinsX = GLZ[iY]->GetNbinsX();
00194       int nBinsY = GLZ[iY]->GetNbinsY();
00195       if( nBinsX != MuScleFitUtils::nbins+1 || nBinsY != MuScleFitUtils::nbins+1 ) {
00196         std::cout << "Error: for histogram \"" << GLZ[iY]->GetName() << "\" bins are not " << MuScleFitUtils::nbins << std::endl;
00197         std::cout<< "nBinsX = " << nBinsX << ", nBinsY = " << nBinsY << std::endl;
00198         exit(1);
00199       }
00200       for (int iy=0; iy<=MuScleFitUtils::nbins; iy++) {
00201         MuScleFitUtils::GLZNorm[iY][iy] = 0.;
00202         for (int ix=0; ix<=MuScleFitUtils::nbins; ix++) {
00203           MuScleFitUtils::GLZValue[iY][ix][iy] = GLZ[iY]->GetBinContent(ix+1, iy+1);
00204           MuScleFitUtils::GLZNorm[iY][iy] += MuScleFitUtils::GLZValue[iY][ix][iy]*(2*MuScleFitUtils::ResHalfWidth[0])/MuScleFitUtils::nbins;
00205         }
00206       }
00207     }
00208   }
00209 
00210   if(MuScleFitUtils::resfind[0] && theMuonType_==2){
00211       int nBinsX = GL[0]->GetNbinsX();
00212       int nBinsY = GL[0]->GetNbinsY();
00213       if( nBinsX != MuScleFitUtils::nbins+1 || nBinsY != MuScleFitUtils::nbins+1 ) {
00214         std::cout << "Error: for histogram \"" << GL[0]->GetName() << "\" bins are not " << MuScleFitUtils::nbins << std::endl;
00215         std::cout<< "nBinsX = " << nBinsX << ", nBinsY = " << nBinsY << std::endl;
00216         exit(1);
00217       }
00218 
00219       for (int iy=0; iy<=MuScleFitUtils::nbins; iy++) {
00220         MuScleFitUtils::GLNorm[0][iy] = 0.;
00221         for (int ix=0; ix<=MuScleFitUtils::nbins; ix++) {
00222           MuScleFitUtils::GLValue[0][ix][iy] = GL[0]->GetBinContent(ix+1, iy+1);
00223           // N.B. approximation: we should compute the integral of the function used to compute the probability (linear
00224           // interpolation of the mass points). This computation could be troublesome because the points have a steep
00225           // variation near the mass peak and the normal integral is not precise in these conditions.
00226           // Furthermore it is slow.
00227           MuScleFitUtils::GLNorm[0][iy] += MuScleFitUtils::GLValue[0][ix][iy]*(2*MuScleFitUtils::ResHalfWidth[0])/MuScleFitUtils::nbins;
00228         }
00229       }
00230     }  
00231   // Extract normalization for each mass slice
00232   // -----------------------------------------
00233   for (int ires=1; ires<6; ires++) {
00234     if(MuScleFitUtils::resfind[ires]){
00235       int nBinsX = GL[ires]->GetNbinsX();
00236       int nBinsY = GL[ires]->GetNbinsY();
00237       if( nBinsX != MuScleFitUtils::nbins+1 || nBinsY != MuScleFitUtils::nbins+1 ) {
00238         std::cout << "Error: for histogram \"" << GL[ires]->GetName() << "\" bins are not " << MuScleFitUtils::nbins << std::endl;
00239         std::cout<< "nBinsX = " << nBinsX << ", nBinsY = " << nBinsY << std::endl;
00240         exit(1);
00241       }
00242 
00243       for (int iy=0; iy<=MuScleFitUtils::nbins; iy++) {
00244         MuScleFitUtils::GLNorm[ires][iy] = 0.;
00245         for (int ix=0; ix<=MuScleFitUtils::nbins; ix++) {
00246           MuScleFitUtils::GLValue[ires][ix][iy] = GL[ires]->GetBinContent(ix+1, iy+1);
00247           // N.B. approximation: we should compute the integral of the function used to compute the probability (linear
00248           // interpolation of the mass points). This computation could be troublesome because the points have a steep
00249           // variation near the mass peak and the normal integral is not precise in these conditions.
00250           // Furthermore it is slow.
00251           MuScleFitUtils::GLNorm[ires][iy] += MuScleFitUtils::GLValue[ires][ix][iy]*(2*MuScleFitUtils::ResHalfWidth[ires])/MuScleFitUtils::nbins;
00252         }
00253       }
00254     }
00255   }
00256   // Free all the memory for the probability histograms.
00257   if(MuScleFitUtils::resfind[0] && theMuonType_!=2) {
00258     for ( int i=0; i<24; i++ ) {
00259       delete GLZ[i];
00260     }
00261   }
00262   if(MuScleFitUtils::resfind[0] && theMuonType_==2)
00263     delete GL[0];
00264   for (int ires=1; ires<6; ires++) {
00265     if(MuScleFitUtils::resfind[ires])
00266       delete GL[ires];
00267   }
00268   delete ProbsFile;
00269 }
00270 
00271 #endif
/data/refman/pasoursint/CMSSW_6_1_1/src/MuonAnalysis/MomentumScaleCalibration/src/MuScleFitBase.cc
