Reference Manual

Go to the documentation of this file.
00001 /*
00002  * \file EcalPreshowerSimHitsValidation.cc
00003  *
00004  * \author C.Rovelli
00005  *
00006 */
00007 
00008 #include <DataFormats/EcalDetId/interface/EEDetId.h>
00009 #include <DataFormats/EcalDetId/interface/ESDetId.h>
00010 #include "FWCore/Utilities/interface/Exception.h"
00011 #include "Validation/EcalHits/interface/EcalPreshowerSimHitsValidation.h"
00012 #include "DQMServices/Core/interface/DQMStore.h"
00013 
00014 using namespace cms;
00015 using namespace edm;
00016 using namespace std;
00017 
00018 EcalPreshowerSimHitsValidation::EcalPreshowerSimHitsValidation(const edm::ParameterSet& ps):
00019 
00020   HepMCLabel(ps.getParameter<std::string>("moduleLabelMC")),
00021   g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
00022   EEHitsCollection(ps.getParameter<std::string>("EEHitsCollection")),
00023   ESHitsCollection(ps.getParameter<std::string>("ESHitsCollection")){
00024 
00025   // verbosity switch
00026   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
00027  
00028   // get hold of back-end interface
00029   dbe_ = 0;
00030   dbe_ = edm::Service<DQMStore>().operator->();
00031   if ( dbe_ ) {
00032     if ( verbose_ ) { dbe_->setVerbose(1); } 
00033     else            { dbe_->setVerbose(0); }
00034   }
00035 
00036   if ( dbe_ ) {
00037     if ( verbose_ ) dbe_->showDirStructure();
00038   }
00039 
00040 
00041   menESHits1zp_ = 0;     
00042   menESHits2zp_ = 0;     
00043   menESHits1zm_ = 0;     
00044   menESHits2zm_ = 0;     
00045                                     
00046   meESEnergyHits1zp_ = 0;
00047   meESEnergyHits2zp_ = 0;
00048   meESEnergyHits1zm_ = 0;
00049   meESEnergyHits2zm_ = 0;
00050 
00051   meEShitLog10Energy_       = 0;
00052   meEShitLog10EnergyNorm_   = 0;
00053 
00054   meE1alphaE2zp_ = 0;
00055   meE1alphaE2zm_ = 0;
00056   meEEoverESzp_  = 0;
00057   meEEoverESzm_  = 0;
00058 
00059   me2eszpOver1eszp_ = 0;
00060   me2eszmOver1eszm_ = 0;
00061 
00062 
00063   Char_t histo[200];
00064  
00065   if ( dbe_ ) {
00066     dbe_->setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
00067   
00068     sprintf (histo, "ES hits layer 1 multiplicity z+" ) ;
00069     menESHits1zp_ = dbe_->book1D(histo, histo, 50, 0., 50. ) ;
00070 
00071     sprintf (histo, "ES hits layer 2 multiplicity z+" ) ;
00072     menESHits2zp_ = dbe_->book1D(histo, histo, 50, 0., 50. ) ;
00073 
00074     sprintf (histo, "ES hits layer 1 multiplicity z-" ) ;
00075     menESHits1zm_ = dbe_->book1D(histo, histo, 50, 0., 50. ) ;
00076 
00077     sprintf (histo, "ES hits layer 2 multiplicity z-" ) ;
00078     menESHits2zm_ = dbe_->book1D(histo, histo, 50, 0., 50. ) ;
00079 
00080     sprintf (histo, "ES hits energy layer 1 z+" ) ;
00081     meESEnergyHits1zp_ = dbe_->book1D(histo, histo, 100, 0., 0.05 ) ;
00082 
00083     sprintf (histo, "ES hits energy layer 2 z+" ) ;
00084     meESEnergyHits2zp_ = dbe_->book1D(histo, histo, 100, 0., 0.05 ) ;
00085 
00086     sprintf (histo, "ES hits energy layer 1 z-" ) ;
00087     meESEnergyHits1zm_ = dbe_->book1D(histo, histo, 100, 0., 0.05 ) ;
00088 
00089     sprintf (histo, "ES hits energy layer 2 z-" ) ;
00090     meESEnergyHits2zm_ = dbe_->book1D(histo, histo, 100, 0., 0.05 ) ;
00091 
00092     sprintf (histo, "ES hits log10energy spectrum" );
00093     meEShitLog10Energy_ = dbe_->book1D(histo, histo, 140, -10., 4.);
00094 
00095     sprintf (histo, "ES hits log10energy spectrum vs normalized energy" );
00096     meEShitLog10EnergyNorm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
00097 
00098     sprintf (histo, "ES E1+07E2 z+" ) ;
00099     meE1alphaE2zp_ = dbe_->book1D(histo, histo, 100, 0., 0.05);
00100 
00101     sprintf (histo, "ES E1+07E2 z-" ) ;
00102     meE1alphaE2zm_ = dbe_->book1D(histo, histo, 100, 0., 0.05);
00103 
00104     sprintf (histo, "EE vs ES z+" ) ;
00105     meEEoverESzp_ = dbe_->bookProfile(histo, histo, 250, 0., 500., 200, 0., 200.);
00106 
00107     sprintf (histo, "EE vs ES z-" ) ;
00108     meEEoverESzm_ = dbe_->bookProfile(histo, histo, 250, 0., 500., 200, 0., 200.);
00109 
00110     sprintf (histo, "ES ene2oEne1 z+" ) ;
00111     me2eszpOver1eszp_ = dbe_->book1D(histo, histo, 50, 0., 10.);
00112 
00113     sprintf (histo, "ES ene2oEne1 z-" ) ;
00114     me2eszmOver1eszm_ = dbe_->book1D(histo, histo, 50, 0., 10.);
00115   }
00116 }
00117 
00118 EcalPreshowerSimHitsValidation::~EcalPreshowerSimHitsValidation(){
00119 
00120 }
00121 
00122 void EcalPreshowerSimHitsValidation::beginJob(){
00123 
00124 }
00125 
00126 void EcalPreshowerSimHitsValidation::endJob(){
00127 
00128 }
00129 
00130 void EcalPreshowerSimHitsValidation::analyze(const edm::Event& e, const edm::EventSetup& c){
00131 
00132   edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
00133   
00134   edm::Handle<edm::HepMCProduct> MCEvt;
00135   e.getByLabel(HepMCLabel, MCEvt);
00136 
00137   edm::Handle<edm::PCaloHitContainer> EcalHitsEE;
00138   e.getByLabel(g4InfoLabel,EEHitsCollection,EcalHitsEE);
00139 
00140   edm::Handle<edm::PCaloHitContainer> EcalHitsES;
00141   e.getByLabel(g4InfoLabel,ESHitsCollection,EcalHitsES);
00142 
00143   std::vector<PCaloHit> theEECaloHits;  
00144   if( EcalHitsEE.isValid() ) {
00145     theEECaloHits.insert(theEECaloHits.end(), EcalHitsEE->begin(), EcalHitsEE->end());
00146   }
00147 
00148   std::vector<PCaloHit> theESCaloHits;
00149   if( EcalHitsES.isValid() ) {
00150     theESCaloHits.insert(theESCaloHits.end(), EcalHitsES->begin(), EcalHitsES->end());
00151   }
00152 
00153   double ESEnergy_ = 0.;
00154   //std::map<unsigned int, std::vector<PCaloHit>,std::less<unsigned int> > CaloHitMap;
00155 
00156   // endcap
00157   double EEetzp_ = 0.;
00158   double EEetzm_ = 0.;
00159   for (std::vector<PCaloHit>::iterator isim = theEECaloHits.begin(); isim != theEECaloHits.end(); ++isim){
00160     EEDetId eeid (isim->id()) ;
00161     if (eeid.zside() > 0 ) EEetzp_ += isim->energy();
00162     if (eeid.zside() < 0 ) EEetzm_ += isim->energy();
00163   }
00164   
00165   
00166   uint32_t nESHits1zp = 0;
00167   uint32_t nESHits1zm = 0;
00168   uint32_t nESHits2zp = 0;
00169   uint32_t nESHits2zm = 0;
00170   double ESet1zp_ = 0.;
00171   double ESet2zp_ = 0.;
00172   double ESet1zm_ = 0.;
00173   double ESet2zm_ = 0.;
00174   std::vector<double> econtr(140, 0. );
00175   
00176   for (std::vector<PCaloHit>::iterator isim = theESCaloHits.begin();
00177        isim != theESCaloHits.end(); ++isim){
00178     //CaloHitMap[(*isim).id()].push_back((*isim));
00179     
00180     ESDetId esid (isim->id()) ;
00181     
00182     LogDebug("HitInfo")
00183       << " CaloHit " << isim->getName() << "\n" 
00184       << " DetID = "<<isim->id()<< " ESDetId: z side " << esid.zside() << "  plane " << esid.plane() << esid.six() << ',' << esid.siy() << ':' << esid.strip() << "\n"
00185       << " Time = " << isim->time() << "\n"
00186       << " Track Id = " << isim->geantTrackId() << "\n"
00187       << " Energy = " << isim->energy();
00188     
00189     ESEnergy_ += isim->energy();
00190     if( isim->energy() > 0 ) {
00191       meEShitLog10Energy_->Fill(log10(isim->energy()));
00192       int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
00193       if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
00194     }
00195 
00196     if (esid.plane() == 1 ) { 
00197       if (esid.zside() > 0 ) {
00198         nESHits1zp++ ;
00199         ESet1zp_ += isim->energy();
00200         if (meESEnergyHits1zp_) meESEnergyHits1zp_->Fill(isim->energy()) ; 
00201       }
00202       else if (esid.zside() < 0 ) {
00203         nESHits1zm++ ; 
00204         ESet1zm_ += isim->energy();
00205         if (meESEnergyHits1zm_) meESEnergyHits1zm_->Fill(isim->energy()) ; 
00206       }
00207     }
00208     else if (esid.plane() == 2 ) {
00209       if (esid.zside() > 0 ) {
00210         nESHits2zp++ ; 
00211         ESet2zp_ += isim->energy();
00212         if (meESEnergyHits2zp_) meESEnergyHits2zp_->Fill(isim->energy()) ; 
00213       }
00214       else if (esid.zside() < 0 ) {
00215         nESHits2zm++ ; 
00216         ESet2zm_ += isim->energy();
00217         if (meESEnergyHits2zm_) meESEnergyHits2zm_->Fill(isim->energy()) ; 
00218       }
00219     }
00220     
00221   }
00222   
00223   if (menESHits1zp_) menESHits1zp_->Fill(nESHits1zp);
00224   if (menESHits1zm_) menESHits1zm_->Fill(nESHits1zm);
00225   
00226   if (menESHits2zp_) menESHits2zp_->Fill(nESHits2zp);
00227   if (menESHits2zm_) menESHits2zm_->Fill(nESHits2zm);
00228   
00229   if( meEShitLog10EnergyNorm_ && ESEnergy_ != 0 ) {
00230     for( int i=0; i<140; i++ ) {
00231       meEShitLog10EnergyNorm_->Fill( -10.+(float(i)+0.5)/10., econtr[i]/ESEnergy_ );
00232     }
00233   }
00234 
00235 
00236   for ( HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin();
00237         p != MCEvt->GetEvent()->particles_end(); ++p ) {
00238     
00239     double htheta = (*p)->momentum().theta();
00240     double heta = -99999.;
00241     if( tan(htheta * 0.5) > 0 ) {
00242       heta = -log(tan(htheta * 0.5));
00243     }
00244 
00245     if ( heta > 1.653 && heta < 2.6 ) {
00246 
00247       if (meE1alphaE2zp_) meE1alphaE2zp_->Fill(ESet1zp_+0.7*ESet2zp_);
00248       if (meEEoverESzp_)  meEEoverESzp_ ->Fill((ESet1zp_+0.7*ESet2zp_)/0.00009, EEetzp_);
00249       if ((me2eszpOver1eszp_) && (ESet1zp_ != 0.)) me2eszpOver1eszp_->Fill(ESet2zp_/ESet1zp_);
00250     }
00251     if ( heta < -1.653 && heta > -2.6 ) {
00252       if (meE1alphaE2zm_) meE1alphaE2zm_->Fill(ESet1zm_+0.7*ESet2zm_);
00253       if (meEEoverESzm_)  meEEoverESzm_ ->Fill((ESet1zm_+0.7*ESet2zm_)/0.00009, EEetzm_);
00254       if ((me2eszmOver1eszm_) && (ESet1zm_ != 0.)) me2eszmOver1eszm_->Fill(ESet2zm_/ESet1zm_);
00255     }
00256   }
00257 }
00258 
00259 
00260 //  LocalWords:  EcalSimHitsValidation
/data/refman/pasoursint/CMSSW_4_2_9_HLT1_bphpatch4/src/Validation/EcalHits/src/EcalPreshowerSimHitsValidation.cc
