CMSSW Reference Manual

00001 /*
00002  * \file EcalBarrelSimHitsValidation.cc
00003  *
00004  * \author C.Rovelli
00005  *
00006 */
00007 
00008 #include "FWCore/Utilities/interface/Exception.h"
00009 #include <DataFormats/EcalDetId/interface/EBDetId.h>
00010 #include "Validation/EcalHits/interface/EcalBarrelSimHitsValidation.h"
00011 #include "DQMServices/Core/interface/DQMStore.h"
00012 
00013 using namespace cms;
00014 using namespace edm;
00015 using namespace std;
00016 
00017 EcalBarrelSimHitsValidation::EcalBarrelSimHitsValidation(const edm::ParameterSet& ps):
00018   g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
00019   EBHitsCollection(ps.getParameter<std::string>("EBHitsCollection")),
00020   ValidationCollection(ps.getParameter<std::string>("ValidationCollection")){   
00021 
00022   // verbosity switch
00023   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
00024  
00025   // get hold of back-end interface
00026   dbe_ = 0;
00027   dbe_ = edm::Service<DQMStore>().operator->();           
00028   if ( dbe_ ) {
00029     if ( verbose_ ) { dbe_->setVerbose(1); } 
00030     else            { dbe_->setVerbose(0); }
00031   }
00032                                                                                
00033   if ( dbe_ ) {
00034     if ( verbose_ ) dbe_->showDirStructure();
00035   }
00036  
00037   menEBHits_       = 0;
00038   menEBCrystals_   = 0;
00039   meEBOccupancy_   = 0;
00040   meEBLongitudinalShower_ = 0;
00041   meEBhitEnergy_   = 0;
00042   meEBhitLog10Energy_       = 0;
00043   meEBhitLog10EnergyNorm_   = 0;
00044   meEBhitLog10Energy25Norm_ = 0;
00045   meEBhitEnergy2_   = 0;
00046   meEBcrystalEnergy_   = 0;
00047   meEBcrystalEnergy2_   = 0;
00048 
00049   meEBe1_  = 0;  
00050   meEBe4_  = 0;  
00051   meEBe9_  = 0;  
00052   meEBe16_ = 0; 
00053   meEBe25_ = 0; 
00054 
00055   meEBe1oe4_   = 0;  
00056   meEBe1oe9_   = 0;  
00057   meEBe4oe9_   = 0;  
00058   meEBe9oe16_  = 0;
00059   meEBe1oe25_  = 0;  
00060   meEBe9oe25_  = 0; 
00061   meEBe16oe25_ = 0;
00062 
00063   myEntries = 0;
00064   for (int myStep = 0; myStep<26; myStep++) { eRLength[myStep] = 0.0; }
00065 
00066   Char_t histo[200];
00067    
00068   if ( dbe_ ) {
00069     dbe_->setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
00070   
00071     sprintf (histo, "EB hits multiplicity" ) ;
00072     menEBHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ;
00073  
00074     sprintf (histo, "EB crystals multiplicity" ) ;
00075     menEBCrystals_ = dbe_->book1D(histo, histo, 200, 0., 2000.) ; 
00076 
00077     sprintf (histo, "EB occupancy" ) ;
00078     meEBOccupancy_ = dbe_->book2D(histo, histo, 360, 0., 360., 170, -85., 85.);
00079   
00080     sprintf (histo, "EB longitudinal shower profile" ) ;
00081     meEBLongitudinalShower_ = dbe_->bookProfile(histo, histo, 26, 0., 26., 100, 0., 20000.);
00082 
00083     sprintf (histo, "EB hits energy spectrum" );
00084     meEBhitEnergy_ = dbe_->book1D(histo, histo, 4000, 0., 400.);
00085 
00086     sprintf (histo, "EB hits log10energy spectrum" );
00087     meEBhitLog10Energy_ = dbe_->book1D(histo, histo, 140, -10., 4.);
00088 
00089     sprintf (histo, "EB hits log10energy spectrum vs normalized energy" );
00090     meEBhitLog10EnergyNorm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
00091 
00092     sprintf (histo, "EB hits log10energy spectrum vs normalized energy25" );
00093     meEBhitLog10Energy25Norm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
00094 
00095     sprintf (histo, "EB hits energy spectrum 2" );
00096     meEBhitEnergy2_ = dbe_->book1D(histo, histo, 1000, 0., 0.001);
00097 
00098     sprintf (histo, "EB crystal energy spectrum" );
00099     meEBcrystalEnergy_ = dbe_->book1D(histo, histo, 5000, 0., 50.);
00100 
00101     sprintf (histo, "EB crystal energy spectrum 2" );
00102     meEBcrystalEnergy2_ = dbe_->book1D(histo, histo, 1000, 0., 0.001);
00103 
00104     sprintf (histo, "EB E1" ) ;
00105     meEBe1_ = dbe_->book1D(histo, histo, 400, 0., 400.);
00106 
00107     sprintf (histo, "EB E4" ) ;
00108     meEBe4_ = dbe_->book1D(histo, histo, 400, 0., 400.);
00109 
00110     sprintf (histo, "EB E9" ) ;
00111     meEBe9_ = dbe_->book1D(histo, histo, 400, 0., 400.);
00112 
00113     sprintf (histo, "EB E16" ) ;
00114     meEBe16_ = dbe_->book1D(histo, histo, 400, 0., 400.);
00115 
00116     sprintf (histo, "EB E25" ) ;
00117     meEBe25_ = dbe_->book1D(histo, histo, 400, 0., 400.);
00118 
00119     sprintf (histo, "EB E1oE4" ) ;
00120     meEBe1oe4_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00121 
00122     sprintf (histo, "EB E1oE9" ) ;
00123     meEBe1oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00124 
00125     sprintf (histo, "EB E4oE9" ) ;
00126     meEBe4oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00127 
00128     sprintf (histo, "EB E9oE16" ) ;
00129     meEBe9oe16_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00130 
00131     sprintf (histo, "EB E1oE25" ) ;
00132     meEBe1oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00133 
00134     sprintf (histo, "EB E9oE25" ) ;
00135     meEBe9oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00136 
00137     sprintf (histo, "EB E16oE25" ) ;
00138     meEBe16oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00139   }
00140 }
00141 
00142 EcalBarrelSimHitsValidation::~EcalBarrelSimHitsValidation(){
00143  
00144 }
00145 
00146 void EcalBarrelSimHitsValidation::beginJob(const edm::EventSetup& c){
00147 
00148 }
00149 
00150 void EcalBarrelSimHitsValidation::endJob(){
00151   
00152   for (int myStep=0; myStep<26; myStep++){
00153     if (meEBLongitudinalShower_) meEBLongitudinalShower_->Fill(float(myStep), eRLength[myStep]/myEntries);
00154   }
00155 }
00156 
00157 void EcalBarrelSimHitsValidation::analyze(const edm::Event& e, const edm::EventSetup& c){
00158 
00159   edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
00160   
00161   edm::Handle<edm::PCaloHitContainer> EcalHitsEB;
00162   e.getByLabel(g4InfoLabel,EBHitsCollection,EcalHitsEB);
00163 
00164   // Do nothing if no Barrel data available
00165   if( ! EcalHitsEB.isValid() ) return;
00166 
00167   edm::Handle<PEcalValidInfo> MyPEcalValidInfo;
00168   e.getByLabel(g4InfoLabel,ValidationCollection,MyPEcalValidInfo);
00169 
00170   std::vector<PCaloHit> theEBCaloHits;
00171   theEBCaloHits.insert(theEBCaloHits.end(), EcalHitsEB->begin(), EcalHitsEB->end());
00172 
00173   myEntries++;
00174 
00175   double EBEnergy_ = 0.;
00176   std::map<unsigned int, std::vector<PCaloHit*>,std::less<unsigned int> > CaloHitMap;
00177   
00178   double eb1  = 0.0;
00179   double eb4  = 0.0;
00180   double eb9  = 0.0;
00181   double eb16 = 0.0;
00182   double eb25 = 0.0;
00183   std::vector<double> econtr(140, 0. );
00184   std::vector<double> econtr25(140, 0. );
00185   
00186   MapType ebmap;
00187   uint32_t nEBHits = 0;
00188   
00189   for (std::vector<PCaloHit>::iterator isim = theEBCaloHits.begin();
00190        isim != theEBCaloHits.end(); ++isim){
00191 
00192     if ( isim->time() > 500. ) { continue; }
00193 
00194     CaloHitMap[ isim->id()].push_back( &(*isim) );
00195     
00196     EBDetId ebid (isim->id()) ;
00197     
00198     LogDebug("HitInfo") 
00199       << " CaloHit "    << isim->getName() << "\n" 
00200       << " DetID = "    << isim->id()   << " EBDetId = " << ebid.ieta() << " " << ebid.iphi() << "\n" 
00201       << " Time = "     << isim->time()    << "\n"
00202       << " Track Id = " << isim->geantTrackId() << "\n"
00203       << " Energy = "   << isim->energy();
00204     
00205     if (meEBOccupancy_) meEBOccupancy_->Fill( ebid.iphi(), ebid.ieta() );
00206     
00207     uint32_t crystid = ebid.rawId();
00208     ebmap[crystid] += isim->energy();
00209     
00210     EBEnergy_ += isim->energy();
00211     nEBHits++;
00212     meEBhitEnergy_->Fill(isim->energy());
00213     meEBhitLog10Energy_->Fill(log10(isim->energy()));
00214     meEBhitEnergy2_->Fill(isim->energy());
00215     
00216     int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
00217     if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
00218 
00219   }
00220 
00221   if (menEBCrystals_) menEBCrystals_->Fill(ebmap.size());
00222   if (meEBcrystalEnergy_) {
00223     for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = ebmap.begin(); it != ebmap.end(); ++it ) meEBcrystalEnergy_->Fill((*it).second);
00224   }
00225   if (meEBcrystalEnergy2_) {
00226     for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = ebmap.begin(); it != ebmap.end(); ++it ) meEBcrystalEnergy2_->Fill((*it).second);
00227   }
00228   
00229   if (menEBHits_) menEBHits_->Fill(nEBHits);
00230   
00231   if (nEBHits > 0 ) {
00232     
00233     uint32_t  ebcenterid = getUnitWithMaxEnergy(ebmap);
00234     EBDetId myEBid(ebcenterid);
00235     int bx = myEBid.ietaAbs();
00236     int by = myEBid.iphi();
00237     int bz = myEBid.zside();
00238     eb1 =  energyInMatrixEB(1,1,bx,by,bz,ebmap);
00239     if (meEBe1_) meEBe1_->Fill(eb1);
00240     eb9 =  energyInMatrixEB(3,3,bx,by,bz,ebmap);
00241     if (meEBe9_) meEBe9_->Fill(eb9);
00242     eb25=  energyInMatrixEB(5,5,bx,by,bz,ebmap);
00243     if (meEBe25_) meEBe25_->Fill(eb25);
00244     
00245     std::vector<uint32_t> ids25; ids25 = getIdsAroundMax(5,5,bx,by,bz,ebmap);
00246 
00247     for( unsigned i=0; i<25; i++ ) {
00248       for( unsigned int j=0; j<CaloHitMap[ids25[i]].size(); j++ ) {
00249         int log10i = int( ( log10( CaloHitMap[ids25[i]][j]->energy()) + 10. ) * 10. );
00250         if( log10i >=0 && log10i < 140 ) econtr25[log10i] += CaloHitMap[ids25[i]][j]->energy();
00251       }
00252     }
00253 
00254     MapType  newebmap;
00255     if( fillEBMatrix(3,3,bx,by,bz,newebmap, ebmap)){
00256       eb4 = eCluster2x2(newebmap);
00257       if (meEBe4_) meEBe4_->Fill(eb4);
00258     }
00259     if( fillEBMatrix(5,5,bx,by,bz,newebmap, ebmap)){
00260       eb16 = eCluster4x4(eb9,newebmap); 
00261       if (meEBe16_) meEBe16_->Fill(eb16);
00262     }
00263     
00264     if (meEBe1oe4_   && eb4  > 0.1 ) meEBe1oe4_  -> Fill(eb1/eb4);
00265     if (meEBe1oe9_   && eb9  > 0.1 ) meEBe1oe9_  -> Fill(eb1/eb9);
00266     if (meEBe4oe9_   && eb9  > 0.1 ) meEBe4oe9_  -> Fill(eb4/eb9);
00267     if (meEBe9oe16_  && eb16 > 0.1 ) meEBe9oe16_ -> Fill(eb9/eb16);
00268     if (meEBe1oe25_  && eb25 > 0.1 ) meEBe1oe25_ -> Fill(eb1/eb25);
00269     if (meEBe9oe25_  && eb25 > 0.1 ) meEBe9oe25_ -> Fill(eb9/eb25);
00270     if (meEBe16oe25_ && eb25 > 0.1 ) meEBe16oe25_-> Fill(eb16/eb25);
00271 
00272     if( meEBhitLog10EnergyNorm_ && EBEnergy_ != 0 ) {
00273       for( int i=0; i<140; i++ ) {
00274         meEBhitLog10EnergyNorm_->Fill( -10.+(float(i)+0.5)/10., econtr[i]/EBEnergy_ );
00275       }
00276     }
00277 
00278     if( meEBhitLog10Energy25Norm_ && eb25 != 0 ) {
00279       for( int i=0; i<140; i++ ) {
00280         meEBhitLog10Energy25Norm_->Fill( -10.+(float(i)+0.5)/10., econtr25[i]/eb25 );
00281       }
00282     }
00283 
00284 
00285   }
00286   
00287   if( MyPEcalValidInfo.isValid() ) {
00288     if ( MyPEcalValidInfo->eb1x1() > 0. ) {
00289       std::vector<float>  BX0 = MyPEcalValidInfo->bX0();
00290       for (int myStep=0; myStep< 26; myStep++ ) { eRLength[myStep] += BX0[myStep]; }
00291     }
00292   }
00293 }
00294 
00295 float EcalBarrelSimHitsValidation::energyInMatrixEB(int nCellInEta, int nCellInPhi,
00296                                               int centralEta, int centralPhi,
00297                                               int centralZ, MapType& themap){
00298 
00299   int   ncristals   = 0;
00300   float totalEnergy = 0.;
00301   
00302   int goBackInEta = nCellInEta/2;
00303   int goBackInPhi = nCellInPhi/2;
00304   int startEta    = centralZ*centralEta-goBackInEta;
00305   int startPhi    = centralPhi-goBackInPhi;
00306   
00307   for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
00308     for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
00309       
00310       uint32_t index ;
00311       if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00312       if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
00313       else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00314       else              { index = EBDetId(ieta,iphi).rawId();     }
00315 
00316       totalEnergy   += themap[index];
00317       ncristals     += 1;
00318     }
00319   }
00320   
00321   LogDebug("GeomInfo")
00322     << nCellInEta << " x " << nCellInPhi 
00323     << " EB matrix energy = " << totalEnergy
00324     << " for " << ncristals << " crystals" ;
00325   return totalEnergy;
00326   
00327 }   
00328 
00329 std::vector<uint32_t> EcalBarrelSimHitsValidation::getIdsAroundMax( int nCellInEta, int nCellInPhi, int centralEta, int centralPhi, int centralZ, MapType& themap){
00330 
00331   int   ncristals   = 0;
00332   std::vector<uint32_t> ids( nCellInEta*nCellInPhi );
00333   
00334   int goBackInEta = nCellInEta/2;
00335   int goBackInPhi = nCellInPhi/2;
00336   int startEta    = centralZ*centralEta-goBackInEta;
00337   int startPhi    = centralPhi-goBackInPhi;
00338   
00339   for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
00340     for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
00341       
00342       uint32_t index ;
00343       if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00344       if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
00345       else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00346       else              { index = EBDetId(ieta,iphi).rawId();     }
00347       ids[ncristals] = index;
00348       ncristals     += 1;
00349     }
00350   }
00351   
00352   return ids;
00353   
00354 }   
00355 
00356 bool  EcalBarrelSimHitsValidation::fillEBMatrix(int nCellInEta, int nCellInPhi,
00357                                     int CentralEta, int CentralPhi,int CentralZ,
00358                                     MapType& fillmap, MapType&  themap)
00359 {
00360   int goBackInEta = nCellInEta/2;
00361   int goBackInPhi = nCellInPhi/2;
00362   
00363   int startEta  =  CentralZ*CentralEta - goBackInEta;
00364   int startPhi  =  CentralPhi - goBackInPhi;
00365   
00366   int i = 0 ;
00367   for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
00368     for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
00369       uint32_t  index;
00370       if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00371       if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
00372       else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00373       else              { index = EBDetId(ieta,iphi).rawId();     }
00374       fillmap[i++] = themap[index];
00375     }
00376   }
00377   
00378   uint32_t  ebcenterid = getUnitWithMaxEnergy(themap);
00379   
00380   if ( fillmap[i/2] == themap[ebcenterid] ) 
00381     return true;
00382   else 
00383     return false;
00384 } 
00385 
00386 float EcalBarrelSimHitsValidation::eCluster2x2( MapType& themap){
00387   float  E22=0.;
00388   float e012  = themap[0]+themap[1]+themap[2];
00389   float e036  = themap[0]+themap[3]+themap[6];
00390   float e678  = themap[6]+themap[7]+themap[8];
00391   float e258  = themap[2]+themap[5]+themap[8];
00392 
00393   if ( (e012>e678 || e012==e678) && (e036>e258  || e036==e258))
00394     return     E22=themap[0]+themap[1]+themap[3]+themap[4];
00395   else if ( (e012>e678 || e012==e678)  && (e036<e258 || e036==e258) )
00396     return    E22=themap[1]+themap[2]+themap[4]+themap[5];
00397   else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
00398     return     E22=themap[3]+themap[4]+themap[6]+themap[7];
00399   else if ( (e012<e678|| e012==e678)  && (e036<e258|| e036==e258) )
00400     return    E22=themap[4]+themap[5]+themap[7]+themap[8];
00401   else {
00402     return E22;
00403   }
00404 }
00405 
00406 float EcalBarrelSimHitsValidation::eCluster4x4(float e33,  MapType&  themap){
00407   float E44=0.;
00408   float e0_4   = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
00409   float e0_20  = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
00410   float e4_24  = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
00411   float e0_24  = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];
00412   
00413   if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
00414     return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
00415   else if ((e0_4>e0_24 || e0_4==e0_24)  && (e0_20<e4_24 || e0_20==e4_24))
00416     return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
00417   else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
00418     return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
00419   else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
00420     return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
00421   else{
00422     return E44;
00423   }
00424 }
00425 
00426 uint32_t EcalBarrelSimHitsValidation::getUnitWithMaxEnergy(MapType& themap) {
00427   
00428   //look for max
00429   uint32_t unitWithMaxEnergy = 0;
00430   float    maxEnergy = 0.;
00431   
00432   MapType::iterator iter;
00433   for (iter = themap.begin(); iter != themap.end(); iter++) {
00434     
00435     if (maxEnergy < (*iter).second) {
00436       maxEnergy = (*iter).second;       
00437       unitWithMaxEnergy = (*iter).first;
00438     }                           
00439   }
00440   
00441   LogDebug("GeomInfo")
00442     << " max energy of " << maxEnergy 
00443     << " GeV in Unit id " << unitWithMaxEnergy;
00444   return unitWithMaxEnergy;
00445 }
00446 
00447