Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:02 

0001 // -*- C++ -*-
0002 //
0003 // Package:    MuonIdentification
0004 // Class:      MuonCaloCompatibility
0005 //
0006 /*
0007 
0008  Description: test track muon hypothesis using energy deposition in ECAL,HCAL,HO
0009 
0010 */
0011 //
0012 // Original Author:  Ingo Bloch
0013 //
0014 //
0015 #include "RecoMuon/MuonIdentification/interface/MuonCaloCompatibility.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 
0019 #include "TFile.h"
0020 
0021 void MuonCaloCompatibility::configure(const edm::ParameterSet& iConfig) {
0022   const std::string muonfileName = (iConfig.getParameter<edm::FileInPath>("MuonTemplateFileName")).fullPath();
0023   const std::string pionfileName = (iConfig.getParameter<edm::FileInPath>("PionTemplateFileName")).fullPath();
0024   TFile muon_templates(muonfileName.c_str(), "READ");
0025   TFile pion_templates(pionfileName.c_str(), "READ");
0026 
0027   pion_em_etaEmi.reset((TH2D*)pion_templates.Get("em_etaEmi"));
0028   pion_had_etaEmi.reset((TH2D*)pion_templates.Get("had_etaEmi"));
0029 
0030   pion_em_etaTmi.reset((TH2D*)pion_templates.Get("em_etaTmi"));
0031   pion_had_etaTmi.reset((TH2D*)pion_templates.Get("had_etaTmi"));
0032 
0033   pion_em_etaB.reset((TH2D*)pion_templates.Get("em_etaB"));
0034   pion_had_etaB.reset((TH2D*)pion_templates.Get("had_etaB"));
0035   pion_ho_etaB.reset((TH2D*)pion_templates.Get("ho_etaB"));
0036 
0037   pion_em_etaTpl.reset((TH2D*)pion_templates.Get("em_etaTpl"));
0038   pion_had_etaTpl.reset((TH2D*)pion_templates.Get("had_etaTpl"));
0039 
0040   pion_em_etaEpl.reset((TH2D*)pion_templates.Get("em_etaEpl"));
0041   pion_had_etaEpl.reset((TH2D*)pion_templates.Get("had_etaEpl"));
0042 
0043   muon_em_etaEmi.reset((TH2D*)muon_templates.Get("em_etaEmi"));
0044   muon_had_etaEmi.reset((TH2D*)muon_templates.Get("had_etaEmi"));
0045 
0046   muon_em_etaTmi.reset((TH2D*)muon_templates.Get("em_etaTmi"));
0047   muon_had_etaTmi.reset((TH2D*)muon_templates.Get("had_etaTmi"));
0048 
0049   muon_em_etaB.reset((TH2D*)muon_templates.Get("em_etaB"));
0050   muon_had_etaB.reset((TH2D*)muon_templates.Get("had_etaB"));
0051   muon_ho_etaB.reset((TH2D*)muon_templates.Get("ho_etaB"));
0052 
0053   muon_em_etaTpl.reset((TH2D*)muon_templates.Get("em_etaTpl"));
0054   muon_had_etaTpl.reset((TH2D*)muon_templates.Get("had_etaTpl"));
0055 
0056   muon_em_etaEpl.reset((TH2D*)muon_templates.Get("em_etaEpl"));
0057   muon_had_etaEpl.reset((TH2D*)muon_templates.Get("had_etaEpl"));
0058 
0059   // Release from the opened file
0060   pion_em_etaEmi->SetDirectory(nullptr);
0061   pion_had_etaEmi->SetDirectory(nullptr);
0062 
0063   pion_em_etaTmi->SetDirectory(nullptr);
0064   pion_had_etaTmi->SetDirectory(nullptr);
0065 
0066   pion_em_etaB->SetDirectory(nullptr);
0067   pion_had_etaB->SetDirectory(nullptr);
0068   pion_ho_etaB->SetDirectory(nullptr);
0069 
0070   pion_em_etaTpl->SetDirectory(nullptr);
0071   pion_had_etaTpl->SetDirectory(nullptr);
0072 
0073   pion_em_etaEpl->SetDirectory(nullptr);
0074   pion_had_etaEpl->SetDirectory(nullptr);
0075 
0076   muon_em_etaEmi->SetDirectory(nullptr);
0077   muon_had_etaEmi->SetDirectory(nullptr);
0078 
0079   muon_em_etaTmi->SetDirectory(nullptr);
0080   muon_had_etaTmi->SetDirectory(nullptr);
0081 
0082   muon_em_etaB->SetDirectory(nullptr);
0083   muon_had_etaB->SetDirectory(nullptr);
0084   muon_ho_etaB->SetDirectory(nullptr);
0085 
0086   muon_em_etaTpl->SetDirectory(nullptr);
0087   muon_had_etaTpl->SetDirectory(nullptr);
0088 
0089   muon_em_etaEpl->SetDirectory(nullptr);
0090   muon_had_etaEpl->SetDirectory(nullptr);
0091 
0092   // change names
0093   const std::string prefixPion = "MuonCaloCompatibility_pion_";
0094   pion_em_etaEmi->SetName((prefixPion + pion_em_etaEmi->GetName()).c_str());
0095   pion_had_etaEmi->SetName((prefixPion + pion_had_etaEmi->GetName()).c_str());
0096 
0097   pion_em_etaTmi->SetName((prefixPion + pion_em_etaTmi->GetName()).c_str());
0098   pion_had_etaTmi->SetName((prefixPion + pion_had_etaTmi->GetName()).c_str());
0099 
0100   pion_em_etaB->SetName((prefixPion + pion_em_etaB->GetName()).c_str());
0101   pion_had_etaB->SetName((prefixPion + pion_had_etaB->GetName()).c_str());
0102   pion_ho_etaB->SetName((prefixPion + pion_ho_etaB->GetName()).c_str());
0103 
0104   pion_em_etaTpl->SetName((prefixPion + pion_em_etaTpl->GetName()).c_str());
0105   pion_had_etaTpl->SetName((prefixPion + pion_had_etaTpl->GetName()).c_str());
0106 
0107   pion_em_etaEpl->SetName((prefixPion + pion_em_etaEpl->GetName()).c_str());
0108   pion_had_etaEpl->SetName((prefixPion + pion_had_etaEpl->GetName()).c_str());
0109 
0110   const std::string prefixMuon = "MuonCaloCompatibility_muon_";
0111   muon_em_etaEmi->SetName((prefixMuon + muon_em_etaEmi->GetName()).c_str());
0112   muon_had_etaEmi->SetName((prefixMuon + muon_had_etaEmi->GetName()).c_str());
0113 
0114   muon_em_etaTmi->SetName((prefixMuon + muon_em_etaTmi->GetName()).c_str());
0115   muon_had_etaTmi->SetName((prefixMuon + muon_had_etaTmi->GetName()).c_str());
0116 
0117   muon_em_etaB->SetName((prefixMuon + muon_em_etaB->GetName()).c_str());
0118   muon_had_etaB->SetName((prefixMuon + muon_had_etaB->GetName()).c_str());
0119   muon_ho_etaB->SetName((prefixMuon + muon_ho_etaB->GetName()).c_str());
0120 
0121   muon_em_etaTpl->SetName((prefixMuon + muon_em_etaTpl->GetName()).c_str());
0122   muon_had_etaTpl->SetName((prefixMuon + muon_had_etaTpl->GetName()).c_str());
0123 
0124   muon_em_etaEpl->SetName((prefixMuon + muon_em_etaEpl->GetName()).c_str());
0125   muon_had_etaEpl->SetName((prefixMuon + muon_had_etaEpl->GetName()).c_str());
0126 
0127   pbx = -1;
0128   pby = -1;
0129   pbz = -1;
0130 
0131   psx = -1;
0132   psy = -1;
0133   psz = -1;
0134 
0135   muon_compatibility = -1;
0136 
0137   use_corrected_hcal = true;
0138   use_em_special = true;
0139   isConfigured_ = true;
0140 }
0141 
0142 bool MuonCaloCompatibility::accessing_overflow(const TH2D& histo, double x, double y) {
0143   bool access = false;
0144 
0145   if (histo.GetXaxis()->FindBin(x) == 0 || histo.GetXaxis()->FindBin(x) > histo.GetXaxis()->GetNbins()) {
0146     access = true;
0147   }
0148   if (histo.GetYaxis()->FindBin(y) == 0 || histo.GetYaxis()->FindBin(y) > histo.GetYaxis()->GetNbins()) {
0149     access = true;
0150   }
0151   return access;
0152 }
0153 
0154 double MuonCaloCompatibility::evaluate(const reco::Muon& amuon) {
0155   if (!isConfigured_) {
0156     edm::LogWarning("MuonIdentification") << "MuonCaloCompatibility is not configured! Nothing is calculated.";
0157     return -9999;
0158   }
0159 
0160   double eta = 0.;
0161   double p = 0.;
0162   double em = 0.;
0163   double had = 0.;
0164   double ho = 0.;
0165 
0166   // had forgotten this reset in previous versions 070409
0167   pbx = 1.;
0168   pby = 1.;
0169   pbz = 1.;
0170 
0171   psx = 1.;
0172   psy = 1.;
0173   psz = 1.;
0174 
0175   muon_compatibility = -1.;
0176 
0177   pion_template_em = nullptr;
0178   muon_template_em = nullptr;
0179 
0180   pion_template_had = nullptr;
0181   muon_template_had = nullptr;
0182 
0183   pion_template_ho = nullptr;
0184   muon_template_ho = nullptr;
0185 
0186   // 071002: Get either tracker track, or SAmuon track.
0187   // CaloCompatibility templates may have to be specialized for
0188   // the use with SAmuons, currently just using the ones produced
0189   // using tracker tracks.
0190   const reco::Track* track = nullptr;
0191   if (!amuon.track().isNull()) {
0192     track = amuon.track().get();
0193   } else {
0194     if (!amuon.standAloneMuon().isNull()) {
0195       track = amuon.standAloneMuon().get();
0196     } else {
0197       throw cms::Exception("FatalError")
0198           << "Failed to fill muon id calo_compatibility information for a muon with undefined references to tracks";
0199     }
0200   }
0201 
0202   if (!use_corrected_hcal) {  // old eta regions, uncorrected energy
0203     eta = track->eta();
0204     p = track->p();
0205 
0206     // new 070904: Set lookup momentum to 1999.9 if larger than 2 TeV.
0207     // Though the templates were produced with p<2TeV, we believe that
0208     // this approximation should be roughly valid. A special treatment
0209     // for >1 TeV muons is advisable anyway :)
0210     if (p >= 2000.)
0211       p = 1999.9;
0212 
0213     //    p   = 10./sin(track->theta()); // use this for templates < 1_5
0214     if (use_em_special) {
0215       if (amuon.calEnergy().em == 0.)
0216         em = -5.;
0217       else
0218         em = amuon.calEnergy().em;
0219     } else {
0220       em = amuon.calEnergy().em;
0221     }
0222     had = amuon.calEnergy().had;
0223     ho = amuon.calEnergy().ho;
0224   } else {
0225     eta = track->eta();
0226     p = track->p();
0227 
0228     // new 070904: Set lookup momentum to 1999.9 if larger than 2 TeV.
0229     // Though the templates were produced with p<2TeV, we believe that
0230     // this approximation should be roughly valid. A special treatment
0231     // for >1 TeV muons is advisable anyway :)
0232     if (p >= 2000.)
0233       p = 1999.9;
0234 
0235     //    p   = 10./sin(track->theta());  // use this for templates < 1_5
0236     // hcal energy is now done where we get the template histograms (to use corrected cal energy)!
0237     //     had = amuon.calEnergy().had;
0238     if (use_em_special) {
0239       if (amuon.calEnergy().em == 0.)
0240         em = -5.;
0241       else
0242         em = amuon.calEnergy().em;
0243     } else {
0244       em = amuon.calEnergy().em;
0245     }
0246     ho = amuon.calEnergy().ho;
0247   }
0248 
0249   // Skip everyting and return "I don't know" (i.e. 0.5) for uncovered regions:
0250   //  if( p < 0. || p > 500.) return 0.5; // removed 500 GeV cutoff 070817 after updating the tempates (v2_0) to have valid entried beyond 500 GeV
0251   if (p < 0.)
0252     return 0.5;  // return "unknown" for unphysical momentum input.
0253   if (fabs(eta) > 2.5)
0254     return 0.5;
0255   // temporary fix for low association efficiency:
0256   // set caloCompatibility to 0.12345 for tracks
0257   // which have 0 energy in BOTH ecal and hcal
0258   if (amuon.calEnergy().had == 0.0 && amuon.calEnergy().em == 0.0)
0259     return 0.12345;
0260 
0261   //  std::cout<<std::endl<<"Input values are: "<<eta <<" "<< p <<" "<< em <<" "<< had <<" "<< ho;
0262 
0263   //  depending on the eta, choose correct histogram: (now all for barrel):
0264   // bad! eta range has to be syncronised with choice for histogram... should be read out from the histo file somehow... 070322
0265   if (42 != 42) {  // old eta ranges and uncorrected hcal energy
0266     if (eta <= -1.4) {
0267       //    std::cout<<"Emi"<<std::endl;
0268       pion_template_em = pion_em_etaEmi;
0269       pion_template_had = pion_had_etaEmi;
0270       muon_template_em = muon_em_etaEmi;
0271       muon_template_had = muon_had_etaEmi;
0272     } else if (eta > -1.4 && eta <= -1.31) {
0273       //    std::cout<<"Tmi"<<std::endl;
0274       pion_template_em = pion_em_etaTmi;
0275       pion_template_had = pion_had_etaTmi;
0276       muon_template_em = muon_em_etaTmi;
0277       muon_template_had = muon_had_etaTmi;
0278     } else if (eta > -1.31 && eta <= 1.31) {
0279       //    std::cout<<"B"<<std::endl;
0280       pion_template_em = pion_em_etaB;
0281       pion_template_had = pion_had_etaB;
0282       pion_template_ho = pion_ho_etaB;
0283       muon_template_em = muon_em_etaB;
0284       muon_template_had = muon_had_etaB;
0285       muon_template_ho = muon_ho_etaB;
0286     } else if (eta > 1.31 && eta <= 1.4) {
0287       //    std::cout<<"Tpl"<<std::endl;
0288       pion_template_em = pion_em_etaTpl;
0289       pion_template_had = pion_had_etaTpl;
0290       muon_template_em = muon_em_etaTpl;
0291       muon_template_had = muon_had_etaTpl;
0292     } else if (eta > 1.4) {
0293       //    std::cout<<"Epl"<<std::endl;
0294       pion_template_em = pion_em_etaEpl;
0295       pion_template_had = pion_had_etaEpl;
0296       muon_template_em = muon_em_etaEpl;
0297       muon_template_had = muon_had_etaEpl;
0298     } else {
0299       LogTrace("MuonIdentification")
0300           << "Some very weird thing happened in MuonCaloCompatibility::evaluate - go figure ;) ";
0301       return -999;
0302     }
0303   } else if (42 == 42) {  // new eta bins, corrected hcal energy
0304     if (track->eta() > 1.27) {
0305       //    had_etaEpl ->Fill(muon->track().get()->p(),1.8/2.2*muon->calEnergy().had );
0306       if (use_corrected_hcal)
0307         had = 1.8 / 2.2 * amuon.calEnergy().had;
0308       else
0309         had = amuon.calEnergy().had;
0310       pion_template_had = pion_had_etaEpl;
0311       muon_template_had = muon_had_etaEpl;
0312     }
0313     if (track->eta() <= 1.27 && track->eta() > 1.1) {
0314       //    had_etaTpl ->Fill(muon->track().get()->p(),(1.8/(-2.2*muon->track().get()->eta()+5.5))*muon->calEnergy().had );
0315       if (use_corrected_hcal)
0316         had = (1.8 / (-2.2 * track->eta() + 5.5)) * amuon.calEnergy().had;
0317       else
0318         had = amuon.calEnergy().had;
0319       pion_template_had = pion_had_etaTpl;
0320       muon_template_had = muon_had_etaTpl;
0321     }
0322     if (track->eta() <= 1.1 && track->eta() > -1.1) {
0323       //    had_etaB   ->Fill(muon->track().get()->p(),sin(muon->track().get()->theta())*muon->calEnergy().had );
0324       if (use_corrected_hcal)
0325         had = sin(track->theta()) * amuon.calEnergy().had;
0326       else
0327         had = amuon.calEnergy().had;
0328       pion_template_had = pion_had_etaB;
0329       muon_template_had = muon_had_etaB;
0330     }
0331     if (track->eta() <= -1.1 && track->eta() > -1.27) {
0332       //    had_etaTmi ->Fill(muon->track().get()->p(),(1.8/(-2.2*muon->track().get()->eta()+5.5))*muon->calEnergy().had );
0333       if (use_corrected_hcal)
0334         had = (1.8 / (2.2 * track->eta() + 5.5)) * amuon.calEnergy().had;
0335       else
0336         had = amuon.calEnergy().had;
0337       pion_template_had = pion_had_etaTmi;
0338       muon_template_had = muon_had_etaTmi;
0339     }
0340     if (track->eta() <= -1.27) {
0341       //    had_etaEmi ->Fill(muon->track().get()->p(),1.8/2.2*muon->calEnergy().had );
0342       if (use_corrected_hcal)
0343         had = 1.8 / 2.2 * amuon.calEnergy().had;
0344       else
0345         had = amuon.calEnergy().had;
0346       pion_template_had = pion_had_etaEmi;
0347       muon_template_had = muon_had_etaEmi;
0348     }
0349 
0350     // just two eta regions for Ecal (+- 1.479 for barrel, else for rest), no correction:
0351 
0352     //    std::cout<<"We have a muon with an eta of: "<<track->eta()<<std::endl;
0353 
0354     if (track->eta() > 1.479) {
0355       //    em_etaEpl  ->Fill(muon->track().get()->p(),muon->calEnergy().em   );
0356       //    //  em_etaTpl  ->Fill(muon->track().get()->p(),muon->calEnergy().em   );
0357       ////      em  = amuon.calEnergy().em;
0358       pion_template_em = pion_em_etaEpl;
0359       muon_template_em = muon_em_etaEpl;
0360     }
0361     if (track->eta() <= 1.479 && track->eta() > -1.479) {
0362       //    em_etaB    ->Fill(muon->track().get()->p(),muon->calEnergy().em   );
0363       ////      em  = amuon.calEnergy().em;
0364       pion_template_em = pion_em_etaB;
0365       muon_template_em = muon_em_etaB;
0366     }
0367     if (track->eta() <= -1.479) {
0368       //    //  em_etaTmi  ->Fill(muon->track().get()->p(),muon->calEnergy().em   );
0369       //    em_etaEmi  ->Fill(muon->track().get()->p(),muon->calEnergy().em   );
0370       ////      em  = amuon.calEnergy().em;
0371       pion_template_em = pion_em_etaEmi;
0372       muon_template_em = muon_em_etaEmi;
0373     }
0374 
0375     // just one barrel eta region for the HO, no correction
0376     //    if( track->eta() < 1.4 && track->eta() > -1.4 ) { // experimenting now...
0377     if (track->eta() < 1.28 && track->eta() > -1.28) {
0378       //    ho_etaB    ->Fill(muon->track().get()->p(),muon->calEnergy().ho   );
0379       ////      ho  = amuon.calEnergy().ho;
0380       pion_template_ho = pion_ho_etaB;
0381       muon_template_ho = muon_ho_etaB;
0382     }
0383   }
0384 
0385   if (42 != 42) {  // check validity of input template histos and input variables"
0386     pion_template_em->ls();
0387     pion_template_had->ls();
0388     if (pion_template_ho)
0389       pion_template_ho->ls();
0390     muon_template_em->ls();
0391     muon_template_had->ls();
0392     if (muon_template_ho)
0393       muon_template_ho->ls();
0394 
0395     LogTrace("MuonIdentification") << "Input variables: eta    p     em     had    ho "
0396                                    << "\n"
0397                                    << eta << " " << p << " " << em << " " << had << " " << ho << " "
0398                                    << "\n"
0399                                    << "cal uncorr:    em     had    ho "
0400                                    << "\n"
0401                                    << eta << " " << p << " " << amuon.calEnergy().em << " " << amuon.calEnergy().had
0402                                    << " " << amuon.calEnergy().ho;
0403   }
0404 
0405   //  Look up Compatibility by, where x is p and y the energy.
0406   //  We have a set of different histograms for different regions of eta.
0407 
0408   // need error meassage in case the template histos are missing / the template file is not present!!! 070412
0409 
0410   if (pion_template_em) {  // access ecal background template
0411     if (accessing_overflow(*pion_template_em, p, em)) {
0412       pbx = 1.;
0413       psx = 1.;
0414       LogTrace("MuonIdentification") << "            // Message: trying to access overflow bin in MuonCompatibility "
0415                                         "template for ecal - defaulting signal and background  ";
0416       LogTrace("MuonIdentification") << "            // template value to 1. " << pion_template_em->GetName()
0417                                      << " e: " << em << " p: " << p;
0418     } else
0419       pbx = pion_template_em->GetBinContent(pion_template_em->GetXaxis()->FindBin(p),
0420                                             pion_template_em->GetYaxis()->FindBin(em));
0421   }
0422   if (pion_template_had) {  // access hcal background template
0423     if (accessing_overflow(*pion_template_had, p, had)) {
0424       pby = 1.;
0425       psy = 1.;
0426       LogTrace("MuonIdentification") << "            // Message: trying to access overflow bin in MuonCompatibility "
0427                                         "template for hcal - defaulting signal and background  ";
0428       LogTrace("MuonIdentification") << "            // template value to 1. " << pion_template_had->GetName()
0429                                      << " e: " << had << " p: " << p;
0430     } else
0431       pby = pion_template_had->GetBinContent(pion_template_had->GetXaxis()->FindBin(p),
0432                                              pion_template_had->GetYaxis()->FindBin(had));
0433   }
0434   if (pion_template_ho) {  // access ho background template
0435     if (accessing_overflow(*pion_template_ho, p, ho)) {
0436       pbz = 1.;
0437       psz = 1.;
0438       LogTrace("MuonIdentification") << "            // Message: trying to access overflow bin in MuonCompatibility "
0439                                         "template for ho   - defaulting signal and background  ";
0440       LogTrace("MuonIdentification") << "            // template value to 1. " << pion_template_ho->GetName()
0441                                      << " e: " << em << " p: " << p;
0442     } else
0443       pbz = pion_template_ho->GetBinContent(pion_template_ho->GetXaxis()->FindBin(p),
0444                                             pion_template_ho->GetYaxis()->FindBin(ho));
0445   }
0446 
0447   if (muon_template_em) {  // access ecal background template
0448     if (accessing_overflow(*muon_template_em, p, em)) {
0449       psx = 1.;
0450       pbx = 1.;
0451       LogTrace("MuonIdentification") << "            // Message: trying to access overflow bin in MuonCompatibility "
0452                                         "template for ecal - defaulting signal and background  ";
0453       LogTrace("MuonIdentification") << "            // template value to 1. " << muon_template_em->GetName()
0454                                      << " e: " << em << " p: " << p;
0455     } else
0456       psx = muon_template_em->GetBinContent(muon_template_em->GetXaxis()->FindBin(p),
0457                                             muon_template_em->GetYaxis()->FindBin(em));
0458   }
0459   if (muon_template_had) {  // access hcal background template
0460     if (accessing_overflow(*muon_template_had, p, had)) {
0461       psy = 1.;
0462       pby = 1.;
0463       LogTrace("MuonIdentification") << "            // Message: trying to access overflow bin in MuonCompatibility "
0464                                         "template for hcal - defaulting signal and background  ";
0465       LogTrace("MuonIdentification") << "            // template value to 1. " << muon_template_had->GetName()
0466                                      << " e: " << had << " p: " << p;
0467     } else
0468       psy = muon_template_had->GetBinContent(muon_template_had->GetXaxis()->FindBin(p),
0469                                              muon_template_had->GetYaxis()->FindBin(had));
0470   }
0471   if (muon_template_ho) {  // access ho background template
0472     if (accessing_overflow(*muon_template_ho, p, ho)) {
0473       psz = 1.;
0474       pbz = 1.;
0475       LogTrace("MuonIdentification") << "            // Message: trying to access overflow bin in MuonCompatibility "
0476                                         "template for ho   - defaulting signal and background  ";
0477       LogTrace("MuonIdentification") << "            // template value to 1. " << muon_template_ho->GetName()
0478                                      << " e: " << ho << " p: " << p;
0479     } else
0480       psz = muon_template_ho->GetBinContent(muon_template_ho->GetXaxis()->FindBin(p),
0481                                             muon_template_ho->GetYaxis()->FindBin(ho));
0482   }
0483 
0484   // erm - what is this?!?! How could the HO probability be less than 0????? Do we want this line!?!?
0485   if (psz <= 0.)
0486     psz = 1.;
0487   if (pbz <= 0.)
0488     pbz = 1.;
0489 
0490   // Protection agains empty bins - set cal part to neutral if the bin of the template is empty
0491   // (temporary fix, a proper extrapolation would be better)
0492   if (psx == 0. || pbx == 0.) {
0493     psx = 1.;
0494     pbx = 1.;
0495   }
0496   if (psy == 0. || pby == 0.) {
0497     psy = 1.;
0498     pby = 1.;
0499   }
0500   if (psz == 0. || pbz == 0.) {
0501     psz = 1.;
0502     pbz = 1.;
0503   }
0504 
0505   // There are two classes of events that deliver 0 for the hcal or ho energy:
0506   // 1) the track momentum is so low that the extrapolation tells us it should not have reached the cal
0507   // 2) the crossed cell had an reading below the readout cuts.
0508   // The 2nd case discriminates between muons and hadrons, the 1st not. Thus for the time being,
0509   // we set the returned ps and pb to 0 for these cases.
0510   // We need to have a return value different from 0 for the 1st case in the long run.
0511   if (had == 0.0) {
0512     psy = 1.;
0513     pby = 1.;
0514   }
0515   if (ho == 0.0) {
0516     psz = 1.;
0517     pbz = 1.;
0518   }
0519 
0520   // Set em to neutral if no energy in em or negative energy measured.
0521   // (These cases might indicate problems in the ecal association or readout?! The only
0522   // hint so far: for critical eta region (eta in [1.52, 1.64]) have negative em values.)
0523   if (em <= 0. && !use_em_special) {
0524     pbx = 1.;
0525     psx = 1.;
0526   }
0527 
0528   if ((psx * psy * psz + pbx * pby * pbz) > 0.)
0529     muon_compatibility = psx * psy * psz / (psx * psy * psz + pbx * pby * pbz);
0530   else {
0531     LogTrace("MuonIdentification") << "Divide by 0 - defaulting consistency to 0.5 (neutral)!!";
0532     muon_compatibility = 0.5;
0533     LogTrace("MuonIdentification") << "Input variables: eta    p     em     had    ho "
0534                                    << "\n"
0535                                    << eta << " " << p << " " << em << " " << had << " " << ho << " "
0536                                    << "\n"
0537                                    << "cal uncorr:    em     had    ho "
0538                                    << "\n"
0539                                    << eta << " " << p << " " << amuon.calEnergy().em << " " << amuon.calEnergy().had
0540                                    << " " << amuon.calEnergy().ho;
0541   }
0542   return muon_compatibility;
0543 }
[ 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!