CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoMuon/MuonIdentification/src/MuonCaloCompatibility.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    MuonIdentification
00004 // Class:      MuonCaloCompatibility
00005 // 
00006 /*
00007 
00008  Description: test track muon hypothesis using energy deposition in ECAL,HCAL,HO
00009 
00010 */
00011 //
00012 // Original Author:  Ingo Bloch
00013 // $Id: MuonCaloCompatibility.cc,v 1.8 2010/06/02 15:50:51 andersj Exp $
00014 //
00015 //
00016 #include "RecoMuon/MuonIdentification/interface/MuonCaloCompatibility.h"
00017 #include "DataFormats/TrackReco/interface/Track.h"
00018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00019 
00020 void MuonCaloCompatibility::configure(const edm::ParameterSet& iConfig)
00021 {
00022    MuonfileName_ = (iConfig.getParameter<edm::FileInPath>("MuonTemplateFileName")).fullPath();
00023    PionfileName_ = (iConfig.getParameter<edm::FileInPath>("PionTemplateFileName")).fullPath();
00024    muon_templates.reset( new TFile(MuonfileName_.c_str(),"READ") );
00025    pion_templates.reset( new TFile(PionfileName_.c_str(),"READ") );
00026 
00027    pion_em_etaEmi  = (TH2D*) pion_templates->Get("em_etaEmi");
00028    pion_had_etaEmi = (TH2D*) pion_templates->Get("had_etaEmi");
00029                 
00030    pion_em_etaTmi  = (TH2D*) pion_templates->Get("em_etaTmi");
00031    pion_had_etaTmi = (TH2D*) pion_templates->Get("had_etaTmi");
00032                 
00033    pion_em_etaB    = (TH2D*) pion_templates->Get("em_etaB");
00034    pion_had_etaB   = (TH2D*) pion_templates->Get("had_etaB");
00035    pion_ho_etaB    = (TH2D*) pion_templates->Get("ho_etaB");
00036    
00037    pion_em_etaTpl  = (TH2D*) pion_templates->Get("em_etaTpl");
00038    pion_had_etaTpl = (TH2D*) pion_templates->Get("had_etaTpl");
00039                 
00040    pion_em_etaEpl  = (TH2D*) pion_templates->Get("em_etaEpl");
00041    pion_had_etaEpl = (TH2D*) pion_templates->Get("had_etaEpl");
00042                 
00043    muon_em_etaEmi  = (TH2D*) muon_templates->Get("em_etaEmi");
00044    muon_had_etaEmi = (TH2D*) muon_templates->Get("had_etaEmi");
00045                 
00046    muon_em_etaTmi  = (TH2D*) muon_templates->Get("em_etaTmi");
00047    muon_had_etaTmi = (TH2D*) muon_templates->Get("had_etaTmi");
00048                 
00049    muon_em_etaB    = (TH2D*) muon_templates->Get("em_etaB");
00050    muon_had_etaB   = (TH2D*) muon_templates->Get("had_etaB");
00051    muon_ho_etaB    = (TH2D*) muon_templates->Get("ho_etaB");
00052                 
00053    muon_em_etaTpl  = (TH2D*) muon_templates->Get("em_etaTpl");
00054    muon_had_etaTpl = (TH2D*) muon_templates->Get("had_etaTpl");
00055                 
00056    muon_em_etaEpl  = (TH2D*) muon_templates->Get("em_etaEpl");
00057    muon_had_etaEpl = (TH2D*) muon_templates->Get("had_etaEpl");
00058 
00059    pbx = -1;
00060    pby = -1;
00061    pbz = -1;
00062 
00063    psx = -1;
00064    psy = -1;
00065    psz = -1;
00066 
00067    muon_compatibility = -1;
00068 
00069    use_corrected_hcal = true;
00070    use_em_special = true;
00071    isConfigured_ = true;
00072 }
00073 
00074 bool MuonCaloCompatibility::accessing_overflow( TH2D* histo, double x, double y ) {
00075   bool access = false;
00076 
00077   if( histo->GetXaxis()->FindBin(x) == 0 || 
00078       histo->GetXaxis()->FindBin(x) > histo->GetXaxis()->GetNbins() ) {
00079     access = true;
00080   }
00081   if( histo->GetYaxis()->FindBin(y) == 0 || 
00082       histo->GetYaxis()->FindBin(y) > histo->GetYaxis()->GetNbins() ) {
00083     access = true;
00084   }
00085   return access;
00086 }
00087 
00088 double MuonCaloCompatibility::evaluate( const reco::Muon& amuon ) {
00089   if (! isConfigured_) {
00090      edm::LogWarning("MuonIdentification") << "MuonCaloCompatibility is not configured! Nothing is calculated.";
00091      return -9999;
00092   }
00093    
00094   double eta = 0.;
00095   double p   = 0.;
00096   double em  = 0.;
00097   double had = 0.;
00098   double ho  = 0.;
00099 
00100   // had forgotten this reset in previous versions 070409
00101   pbx = 1.;
00102   pby = 1.;
00103   pbz = 1.;
00104 
00105   psx = 1.;
00106   psy = 1.;
00107   psz = 1.;
00108 
00109   muon_compatibility = -1.;
00110   
00111   pion_template_em   = NULL;
00112   muon_template_em   = NULL;
00113   
00114   pion_template_had  = NULL;
00115   muon_template_had  = NULL;
00116   
00117   pion_template_ho   = NULL;
00118   muon_template_ho   = NULL;
00119   
00120   // 071002: Get either tracker track, or SAmuon track.
00121   // CaloCompatibility templates may have to be specialized for 
00122   // the use with SAmuons, currently just using the ones produced
00123   // using tracker tracks. 
00124   const reco::Track* track = 0;
00125   if ( ! amuon.track().isNull() ) {
00126     track = amuon.track().get();
00127   }
00128   else {
00129     if ( ! amuon.standAloneMuon().isNull() ) {
00130       track = amuon.standAloneMuon().get();
00131     }
00132     else {
00133       throw cms::Exception("FatalError") << "Failed to fill muon id calo_compatibility information for a muon with undefined references to tracks"; 
00134     }
00135   }
00136 
00137   if( !use_corrected_hcal ) { // old eta regions, uncorrected energy
00138     eta = track->eta();
00139     p   = track->p(); 
00140 
00141     // new 070904: Set lookup momentum to 1999.9 if larger than 2 TeV. 
00142     // Though the templates were produced with p<2TeV, we believe that
00143     // this approximation should be roughly valid. A special treatment
00144     // for >1 TeV muons is advisable anyway :)
00145     if( p>=2000. ) p = 1999.9;
00146 
00147     //    p   = 10./sin(track->theta()); // use this for templates < 1_5
00148     if( use_em_special ) {
00149       if( amuon.calEnergy().em == 0. )    em  = -5.;
00150       else em  = amuon.calEnergy().em;
00151     }
00152     else {
00153       em  = amuon.calEnergy().em;
00154     }
00155     had = amuon.calEnergy().had;
00156     ho  = amuon.calEnergy().ho;
00157   }
00158   else {
00159     eta = track->eta();
00160     p   = track->p();
00161     
00162     // new 070904: Set lookup momentum to 1999.9 if larger than 2 TeV. 
00163     // Though the templates were produced with p<2TeV, we believe that
00164     // this approximation should be roughly valid. A special treatment
00165     // for >1 TeV muons is advisable anyway :)
00166     if( p>=2000. ) p = 1999.9;
00167 
00168     //    p   = 10./sin(track->theta());  // use this for templates < 1_5
00169     // hcal energy is now done where we get the template histograms (to use corrected cal energy)!
00170     //     had = amuon.calEnergy().had;
00171     if( use_em_special ) {
00172       if( amuon.calEnergy().em == 0. )    em  = -5.;
00173       else em  = amuon.calEnergy().em;
00174     }
00175     else {
00176       em  = amuon.calEnergy().em;
00177     }
00178     ho  = amuon.calEnergy().ho;
00179   }
00180 
00181 
00182   // Skip everyting and return "I don't know" (i.e. 0.5) for uncovered regions:
00183   //  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
00184   if( p < 0. ) return 0.5; // return "unknown" for unphysical momentum input.
00185   if( fabs(eta) >  2.5 ) return 0.5; 
00186   // temporary fix for low association efficiency:
00187   // set caloCompatibility to 0.12345 for tracks
00188   // which have 0 energy in BOTH ecal and hcal
00189   if( amuon.calEnergy().had == 0.0 && amuon.calEnergy().em == 0.0 ) return 0.12345; 
00190 
00191   //  std::cout<<std::endl<<"Input values are: "<<eta <<" "<< p <<" "<< em <<" "<< had <<" "<< ho;
00192 
00193   //  depending on the eta, choose correct histogram: (now all for barrel):
00194   // bad! eta range has to be syncronised with choice for histogram... should be read out from the histo file somehow... 070322
00195   if(42 != 42) { // old eta ranges and uncorrected hcal energy
00196     if(eta <= -1.4) {
00197       //    std::cout<<"Emi"<<std::endl;
00198       pion_template_em  = pion_em_etaEmi;
00199       pion_template_had = pion_had_etaEmi;
00200       muon_template_em  = muon_em_etaEmi;
00201       muon_template_had = muon_had_etaEmi;
00202     }
00203     else if(eta > -1.4 && eta <= -1.31) {
00204       //    std::cout<<"Tmi"<<std::endl;
00205       pion_template_em  = pion_em_etaTmi;
00206       pion_template_had = pion_had_etaTmi;
00207       muon_template_em  = muon_em_etaTmi;
00208       muon_template_had = muon_had_etaTmi;
00209     }
00210     else if(eta > -1.31 && eta <= 1.31) {
00211       //    std::cout<<"B"<<std::endl;
00212       pion_template_em  = pion_em_etaB;
00213       pion_template_had = pion_had_etaB;
00214       pion_template_ho  = pion_ho_etaB;
00215       muon_template_em  = muon_em_etaB;
00216       muon_template_had = muon_had_etaB;
00217       muon_template_ho  = muon_ho_etaB;
00218     }
00219     else if(eta > 1.31 && eta <= 1.4) {
00220       //    std::cout<<"Tpl"<<std::endl;
00221       pion_template_em  = pion_em_etaTpl;
00222       pion_template_had = pion_had_etaTpl;
00223       muon_template_em  = muon_em_etaTpl;
00224       muon_template_had = muon_had_etaTpl;
00225     }
00226     else if(eta > 1.4) {
00227       //    std::cout<<"Epl"<<std::endl;
00228       pion_template_em  = pion_em_etaEpl;
00229       pion_template_had = pion_had_etaEpl;
00230       muon_template_em  = muon_em_etaEpl;
00231       muon_template_had = muon_had_etaEpl;
00232     }
00233     else {
00234       LogTrace("MuonIdentification")<<"Some very weird thing happened in MuonCaloCompatibility::evaluate - go figure ;) ";
00235       return -999;
00236     }
00237   }
00238   else if( 42 == 42 ) { // new eta bins, corrected hcal energy
00239     if(  track->eta() >  1.27  ) {
00240       //        had_etaEpl ->Fill(muon->track().get()->p(),1.8/2.2*muon->calEnergy().had );
00241       if(use_corrected_hcal)    had = 1.8/2.2*amuon.calEnergy().had;
00242       else      had = amuon.calEnergy().had;
00243       pion_template_had = pion_had_etaEpl;
00244       muon_template_had = muon_had_etaEpl;
00245     }
00246     if( track->eta() <=  1.27  && track->eta() >  1.1 ) {
00247       //        had_etaTpl ->Fill(muon->track().get()->p(),(1.8/(-2.2*muon->track().get()->eta()+5.5))*muon->calEnergy().had );
00248       if(use_corrected_hcal)    had = (1.8/(-2.2*track->eta()+5.5))*amuon.calEnergy().had;
00249       else      had = amuon.calEnergy().had;
00250       pion_template_had  = pion_had_etaTpl;
00251       muon_template_had  = muon_had_etaTpl;
00252     }
00253     if( track->eta() <=  1.1 && track->eta() > -1.1 ) {
00254       //        had_etaB   ->Fill(muon->track().get()->p(),sin(muon->track().get()->theta())*muon->calEnergy().had );
00255       if(use_corrected_hcal)    had = sin(track->theta())*amuon.calEnergy().had;
00256       else      had = amuon.calEnergy().had;
00257       pion_template_had  = pion_had_etaB;
00258       muon_template_had  = muon_had_etaB;
00259     }
00260     if( track->eta() <= -1.1 && track->eta() > -1.27 ) {
00261       //        had_etaTmi ->Fill(muon->track().get()->p(),(1.8/(-2.2*muon->track().get()->eta()+5.5))*muon->calEnergy().had );
00262       if(use_corrected_hcal)    had = (1.8/(2.2*track->eta()+5.5))*amuon.calEnergy().had;
00263       else      had = amuon.calEnergy().had;
00264       pion_template_had = pion_had_etaTmi;
00265       muon_template_had = muon_had_etaTmi;
00266     }
00267     if( track->eta() <= -1.27 ) {
00268       //        had_etaEmi ->Fill(muon->track().get()->p(),1.8/2.2*muon->calEnergy().had );
00269       if(use_corrected_hcal)    had = 1.8/2.2*amuon.calEnergy().had;
00270       else      had = amuon.calEnergy().had;
00271       pion_template_had = pion_had_etaEmi;
00272       muon_template_had = muon_had_etaEmi;
00273     }
00274     
00275     // just two eta regions for Ecal (+- 1.479 for barrel, else for rest), no correction:
00276 
00277     //    std::cout<<"We have a muon with an eta of: "<<track->eta()<<std::endl;
00278 
00279     if(  track->eta() >  1.479  ) {
00280       //        em_etaEpl  ->Fill(muon->track().get()->p(),muon->calEnergy().em   );
00281       //        //      em_etaTpl  ->Fill(muon->track().get()->p(),muon->calEnergy().em   );
00283       pion_template_em  = pion_em_etaEpl;
00284       muon_template_em  = muon_em_etaEpl;
00285     }
00286     if( track->eta() <=  1.479 && track->eta() > -1.479 ) {
00287       //        em_etaB    ->Fill(muon->track().get()->p(),muon->calEnergy().em   );
00289       pion_template_em  = pion_em_etaB;
00290       muon_template_em  = muon_em_etaB;
00291     }
00292     if( track->eta() <= -1.479 ) {
00293       //        //      em_etaTmi  ->Fill(muon->track().get()->p(),muon->calEnergy().em   );
00294       //        em_etaEmi  ->Fill(muon->track().get()->p(),muon->calEnergy().em   );
00296       pion_template_em  = pion_em_etaEmi;
00297       muon_template_em  = muon_em_etaEmi;
00298     }
00299     
00300     // just one barrel eta region for the HO, no correction
00301     //    if( track->eta() < 1.4 && track->eta() > -1.4 ) { // experimenting now...
00302     if( track->eta() < 1.28 && track->eta() > -1.28 ) {
00303       //        ho_etaB    ->Fill(muon->track().get()->p(),muon->calEnergy().ho   );
00305       pion_template_ho  = pion_ho_etaB;
00306       muon_template_ho  = muon_ho_etaB;
00307     }
00308 
00309     
00310   }
00311 
00312 
00313   if( 42 != 42 ) { // check validity of input template histos and input variables"
00314     pion_template_em ->ls();
00315     pion_template_had->ls();
00316     if(pion_template_ho)   pion_template_ho ->ls();
00317     muon_template_em ->ls();
00318     muon_template_had->ls();
00319     if(muon_template_ho)   muon_template_ho ->ls();
00320 
00321     LogTrace("MuonIdentification")<<"Input variables: eta    p     em     had    ho "<<"\n"
00322              <<eta<<" "<<p<<" "<<em<<" "<<had<<" "<<ho<<" "<<"\n"
00323              <<"cal uncorr:    em     had    ho "<<"\n"
00324              <<eta<<" "<<p<<" "<<amuon.calEnergy().em<<" "<<amuon.calEnergy().had<<" "<<amuon.calEnergy().ho;
00325   }
00326   
00327 
00328   //  Look up Compatibility by, where x is p and y the energy. 
00329   //  We have a set of different histograms for different regions of eta.
00330 
00331   // need error meassage in case the template histos are missing / the template file is not present!!! 070412
00332 
00333   if( pion_template_em )  { // access ecal background template
00334     if( accessing_overflow( pion_template_em, p, em ) ) {
00335       pbx = 1.;
00336       psx = 1.;
00337       LogTrace("MuonIdentification")<<"            // Message: trying to access overflow bin in MuonCompatibility template for ecal - defaulting signal and background  ";
00338       LogTrace("MuonIdentification")<<"            // template value to 1. "<<pion_template_em->GetName()<<" e: "<<em<<" p: "<<p;
00339     }
00340     else pbx =  pion_template_em->GetBinContent(  pion_template_em->GetXaxis()->FindBin(p), pion_template_em->GetYaxis()->FindBin(em) );
00341   }
00342   if( pion_template_had ) { // access hcal background template
00343     if( accessing_overflow( pion_template_had, p, had ) ) {
00344       pby = 1.;
00345       psy = 1.;
00346       LogTrace("MuonIdentification")<<"            // Message: trying to access overflow bin in MuonCompatibility template for hcal - defaulting signal and background  ";
00347       LogTrace("MuonIdentification")<<"            // template value to 1. "<<pion_template_had->GetName()<<" e: "<<had<<" p: "<<p;
00348     }
00349     else pby =  pion_template_had->GetBinContent(  pion_template_had->GetXaxis()->FindBin(p), pion_template_had->GetYaxis()->FindBin(had) );
00350   }
00351   if(pion_template_ho) { // access ho background template
00352     if( accessing_overflow( pion_template_ho, p, ho ) ) {
00353       pbz = 1.;
00354       psz = 1.;
00355       LogTrace("MuonIdentification")<<"            // Message: trying to access overflow bin in MuonCompatibility template for ho   - defaulting signal and background  ";
00356       LogTrace("MuonIdentification")<<"            // template value to 1. "<<pion_template_ho->GetName()<<" e: "<<em<<" p: "<<p; 
00357     }
00358     else pbz =  pion_template_ho->GetBinContent(  pion_template_ho->GetXaxis()->FindBin(p), pion_template_ho->GetYaxis()->FindBin(ho) );
00359   }
00360 
00361 
00362   if( muon_template_em )  { // access ecal background template
00363     if( accessing_overflow( muon_template_em, p, em ) ) {
00364       psx = 1.;
00365       pbx = 1.;
00366       LogTrace("MuonIdentification")<<"            // Message: trying to access overflow bin in MuonCompatibility template for ecal - defaulting signal and background  ";
00367       LogTrace("MuonIdentification")<<"            // template value to 1. "<<muon_template_em->GetName()<<" e: "<<em<<" p: "<<p;
00368     }
00369     else psx =  muon_template_em->GetBinContent(  muon_template_em->GetXaxis()->FindBin(p), muon_template_em->GetYaxis()->FindBin(em) );
00370   }
00371   if( muon_template_had ) { // access hcal background template
00372     if( accessing_overflow( muon_template_had, p, had ) ) {
00373       psy = 1.;
00374       pby = 1.;
00375       LogTrace("MuonIdentification")<<"            // Message: trying to access overflow bin in MuonCompatibility template for hcal - defaulting signal and background  ";
00376       LogTrace("MuonIdentification")<<"            // template value to 1. "<<muon_template_had->GetName()<<" e: "<<had<<" p: "<<p;
00377     }
00378     else psy =  muon_template_had->GetBinContent(  muon_template_had->GetXaxis()->FindBin(p), muon_template_had->GetYaxis()->FindBin(had) );
00379   }
00380   if(muon_template_ho) { // access ho background template
00381     if( accessing_overflow( muon_template_ho, p, ho ) ) {
00382       psz = 1.;
00383       pbz = 1.;
00384        LogTrace("MuonIdentification")<<"            // Message: trying to access overflow bin in MuonCompatibility template for ho   - defaulting signal and background  ";
00385        LogTrace("MuonIdentification")<<"            // template value to 1. "<<muon_template_ho->GetName()<<" e: "<<ho<<" p: "<<p;
00386     }
00387     else psz =  muon_template_ho->GetBinContent(  muon_template_ho->GetXaxis()->FindBin(p), muon_template_ho->GetYaxis()->FindBin(ho) );
00388   }
00389 
00390   // erm - what is this?!?! How could the HO probability be less than 0????? Do we want this line!?!?
00391   if(psz <= 0.) psz = 1.;
00392   if(pbz <= 0.) pbz = 1.;
00393 
00394   // Protection agains empty bins - set cal part to neutral if the bin of the template is empty 
00395   // (temporary fix, a proper extrapolation would be better)
00396   if (psx == 0. || pbx == 0.) {
00397     psx = 1.;
00398     pbx = 1.;
00399   } 
00400   if (psy == 0. || pby == 0.) {
00401     psy = 1.;
00402     pby = 1.;
00403   } 
00404   if (psz == 0. || pbz == 0.) {
00405     psz = 1.;
00406     pbz = 1.;
00407   }
00408 
00409   // There are two classes of events that deliver 0 for the hcal or ho energy:
00410   // 1) the track momentum is so low that the extrapolation tells us it should not have reached the cal
00411   // 2) the crossed cell had an reading below the readout cuts.
00412   // The 2nd case discriminates between muons and hadrons, the 1st not. Thus for the time being, 
00413   // we set the returned ps and pb to 0 for these cases.
00414   // We need to have a return value different from 0 for the 1st case in the long run.
00415   if ( had == 0.0 ) {
00416     psy = 1.;
00417     pby = 1.;
00418   } 
00419   if ( ho == 0.0 ) {
00420     psz = 1.;
00421     pbz = 1.;
00422   }
00423   
00424 
00425   // Set em to neutral if no energy in em or negative energy measured. 
00426   // (These cases might indicate problems in the ecal association or readout?! The only 
00427   // hint so far: for critical eta region (eta in [1.52, 1.64]) have negative em values.)
00428   if( em <= 0. && !use_em_special ) {
00429     pbx = 1.;
00430     psx = 1.;
00431   }
00432 
00433   if( (psx*psy*psz+pbx*pby*pbz) > 0. ) muon_compatibility = psx*psy*psz / (psx*psy*psz+pbx*pby*pbz);
00434   else {
00435     LogTrace("MuonIdentification")<<"Divide by 0 - defaulting consistency to 0.5 (neutral)!!";
00436     muon_compatibility = 0.5;
00437     LogTrace("MuonIdentification")<<"Input variables: eta    p     em     had    ho "<<"\n"
00438              <<eta<<" "<<p<<" "<<em<<" "<<had<<" "<<ho<<" "<<"\n"
00439              <<"cal uncorr:    em     had    ho "<<"\n"
00440              <<eta<<" "<<p<<" "<<amuon.calEnergy().em<<" "<<amuon.calEnergy().had<<" "<<amuon.calEnergy().ho;
00441   }
00442   return muon_compatibility;
00443 }