Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:22:39 

0001 #ifndef MUSCLEFITBASE_C
0002 #define MUSCLEFITBASE_C
0003 
0004 #include "MuScleFitBase.h"
0005 #include "FWCore/ParameterSet/interface/FileInPath.h"
0006 #include <memory>
0007 #include <array>
0008 
0009 void MuScleFitBase::fillHistoMap(TFile* outputFile, unsigned int iLoop) {
0010   //Reconstructed muon kinematics
0011   //-----------------------------
0012   outputFile->cd();
0013   // If no Z is required, use a smaller mass range.
0014   double minMass = 0.;
0015   double maxMass = 200.;
0016   double maxPt = 100.;
0017   double yMaxEta = 4.;
0018   double yMaxPt = 2.;
0019   if (MuScleFitUtils::resfind[0] != 1) {
0020     maxMass = 20.;
0021     maxPt = 20.;
0022     yMaxEta = 0.2;
0023     yMaxPt = 0.2;
0024     // If running on standalone muons we need to expand the window range
0025     if (theMuonType_ == 2) {
0026       yMaxEta = 20.;
0027     }
0028   }
0029 
0030   LogDebug("MuScleFitBase") << "Creating new histograms" << std::endl;
0031 
0032   mapHisto_["hRecBestMu"] = new HParticle("hRecBestMu", minMass, maxMass, maxPt);
0033   mapHisto_["hRecBestMuVSEta"] = new HPartVSEta("hRecBestMuVSEta", minMass, maxMass, maxPt);
0034   //mapHisto_["hRecBestMuVSPhi"] = new HPartVSPhi ("hRecBestMuVSPhi");
0035   //mapHisto_["hRecBestMu_Acc"]  = new HParticle ("hRecBestMu_Acc", minMass, maxMass);
0036   mapHisto_["hDeltaRecBestMu"] = new HDelta("hDeltaRecBestMu");
0037 
0038   mapHisto_["hRecBestRes"] = new HParticle("hRecBestRes", minMass, maxMass, maxPt);
0039   mapHisto_["hRecBestResAllEvents"] = new HParticle("hRecBestResAllEvents", minMass, maxMass, maxPt);
0040   //mapHisto_["hRecBestRes_Acc"] = new HParticle   ("hRecBestRes_Acc", minMass, maxMass);
0041   // If not finding Z, use a smaller mass window
0042   mapHisto_["hRecBestResVSMu"] = new HMassVSPart("hRecBestResVSMu", minMass, maxMass, maxPt);
0043   mapHisto_["hRecBestResVSRes"] = new HMassVSPart("hRecBestResVSRes", minMass, maxMass, maxPt);
0044   //Generated Mass versus pt
0045   mapHisto_["hGenResVSMu"] = new HMassVSPart("hGenResVSMu", minMass, maxMass, maxPt);
0046   // Likelihood values VS muon variables
0047   // -------------------------------------
0048   mapHisto_["hLikeVSMu"] = new HLikelihoodVSPart("hLikeVSMu");
0049   mapHisto_["hLikeVSMuMinus"] = new HLikelihoodVSPart("hLikeVSMuMinus");
0050   mapHisto_["hLikeVSMuPlus"] = new HLikelihoodVSPart("hLikeVSMuPlus");
0051 
0052   //Resolution VS muon kinematic
0053   //----------------------------
0054   mapHisto_["hResolMassVSMu"] =
0055       new HResolutionVSPart(outputFile, "hResolMassVSMu", maxPt, 0., yMaxEta, 0., yMaxPt, true);
0056   mapHisto_["hFunctionResolMassVSMu"] =
0057       new HResolutionVSPart(outputFile, "hFunctionResolMassVSMu", maxPt, 0, 0.1, 0, 0.1, true);
0058   mapHisto_["hResolPtGenVSMu"] = new HResolutionVSPart(outputFile, "hResolPtGenVSMu", maxPt, -0.1, 0.1, -0.1, 0.1);
0059   mapHisto_["hResolPtSimVSMu"] = new HResolutionVSPart(outputFile, "hResolPtSimVSMu", maxPt, -0.1, 0.1, -0.1, 0.1);
0060   mapHisto_["hResolEtaGenVSMu"] =
0061       new HResolutionVSPart(outputFile, "hResolEtaGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02);
0062   mapHisto_["hResolEtaSimVSMu"] =
0063       new HResolutionVSPart(outputFile, "hResolEtaSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02);
0064   mapHisto_["hResolThetaGenVSMu"] =
0065       new HResolutionVSPart(outputFile, "hResolThetaGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02);
0066   mapHisto_["hResolThetaSimVSMu"] =
0067       new HResolutionVSPart(outputFile, "hResolThetaSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02);
0068   mapHisto_["hResolCotgThetaGenVSMu"] =
0069       new HResolutionVSPart(outputFile, "hResolCotgThetaGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02);
0070   mapHisto_["hResolCotgThetaSimVSMu"] =
0071       new HResolutionVSPart(outputFile, "hResolCotgThetaSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02);
0072   mapHisto_["hResolPhiGenVSMu"] =
0073       new HResolutionVSPart(outputFile, "hResolPhiGenVSMu", maxPt, -0.02, 0.02, -0.02, 0.02);
0074   mapHisto_["hResolPhiSimVSMu"] =
0075       new HResolutionVSPart(outputFile, "hResolPhiSimVSMu", maxPt, -0.02, 0.02, -0.02, 0.02);
0076 
0077   if (MuScleFitUtils::debugMassResol_) {
0078     mapHisto_["hdMdPt1"] = new HResolutionVSPart(outputFile, "hdMdPt1", maxPt, 0, 100, -3.2, 3.2, true);
0079     mapHisto_["hdMdPt2"] = new HResolutionVSPart(outputFile, "hdMdPt2", maxPt, 0, 100, -3.2, 3.2, true);
0080     mapHisto_["hdMdPhi1"] = new HResolutionVSPart(outputFile, "hdMdPhi1", maxPt, 0, 100, -3.2, 3.2, true);
0081     mapHisto_["hdMdPhi2"] = new HResolutionVSPart(outputFile, "hdMdPhi2", maxPt, 0, 100, -3.2, 3.2, true);
0082     mapHisto_["hdMdCotgTh1"] = new HResolutionVSPart(outputFile, "hdMdCotgTh1", maxPt, 0, 100, -3.2, 3.2, true);
0083     mapHisto_["hdMdCotgTh2"] = new HResolutionVSPart(outputFile, "hdMdCotgTh2", maxPt, 0, 100, -3.2, 3.2, true);
0084   }
0085 
0086   HTH2D* recoGenHisto =
0087       new HTH2D(outputFile, "hPtRecoVsPtGen", "Pt reco vs Pt gen", "hPtRecoVsPtGen", 120, 0., 120., 120, 0, 120.);
0088   (*recoGenHisto)->SetXTitle("Pt gen (GeV)");
0089   (*recoGenHisto)->SetYTitle("Pt reco (GeV)");
0090   mapHisto_["hPtRecoVsPtGen"] = recoGenHisto;
0091   HTH2D* recoSimHisto =
0092       new HTH2D(outputFile, "hPtRecoVsPtSim", "Pt reco vs Pt sim", "hPtRecoVsPtSim", 120, 0., 120., 120, 0, 120.);
0093   (*recoSimHisto)->SetXTitle("Pt sim (GeV)");
0094   (*recoSimHisto)->SetYTitle("Pt reco (GeV)");
0095   mapHisto_["hPtRecoVsPtSim"] = recoSimHisto;
0096   // Resolutions from resolution functions
0097   // -------------------------------------
0098   mapHisto_["hFunctionResolPt"] = new HFunctionResolution(outputFile, "hFunctionResolPt", maxPt);
0099   mapHisto_["hFunctionResolCotgTheta"] = new HFunctionResolution(outputFile, "hFunctionResolCotgTheta", maxPt);
0100   mapHisto_["hFunctionResolPhi"] = new HFunctionResolution(outputFile, "hFunctionResolPhi", maxPt);
0101 
0102   // Mass probability histograms
0103   // ---------------------------
0104   // The word "profile" is added to the title automatically
0105   mapHisto_["hMass_P"] = new HTProfile(outputFile, "Mass_P", "Mass probability", 4000, 0., 200., 0., 50.);
0106   mapHisto_["hMass_fine_P"] = new HTProfile(outputFile, "Mass_fine_P", "Mass probability", 4000, 0., 20., 0., 50.);
0107   mapHisto_["hMass_Probability"] = new HTH1D(outputFile, "Mass_Probability", "Mass probability", 4000, 0., 200.);
0108   mapHisto_["hMass_fine_Probability"] =
0109       new HTH1D(outputFile, "Mass_fine_Probability", "Mass probability", 4000, 0., 20.);
0110   mapHisto_["hMassProbVsMu"] = new HMassVSPartProfile("hMassProbVsMu", minMass, maxMass, maxPt);
0111   mapHisto_["hMassProbVsRes"] = new HMassVSPartProfile("hMassProbVsRes", minMass, maxMass, maxPt);
0112   mapHisto_["hMassProbVsMu_fine"] = new HMassVSPartProfile("hMassProbVsMu_fine", minMass, maxMass, maxPt);
0113   mapHisto_["hMassProbVsRes_fine"] = new HMassVSPartProfile("hMassProbVsRes_fine", minMass, maxMass, maxPt);
0114 
0115   // (M_reco-M_gen)/M_gen vs (pt, eta) of the muons from MC
0116   mapHisto_["hDeltaMassOverGenMassVsPt"] = new HTH2D(outputFile,
0117                                                      "DeltaMassOverGenMassVsPt",
0118                                                      "DeltaMassOverGenMassVsPt",
0119                                                      "DeltaMassOverGenMass",
0120                                                      200,
0121                                                      0,
0122                                                      maxPt,
0123                                                      200,
0124                                                      -0.2,
0125                                                      0.2);
0126   mapHisto_["hDeltaMassOverGenMassVsEta"] = new HTH2D(outputFile,
0127                                                       "DeltaMassOverGenMassVsEta",
0128                                                       "DeltaMassOverGenMassVsEta",
0129                                                       "DeltaMassOverGenMass",
0130                                                       200,
0131                                                       -3.,
0132                                                       3.,
0133                                                       200,
0134                                                       -0.2,
0135                                                       0.2);
0136 
0137   // Square of mass resolution vs (pt, eta) of the muons from MC
0138   // EM 2012.12.19  mapHisto_["hMassResolutionVsPtEta"] = new HCovarianceVSxy( "Mass", "Mass", 100, 0., maxPt, 60, -3, 3, outputFile->mkdir("MassCovariance") );
0139   // Mass resolution vs (pt, eta) from resolution function
0140   mapHisto_["hFunctionResolMass"] = new HFunctionResolution(outputFile, "hFunctionResolMass", maxPt);
0141 }
0142 
0143 void MuScleFitBase::clearHistoMap() {
0144   for (std::map<std::string, Histograms*>::const_iterator histo = mapHisto_.begin(); histo != mapHisto_.end();
0145        histo++) {
0146     delete (*histo).second;
0147   }
0148 }
0149 
0150 void MuScleFitBase::writeHistoMap(const unsigned int iLoop) {
0151   for (std::map<std::string, Histograms*>::const_iterator histo = mapHisto_.begin(); histo != mapHisto_.end();
0152        histo++) {
0153     // This is to avoid writing into subdirs. Need a workaround.
0154     theFiles_[iLoop]->cd();
0155     (*histo).second->Write();
0156   }
0157 }
0158 
0159 void MuScleFitBase::readProbabilityDistributionsFromFile() {
0160   std::array<TH2D*, 6> GLZ = {{nullptr}};
0161   std::array<TH2D*, 6> GL = {{nullptr}};
0162   std::unique_ptr<TFile> ProbsFile;
0163   if (!probabilitiesFile_.empty()) {
0164     ProbsFile = std::make_unique<TFile>(probabilitiesFile_.c_str());
0165     std::cout << "[MuScleFit-Constructor]: Reading TH2D probabilities from " << probabilitiesFile_ << std::endl;
0166   } else {
0167     // edm::FileInPath file("MuonAnalysis/MomentumScaleCalibration/test/Probs_new_1000_CTEQ.root");
0168     // edm::FileInPath file("MuonAnalysis/MomentumScaleCalibration/test/Probs_new_Horace_CTEQ_1000.root");
0169     // edm::FileInPath file("MuonAnalysis/MomentumScaleCalibration/test/Probs_merge.root");
0170     edm::FileInPath file(probabilitiesFileInPath_.c_str());
0171     ProbsFile = std::make_unique<TFile>(file.fullPath().c_str());
0172     std::cout << "[MuScleFit-Constructor]: Reading TH2D probabilities from " << probabilitiesFileInPath_ << std::endl;
0173   }
0174 
0175   ProbsFile->cd();
0176   bool resfindEmpty = true;
0177   if (MuScleFitUtils::rapidityBinsForZ_ && MuScleFitUtils::resfind[0] && theMuonType_ != 2) {
0178     resfindEmpty = false;
0179     for (unsigned char i = 0; i < GLZ.size(); i++) {
0180       char nameh[6];
0181       snprintf(nameh, 6, "GLZ%hhu", i);
0182       GLZ[i] = dynamic_cast<TH2D*>(ProbsFile->Get(nameh));
0183     }
0184   } else if (MuScleFitUtils::resfind[0] && theMuonType_ == 2) {
0185     GL[0] = dynamic_cast<TH2D*>(ProbsFile->Get("GL0"));
0186     resfindEmpty = false;
0187   }
0188   for (unsigned char i = 1; i < GL.size(); ++i) {
0189     if (MuScleFitUtils::resfind[i]) {
0190       char nameh[6];
0191       snprintf(nameh, 6, "GL%hhu", i);
0192       GL[i] = dynamic_cast<TH2D*>(ProbsFile->Get(nameh));
0193       resfindEmpty = false;
0194     }
0195   }
0196   if (resfindEmpty) {
0197     std::cout << "[MuScleFit-Constructor]: No resonance selected, please fill the resfind array" << std::endl;
0198     exit(1);
0199   }
0200 
0201   // Read the limits for M and Sigma axis for each pdf
0202   // Note: we assume all the Z histograms to have the same limits
0203   // x is mass, y is sigma
0204   if (MuScleFitUtils::rapidityBinsForZ_ && MuScleFitUtils::resfind[0] && theMuonType_ != 2) {
0205     MuScleFitUtils::ResHalfWidth[0] = (GLZ[0]->GetXaxis()->GetXmax() - GLZ[0]->GetXaxis()->GetXmin()) / 2.;
0206     MuScleFitUtils::ResMaxSigma[0] = (GLZ[0]->GetYaxis()->GetXmax() - GLZ[0]->GetYaxis()->GetXmin());
0207     MuScleFitUtils::ResMinMass[0] = (GLZ[0]->GetXaxis()->GetXmin());
0208   }
0209   if (MuScleFitUtils::resfind[0] && theMuonType_ == 2) {
0210     MuScleFitUtils::ResHalfWidth[0] = (GL[0]->GetXaxis()->GetXmax() - GL[0]->GetXaxis()->GetXmin()) / 2.;
0211     MuScleFitUtils::ResMaxSigma[0] = (GL[0]->GetYaxis()->GetXmax() - GL[0]->GetYaxis()->GetXmin());
0212     MuScleFitUtils::ResMinMass[0] = (GL[0]->GetXaxis()->GetXmin());
0213   }
0214   for (unsigned int i = 1; i < GL.size(); ++i) {
0215     if (MuScleFitUtils::resfind[i]) {
0216       MuScleFitUtils::ResHalfWidth[i] = (GL[i]->GetXaxis()->GetXmax() - GL[i]->GetXaxis()->GetXmin()) / 2.;
0217       MuScleFitUtils::ResMaxSigma[i] = (GL[i]->GetYaxis()->GetXmax() - GL[i]->GetYaxis()->GetXmin());
0218       MuScleFitUtils::ResMinMass[i] = (GL[i]->GetXaxis()->GetXmin());
0219       // if( debug_>2 ) {
0220       std::cout << "MuScleFitUtils::ResHalfWidth[" << i << "] = " << MuScleFitUtils::ResHalfWidth[i] << std::endl;
0221       std::cout << "MuScleFitUtils::ResMaxSigma[" << i << "] = " << MuScleFitUtils::ResMaxSigma[i] << std::endl;
0222       // }
0223     }
0224   }
0225 
0226   // Extract normalization for mass slice in Y bins of Z
0227   // ---------------------------------------------------
0228   if (MuScleFitUtils::rapidityBinsForZ_ && MuScleFitUtils::resfind[0] && theMuonType_ != 2) {
0229     for (unsigned int iY = 0; iY < GLZ.size(); iY++) {
0230       int nBinsX = GLZ[iY]->GetNbinsX();
0231       int nBinsY = GLZ[iY]->GetNbinsY();
0232       if (nBinsX != MuScleFitUtils::nbins + 1 || nBinsY != MuScleFitUtils::nbins + 1) {
0233         std::cout << "Error: for histogram \"" << GLZ[iY]->GetName() << "\" bins are not " << MuScleFitUtils::nbins
0234                   << std::endl;
0235         std::cout << "nBinsX = " << nBinsX << ", nBinsY = " << nBinsY << std::endl;
0236         exit(1);
0237       }
0238       for (int iy = 0; iy <= MuScleFitUtils::nbins; iy++) {
0239         MuScleFitUtils::GLZNorm[iY][iy] = 0.;
0240         for (int ix = 0; ix <= MuScleFitUtils::nbins; ix++) {
0241           MuScleFitUtils::GLZValue[iY][ix][iy] = GLZ[iY]->GetBinContent(ix + 1, iy + 1);
0242           MuScleFitUtils::GLZNorm[iY][iy] +=
0243               MuScleFitUtils::GLZValue[iY][ix][iy] * (2 * MuScleFitUtils::ResHalfWidth[0]) / MuScleFitUtils::nbins;
0244         }
0245       }
0246     }
0247   }
0248 
0249   if (MuScleFitUtils::resfind[0] && theMuonType_ == 2) {
0250     int nBinsX = GL[0]->GetNbinsX();
0251     int nBinsY = GL[0]->GetNbinsY();
0252     if (nBinsX != MuScleFitUtils::nbins + 1 || nBinsY != MuScleFitUtils::nbins + 1) {
0253       std::cout << "Error: for histogram \"" << GL[0]->GetName() << "\" bins are not " << MuScleFitUtils::nbins
0254                 << std::endl;
0255       std::cout << "nBinsX = " << nBinsX << ", nBinsY = " << nBinsY << std::endl;
0256       exit(1);
0257     }
0258 
0259     for (int iy = 0; iy <= MuScleFitUtils::nbins; iy++) {
0260       MuScleFitUtils::GLNorm[0][iy] = 0.;
0261       for (int ix = 0; ix <= MuScleFitUtils::nbins; ix++) {
0262         MuScleFitUtils::GLValue[0][ix][iy] = GL[0]->GetBinContent(ix + 1, iy + 1);
0263         // N.B. approximation: we should compute the integral of the function used to compute the probability (linear
0264         // interpolation of the mass points). This computation could be troublesome because the points have a steep
0265         // variation near the mass peak and the normal integral is not precise in these conditions.
0266         // Furthermore it is slow.
0267         MuScleFitUtils::GLNorm[0][iy] +=
0268             MuScleFitUtils::GLValue[0][ix][iy] * (2 * MuScleFitUtils::ResHalfWidth[0]) / MuScleFitUtils::nbins;
0269       }
0270     }
0271   }
0272   // Extract normalization for each mass slice
0273   // -----------------------------------------
0274   for (unsigned int ires = 1; ires < GL.size(); ires++) {
0275     if (MuScleFitUtils::resfind[ires]) {
0276       int nBinsX = GL[ires]->GetNbinsX();
0277       int nBinsY = GL[ires]->GetNbinsY();
0278       if (nBinsX != MuScleFitUtils::nbins + 1 || nBinsY != MuScleFitUtils::nbins + 1) {
0279         std::cout << "Error: for histogram \"" << GL[ires]->GetName() << "\" bins are not " << MuScleFitUtils::nbins
0280                   << std::endl;
0281         std::cout << "nBinsX = " << nBinsX << ", nBinsY = " << nBinsY << std::endl;
0282         exit(1);
0283       }
0284 
0285       for (int iy = 0; iy <= MuScleFitUtils::nbins; iy++) {
0286         MuScleFitUtils::GLNorm[ires][iy] = 0.;
0287         for (int ix = 0; ix <= MuScleFitUtils::nbins; ix++) {
0288           MuScleFitUtils::GLValue[ires][ix][iy] = GL[ires]->GetBinContent(ix + 1, iy + 1);
0289           // N.B. approximation: we should compute the integral of the function used to compute the probability (linear
0290           // interpolation of the mass points). This computation could be troublesome because the points have a steep
0291           // variation near the mass peak and the normal integral is not precise in these conditions.
0292           // Furthermore it is slow.
0293           MuScleFitUtils::GLNorm[ires][iy] +=
0294               MuScleFitUtils::GLValue[ires][ix][iy] * (2 * MuScleFitUtils::ResHalfWidth[ires]) / MuScleFitUtils::nbins;
0295         }
0296       }
0297     }
0298   }
0299   // Free all the memory for the probability histograms.
0300   if (MuScleFitUtils::rapidityBinsForZ_ && MuScleFitUtils::resfind[0] && theMuonType_ != 2) {
0301     for (unsigned int i = 0; i < GLZ.size(); i++) {
0302       delete GLZ[i];
0303     }
0304   }
0305   if (MuScleFitUtils::resfind[0] && theMuonType_ == 2)
0306     delete GL[0];
0307   for (unsigned int ires = 1; ires < GL.size(); ires++) {
0308     if (MuScleFitUtils::resfind[ires])
0309       delete GL[ires];
0310   }
0311 }
0312 
0313 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!