Reference Manual

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