Reference Manual

Go to the documentation of this file.
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(){
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     if( isim->energy() > 0 ) {
00214       meEBhitLog10Energy_->Fill(log10(isim->energy()));
00215       int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
00216       if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
00217     }
00218     meEBhitEnergy2_->Fill(isim->energy());
00219     
00220   }
00221 
00222   if (menEBCrystals_) menEBCrystals_->Fill(ebmap.size());
00223   if (meEBcrystalEnergy_) {
00224     for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = ebmap.begin(); it != ebmap.end(); ++it ) meEBcrystalEnergy_->Fill((*it).second);
00225   }
00226   if (meEBcrystalEnergy2_) {
00227     for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = ebmap.begin(); it != ebmap.end(); ++it ) meEBcrystalEnergy2_->Fill((*it).second);
00228   }
00229   
00230   if (menEBHits_) menEBHits_->Fill(nEBHits);
00231   
00232   if (nEBHits > 0 ) {
00233     
00234     uint32_t  ebcenterid = getUnitWithMaxEnergy(ebmap);
00235     EBDetId myEBid(ebcenterid);
00236     int bx = myEBid.ietaAbs();
00237     int by = myEBid.iphi();
00238     int bz = myEBid.zside();
00239     eb1 =  energyInMatrixEB(1,1,bx,by,bz,ebmap);
00240     if (meEBe1_) meEBe1_->Fill(eb1);
00241     eb9 =  energyInMatrixEB(3,3,bx,by,bz,ebmap);
00242     if (meEBe9_) meEBe9_->Fill(eb9);
00243     eb25=  energyInMatrixEB(5,5,bx,by,bz,ebmap);
00244     if (meEBe25_) meEBe25_->Fill(eb25);
00245     
00246     std::vector<uint32_t> ids25; ids25 = getIdsAroundMax(5,5,bx,by,bz,ebmap);
00247 
00248     for( unsigned i=0; i<25; i++ ) {
00249       for( unsigned int j=0; j<CaloHitMap[ids25[i]].size(); j++ ) {
00250         if( CaloHitMap[ids25[i]][j]->energy() > 0 ) {
00251           int log10i = int( ( log10( CaloHitMap[ids25[i]][j]->energy()) + 10. ) * 10. );
00252           if( log10i >=0 && log10i < 140 ) econtr25[log10i] += CaloHitMap[ids25[i]][j]->energy();
00253         }
00254       }
00255     }
00256 
00257     MapType  newebmap;
00258     if( fillEBMatrix(3,3,bx,by,bz,newebmap, ebmap)){
00259       eb4 = eCluster2x2(newebmap);
00260       if (meEBe4_) meEBe4_->Fill(eb4);
00261     }
00262     if( fillEBMatrix(5,5,bx,by,bz,newebmap, ebmap)){
00263       eb16 = eCluster4x4(eb9,newebmap); 
00264       if (meEBe16_) meEBe16_->Fill(eb16);
00265     }
00266     
00267     if (meEBe1oe4_   && eb4  > 0.1 ) meEBe1oe4_  -> Fill(eb1/eb4);
00268     if (meEBe1oe9_   && eb9  > 0.1 ) meEBe1oe9_  -> Fill(eb1/eb9);
00269     if (meEBe4oe9_   && eb9  > 0.1 ) meEBe4oe9_  -> Fill(eb4/eb9);
00270     if (meEBe9oe16_  && eb16 > 0.1 ) meEBe9oe16_ -> Fill(eb9/eb16);
00271     if (meEBe1oe25_  && eb25 > 0.1 ) meEBe1oe25_ -> Fill(eb1/eb25);
00272     if (meEBe9oe25_  && eb25 > 0.1 ) meEBe9oe25_ -> Fill(eb9/eb25);
00273     if (meEBe16oe25_ && eb25 > 0.1 ) meEBe16oe25_-> Fill(eb16/eb25);
00274 
00275     if( meEBhitLog10EnergyNorm_ && EBEnergy_ != 0 ) {
00276       for( int i=0; i<140; i++ ) {
00277         meEBhitLog10EnergyNorm_->Fill( -10.+(float(i)+0.5)/10., econtr[i]/EBEnergy_ );
00278       }
00279     }
00280 
00281     if( meEBhitLog10Energy25Norm_ && eb25 != 0 ) {
00282       for( int i=0; i<140; i++ ) {
00283         meEBhitLog10Energy25Norm_->Fill( -10.+(float(i)+0.5)/10., econtr25[i]/eb25 );
00284       }
00285     }
00286 
00287 
00288   }
00289   
00290 
00291   if( MyPEcalValidInfo.isValid() ) {
00292     if ( MyPEcalValidInfo->eb1x1() > 0. ) {
00293       std::vector<float>  BX0 = MyPEcalValidInfo->bX0();
00294       if (meEBLongitudinalShower_) meEBLongitudinalShower_->Reset();
00295       for (int myStep=0; myStep< 26; myStep++ ) { 
00296         eRLength[myStep] += BX0[myStep]; 
00297         if (meEBLongitudinalShower_) meEBLongitudinalShower_->Fill(float(myStep), eRLength[myStep]/myEntries);
00298       }
00299     }
00300   }
00301 }
00302 
00303 float EcalBarrelSimHitsValidation::energyInMatrixEB(int nCellInEta, int nCellInPhi,
00304                                               int centralEta, int centralPhi,
00305                                               int centralZ, MapType& themap){
00306 
00307   int   ncristals   = 0;
00308   float totalEnergy = 0.;
00309   
00310   int goBackInEta = nCellInEta/2;
00311   int goBackInPhi = nCellInPhi/2;
00312   int startEta    = centralZ*centralEta-goBackInEta;
00313   int startPhi    = centralPhi-goBackInPhi;
00314   
00315   for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
00316     for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
00317       
00318       uint32_t index ;
00319       if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00320       if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
00321       else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00322       else              { index = EBDetId(ieta,iphi).rawId();     }
00323 
00324       totalEnergy   += themap[index];
00325       ncristals     += 1;
00326     }
00327   }
00328   
00329   LogDebug("GeomInfo")
00330     << nCellInEta << " x " << nCellInPhi 
00331     << " EB matrix energy = " << totalEnergy
00332     << " for " << ncristals << " crystals" ;
00333   return totalEnergy;
00334   
00335 }   
00336 
00337 std::vector<uint32_t> EcalBarrelSimHitsValidation::getIdsAroundMax( int nCellInEta, int nCellInPhi, int centralEta, int centralPhi, int centralZ, MapType& themap){
00338 
00339   int   ncristals   = 0;
00340   std::vector<uint32_t> ids( nCellInEta*nCellInPhi );
00341   
00342   int goBackInEta = nCellInEta/2;
00343   int goBackInPhi = nCellInPhi/2;
00344   int startEta    = centralZ*centralEta-goBackInEta;
00345   int startPhi    = centralPhi-goBackInPhi;
00346   
00347   for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
00348     for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
00349       
00350       uint32_t index ;
00351       if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00352       if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
00353       else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00354       else              { index = EBDetId(ieta,iphi).rawId();     }
00355       ids[ncristals] = index;
00356       ncristals     += 1;
00357     }
00358   }
00359   
00360   return ids;
00361   
00362 }   
00363 
00364 bool  EcalBarrelSimHitsValidation::fillEBMatrix(int nCellInEta, int nCellInPhi,
00365                                     int CentralEta, int CentralPhi,int CentralZ,
00366                                     MapType& fillmap, MapType&  themap)
00367 {
00368   int goBackInEta = nCellInEta/2;
00369   int goBackInPhi = nCellInPhi/2;
00370   
00371   int startEta  =  CentralZ*CentralEta - goBackInEta;
00372   int startPhi  =  CentralPhi - goBackInPhi;
00373   
00374   int i = 0 ;
00375   for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
00376     for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
00377       uint32_t  index;
00378       if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00379       if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
00380       else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00381       else              { index = EBDetId(ieta,iphi).rawId();     }
00382       fillmap[i++] = themap[index];
00383     }
00384   }
00385   
00386   uint32_t  ebcenterid = getUnitWithMaxEnergy(themap);
00387   
00388   if ( fillmap[i/2] == themap[ebcenterid] ) 
00389     return true;
00390   else 
00391     return false;
00392 } 
00393 
00394 float EcalBarrelSimHitsValidation::eCluster2x2( MapType& themap){
00395   float  E22=0.;
00396   float e012  = themap[0]+themap[1]+themap[2];
00397   float e036  = themap[0]+themap[3]+themap[6];
00398   float e678  = themap[6]+themap[7]+themap[8];
00399   float e258  = themap[2]+themap[5]+themap[8];
00400 
00401   if ( (e012>e678 || e012==e678) && (e036>e258  || e036==e258))
00402     return     E22=themap[0]+themap[1]+themap[3]+themap[4];
00403   else if ( (e012>e678 || e012==e678)  && (e036<e258 || e036==e258) )
00404     return    E22=themap[1]+themap[2]+themap[4]+themap[5];
00405   else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
00406     return     E22=themap[3]+themap[4]+themap[6]+themap[7];
00407   else if ( (e012<e678|| e012==e678)  && (e036<e258|| e036==e258) )
00408     return    E22=themap[4]+themap[5]+themap[7]+themap[8];
00409   else {
00410     return E22;
00411   }
00412 }
00413 
00414 float EcalBarrelSimHitsValidation::eCluster4x4(float e33,  MapType&  themap){
00415   float E44=0.;
00416   float e0_4   = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
00417   float e0_20  = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
00418   float e4_24  = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
00419   float e0_24  = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];
00420   
00421   if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
00422     return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
00423   else if ((e0_4>e0_24 || e0_4==e0_24)  && (e0_20<e4_24 || e0_20==e4_24))
00424     return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
00425   else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
00426     return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
00427   else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
00428     return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
00429   else{
00430     return E44;
00431   }
00432 }
00433 
00434 uint32_t EcalBarrelSimHitsValidation::getUnitWithMaxEnergy(MapType& themap) {
00435   
00436   //look for max
00437   uint32_t unitWithMaxEnergy = 0;
00438   float    maxEnergy = 0.;
00439   
00440   MapType::iterator iter;
00441   for (iter = themap.begin(); iter != themap.end(); iter++) {
00442     
00443     if (maxEnergy < (*iter).second) {
00444       maxEnergy = (*iter).second;       
00445       unitWithMaxEnergy = (*iter).first;
00446     }                           
00447   }
00448   
00449   LogDebug("GeomInfo")
00450     << " max energy of " << maxEnergy 
00451     << " GeV in Unit id " << unitWithMaxEnergy;
00452   return unitWithMaxEnergy;
00453 }
00454 
00455
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/Validation/EcalHits/src/EcalBarrelSimHitsValidation.cc
