CMSSW Reference Manual

00001 /*
00002  * \file EcalEndcapSimHitsValidation.cc
00003  *
00004  * \author C.Rovelli
00005  *
00006 */
00007 
00008 #include <DataFormats/EcalDetId/interface/EEDetId.h>
00009 #include "FWCore/Utilities/interface/Exception.h"
00010 #include "Validation/EcalHits/interface/EcalEndcapSimHitsValidation.h"
00011 #include "DQMServices/Core/interface/DQMStore.h"
00012 
00013 using namespace cms;
00014 using namespace edm;
00015 using namespace std;
00016 
00017 EcalEndcapSimHitsValidation::EcalEndcapSimHitsValidation(const edm::ParameterSet& ps):
00018   g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
00019   EEHitsCollection(ps.getParameter<std::string>("EEHitsCollection")),
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 
00038   meEEzpHits_             = 0;
00039   meEEzmHits_             = 0;
00040   meEEzpCrystals_         = 0;
00041   meEEzmCrystals_         = 0;
00042   meEEzpOccupancy_        = 0;
00043   meEEzmOccupancy_        = 0;
00044   meEELongitudinalShower_ = 0;
00045   meEEHitEnergy_          = 0;
00046   meEEhitLog10Energy_       = 0;
00047   meEEhitLog10EnergyNorm_   = 0;
00048   meEEhitLog10Energy25Norm_ = 0;
00049   meEEHitEnergy2_         = 0;
00050   meEEcrystalEnergy_      = 0;
00051   meEEcrystalEnergy2_     = 0;
00052 
00053   meEEe1_  = 0;  
00054   meEEe4_  = 0;  
00055   meEEe9_  = 0;  
00056   meEEe16_ = 0; 
00057   meEEe25_ = 0; 
00058 
00059   meEEe1oe4_   = 0;  
00060   meEEe1oe9_   = 0;  
00061   meEEe4oe9_   = 0;  
00062   meEEe9oe16_  = 0;
00063   meEEe1oe25_  = 0;  
00064   meEEe9oe25_  = 0; 
00065   meEEe16oe25_ = 0;
00066 
00067   myEntries = 0;
00068   for ( int myStep = 0; myStep<26; myStep++) { eRLength[myStep] = 0.0; }
00069 
00070   Char_t histo[200];
00071  
00072   if ( dbe_ ) {
00073     dbe_->setCurrentFolder("EcalHitsV/EcalSimHitsValidation");
00074   
00075     sprintf (histo, "EE+ hits multiplicity" ) ;
00076     meEEzpHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ; 
00077 
00078     sprintf (histo, "EE- hits multiplicity" ) ;
00079     meEEzmHits_ = dbe_->book1D(histo, histo, 50, 0., 5000.) ; 
00080 
00081     sprintf (histo, "EE+ crystals multiplicity" ) ;
00082     meEEzpCrystals_ = dbe_->book1D(histo, histo, 200, 0., 2000.) ; 
00083 
00084     sprintf (histo, "EE- crystals multiplicity" ) ;
00085     meEEzmCrystals_ = dbe_->book1D(histo, histo, 200, 0., 2000.) ; 
00086 
00087     sprintf (histo, "EE+ occupancy" ) ;
00088     meEEzpOccupancy_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
00089   
00090     sprintf (histo, "EE- occupancy" ) ;
00091     meEEzmOccupancy_ = dbe_->book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
00092   
00093     sprintf (histo, "EE longitudinal shower profile" ) ;
00094     meEELongitudinalShower_ = dbe_->bookProfile(histo, histo, 26,0,26, 100, 0, 20000);
00095 
00096     sprintf (histo, "EE hits energy spectrum" );
00097     meEEHitEnergy_ = dbe_->book1D(histo, histo, 4000, 0., 400.);
00098 
00099     sprintf (histo, "EE hits log10energy spectrum" );
00100     meEEhitLog10Energy_ = dbe_->book1D(histo, histo, 140, -10., 4.);
00101 
00102     sprintf (histo, "EE hits log10energy spectrum vs normalized energy" );
00103     meEEhitLog10EnergyNorm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
00104 
00105     sprintf (histo, "EE hits log10energy spectrum vs normalized energy25" );
00106     meEEhitLog10Energy25Norm_ = dbe_->bookProfile(histo, histo, 140, -10., 4., 100, 0., 1.);
00107 
00108     sprintf (histo, "EE hits energy spectrum 2" );
00109     meEEHitEnergy2_ = dbe_->book1D(histo, histo, 1000, 0., 0.001);
00110 
00111     sprintf (histo, "EE crystal energy spectrum" );
00112     meEEcrystalEnergy_ = dbe_->book1D(histo, histo, 5000, 0., 50.);
00113 
00114     sprintf (histo, "EE crystal energy spectrum 2" );
00115     meEEcrystalEnergy2_ = dbe_->book1D(histo, histo, 1000, 0., 0.001);
00116 
00117     sprintf (histo, "EE E1" ) ;
00118     meEEe1_ = dbe_->book1D(histo, histo, 400, 0., 400.);
00119 
00120     sprintf (histo, "EE E4" ) ;
00121     meEEe4_ = dbe_->book1D(histo, histo, 400, 0., 400.);
00122 
00123     sprintf (histo, "EE E9" ) ;
00124     meEEe9_ = dbe_->book1D(histo, histo, 400, 0., 400.);
00125 
00126     sprintf (histo, "EE E16" ) ;
00127     meEEe16_ = dbe_->book1D(histo, histo, 400, 0., 400.);
00128 
00129     sprintf (histo, "EE E25" ) ;
00130     meEEe25_ = dbe_->book1D(histo, histo, 400, 0., 400.);
00131 
00132     sprintf (histo, "EE E1oE4" ) ;
00133     meEEe1oe4_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00134 
00135     sprintf (histo, "EE E1oE9" ) ;
00136     meEEe1oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00137 
00138     sprintf (histo, "EE E4oE9" ) ;
00139     meEEe4oe9_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00140 
00141     sprintf (histo, "EE E9oE16" ) ;
00142     meEEe9oe16_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00143 
00144     sprintf (histo, "EE E1oE25" ) ;
00145     meEEe1oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00146 
00147     sprintf (histo, "EE E9oE25" ) ;
00148     meEEe9oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00149 
00150     sprintf (histo, "EE E16oE25" ) ;
00151     meEEe16oe25_ = dbe_->book1D(histo, histo, 100, 0.4, 1.1);
00152   }
00153 }
00154 
00155 EcalEndcapSimHitsValidation::~EcalEndcapSimHitsValidation(){
00156 
00157 }
00158 
00159 void EcalEndcapSimHitsValidation::beginJob(const edm::EventSetup& c){
00160 
00161 }
00162 
00163 void EcalEndcapSimHitsValidation::endJob(){
00164 
00165   for ( int myStep = 0; myStep<26; myStep++){
00166     if (meEELongitudinalShower_) meEELongitudinalShower_->Fill(float(myStep), eRLength[myStep]/myEntries);
00167   }
00168 
00169 }
00170 
00171 void EcalEndcapSimHitsValidation::analyze(const edm::Event& e, const edm::EventSetup& c){
00172 
00173   edm::LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
00174   
00175   edm::Handle<edm::PCaloHitContainer> EcalHitsEE;
00176   e.getByLabel(g4InfoLabel,EEHitsCollection,EcalHitsEE);
00177 
00178   // Do nothing if no EndCap data available
00179   if( ! EcalHitsEE.isValid() ) return;
00180 
00181   edm::Handle<PEcalValidInfo> MyPEcalValidInfo;
00182   e.getByLabel(g4InfoLabel,ValidationCollection,MyPEcalValidInfo);
00183 
00184   std::vector<PCaloHit> theEECaloHits;
00185   theEECaloHits.insert(theEECaloHits.end(), EcalHitsEE->begin(), EcalHitsEE->end());
00186   
00187   myEntries++;
00188 
00189   std::map<unsigned int, std::vector<PCaloHit*>,std::less<unsigned int> > CaloHitMap;
00190   
00191   double EEetzp_ = 0.;
00192   double EEetzm_ = 0.;
00193   
00194   double ee1  = 0.0;
00195   double ee4  = 0.0;
00196   double ee9  = 0.0;
00197   double ee16 = 0.0;
00198   double ee25 = 0.0;
00199   std::vector<double> econtr(140, 0. );
00200   std::vector<double> econtr25(140, 0. );
00201   
00202   MapType eemap;
00203   MapType eemapzp;
00204   MapType eemapzm;
00205   uint32_t nEEzpHits = 0;
00206   uint32_t nEEzmHits = 0;
00207   
00208   for (std::vector<PCaloHit>::iterator isim = theEECaloHits.begin();
00209        isim != theEECaloHits.end(); ++isim){
00210 
00211     if ( isim->time() > 500. ) { continue; }
00212 
00213     CaloHitMap[ isim->id()].push_back(&(*isim));
00214     
00215     EEDetId eeid (isim->id()) ;
00216     
00217     LogDebug("HitInfo") 
00218       << " CaloHit " << isim->getName() << "\n" 
00219       << " DetID = "<<isim->id()<< " EEDetId = " << eeid.ix() << " " << eeid.iy() << "\n"       
00220       << " Time = " << isim->time() << "\n"
00221       << " Track Id = " << isim->geantTrackId() << "\n"
00222       << " Energy = " << isim->energy();
00223     
00224     uint32_t crystid = eeid.rawId();
00225      
00226     if (eeid.zside() > 0 ) {
00227       nEEzpHits++;
00228       EEetzp_ += isim->energy();
00229       eemapzp[crystid] += isim->energy();
00230       if (meEEzpOccupancy_) meEEzpOccupancy_->Fill(eeid.ix(), eeid.iy());
00231     }
00232     else if (eeid.zside() < 0 ) {
00233       nEEzmHits++;
00234       EEetzm_ += isim->energy();
00235       eemapzm[crystid] += isim->energy();
00236       if (meEEzmOccupancy_) meEEzmOccupancy_->Fill(eeid.ix(), eeid.iy());
00237     }
00238 
00239     if (meEEHitEnergy_) meEEHitEnergy_->Fill(isim->energy());
00240     if( meEEhitLog10Energy_ ) meEEhitLog10Energy_->Fill(log10(isim->energy()));
00241     if (meEEHitEnergy2_) meEEHitEnergy2_->Fill(isim->energy());
00242     eemap[crystid] += isim->energy();
00243 
00244     int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
00245     if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
00246 
00247   }
00248   
00249   if (meEEzpCrystals_) meEEzpCrystals_->Fill(eemapzp.size());
00250   if (meEEzmCrystals_) meEEzmCrystals_->Fill(eemapzm.size());
00251   
00252   if (meEEcrystalEnergy_) {
00253     for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = eemap.begin(); it != eemap.end(); ++it ) meEEcrystalEnergy_->Fill((*it).second);
00254   }
00255   if (meEEcrystalEnergy2_) {
00256     for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = eemap.begin(); it != eemap.end(); ++it ) meEEcrystalEnergy2_->Fill((*it).second);
00257   }
00258     
00259   if (meEEzpHits_)    meEEzpHits_->Fill(nEEzpHits);
00260   if (meEEzmHits_)    meEEzmHits_->Fill(nEEzmHits);
00261     
00262   
00263   int nEEHits = nEEzmHits + nEEzpHits;
00264   if (nEEHits > 0) {
00265     
00266     uint32_t  eecenterid = getUnitWithMaxEnergy(eemap);
00267     EEDetId myEEid(eecenterid);
00268     int bx = myEEid.ix();
00269     int by = myEEid.iy();
00270     int bz = myEEid.zside();
00271     ee1 =  energyInMatrixEE(1,1,bx,by,bz,eemap);
00272     if (meEEe1_) meEEe1_->Fill(ee1);
00273     ee9 =  energyInMatrixEE(3,3,bx,by,bz,eemap);
00274     if (meEEe9_) meEEe9_->Fill(ee9);
00275     ee25=  energyInMatrixEE(5,5,bx,by,bz,eemap);
00276     if (meEEe25_) meEEe25_->Fill(ee25);
00277     
00278     std::vector<uint32_t> ids25; ids25 = getIdsAroundMax(5,5,bx,by,bz,eemap);
00279 
00280     for( unsigned i=0; i<25; i++ ) {
00281       for( unsigned int j=0; j<CaloHitMap[ids25[i]].size(); j++ ) {
00282         int log10i = int( ( log10( CaloHitMap[ids25[i]][j]->energy()) + 10. ) * 10. );
00283         if( log10i >=0 && log10i < 140 ) econtr25[log10i] += CaloHitMap[ids25[i]][j]->energy();
00284       }
00285     }
00286 
00287     MapType  neweemap;
00288     if( fillEEMatrix(3,3,bx,by,bz,neweemap, eemap)){
00289       ee4 = eCluster2x2(neweemap);
00290       if (meEEe4_) meEEe4_->Fill(ee4);
00291     }
00292     if( fillEEMatrix(5,5,bx,by,bz,neweemap, eemap)){
00293       ee16 = eCluster4x4(ee9,neweemap); 
00294       if (meEEe16_) meEEe16_->Fill(ee16);
00295     }
00296     
00297     if (meEEe1oe4_   && ee4  > 0.1 ) meEEe1oe4_  ->Fill(ee1/ee4);
00298     if (meEEe1oe9_   && ee9  > 0.1 ) meEEe1oe9_  ->Fill(ee1/ee9);
00299     if (meEEe4oe9_   && ee9  > 0.1 ) meEEe4oe9_  ->Fill(ee4/ee9);
00300     if (meEEe9oe16_  && ee16 > 0.1 ) meEEe9oe16_ ->Fill(ee9/ee16);
00301     if (meEEe16oe25_ && ee25 > 0.1 ) meEEe16oe25_->Fill(ee16/ee25);
00302     if (meEEe1oe25_  && ee25 > 0.1 ) meEEe1oe25_ ->Fill(ee1/ee25);
00303     if (meEEe9oe25_  && ee25 > 0.1 ) meEEe9oe25_ ->Fill(ee9/ee25);
00304 
00305     if( meEEhitLog10EnergyNorm_ && (EEetzp_+EEetzm_) != 0 ) {
00306       for( int i=0; i<140; i++ ) {
00307         meEEhitLog10EnergyNorm_->Fill( -10.+(float(i)+0.5)/10., econtr[i]/(EEetzp_+EEetzm_) );
00308       }
00309     }
00310 
00311     if( meEEhitLog10Energy25Norm_ && ee25 != 0 ) {
00312       for( int i=0; i<140; i++ ) {
00313         meEEhitLog10Energy25Norm_->Fill( -10.+(float(i)+0.5)/10., econtr25[i]/ee25 );
00314       }
00315     }
00316     
00317   }
00318     
00319   if( MyPEcalValidInfo.isValid() ) {
00320     if ( MyPEcalValidInfo->ee1x1() > 0. ) {
00321       std::vector<float>  EX0 = MyPEcalValidInfo->eX0();
00322       for (int myStep=0; myStep< 26; myStep++ ) { eRLength[myStep] += EX0[myStep]; }
00323     }
00324   }
00325 }
00326 
00327 float EcalEndcapSimHitsValidation::energyInMatrixEE(int nCellInX, int nCellInY,
00328                                         int centralX, int centralY,
00329                                         int centralZ, MapType& themap){
00330   
00331   int   ncristals   = 0;
00332   float totalEnergy = 0.;
00333         
00334   int goBackInX = nCellInX/2;
00335   int goBackInY = nCellInY/2;
00336   int startX    = centralX-goBackInX;
00337   int startY    = centralY-goBackInY;
00338   
00339   for (int ix=startX; ix<startX+nCellInX; ix++) {
00340     for (int iy=startY; iy<startY+nCellInY; iy++) {
00341       uint32_t index ;
00342       
00343       if(EEDetId::validDetId(ix,iy,centralZ)) {
00344         index = EEDetId(ix,iy,centralZ).rawId();
00345       }
00346       else { continue; }
00347 
00348       totalEnergy   += themap[index];
00349       ncristals     += 1;
00350     }
00351   }
00352   
00353   LogDebug("GeomInfo")
00354     << nCellInX << " x " << nCellInY 
00355     << " EE matrix energy = " << totalEnergy
00356     << " for " << ncristals << " crystals";
00357   return totalEnergy;
00358   
00359 }
00360 
00361 std::vector<uint32_t> EcalEndcapSimHitsValidation::getIdsAroundMax( int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType& themap){
00362   
00363   int   ncristals   = 0;
00364   std::vector<uint32_t> ids( nCellInX*nCellInY );
00365         
00366   int goBackInX = nCellInX/2;
00367   int goBackInY = nCellInY/2;
00368   int startX    = centralX-goBackInX;
00369   int startY    = centralY-goBackInY;
00370   
00371   for (int ix=startX; ix<startX+nCellInX; ix++) {
00372     for (int iy=startY; iy<startY+nCellInY; iy++) {
00373       uint32_t index ;
00374       
00375       if(EEDetId::validDetId(ix,iy,centralZ)) {
00376         index = EEDetId(ix,iy,centralZ).rawId();
00377       }
00378       else { continue; }
00379 
00380       ids[ncristals] = index;
00381       ncristals     += 1;
00382     }
00383   }
00384   
00385   return ids;
00386 }
00387 
00388 bool  EcalEndcapSimHitsValidation::fillEEMatrix(int nCellInX, int nCellInY,
00389                                     int CentralX, int CentralY,int CentralZ,
00390                                     MapType& fillmap, MapType&  themap)
00391 {
00392    int goBackInX = nCellInX/2;
00393    int goBackInY = nCellInY/2;
00394 
00395    int startX  =  CentralX - goBackInX;
00396    int startY  =  CentralY - goBackInY;
00397 
00398    int i = 0 ;
00399    for ( int ix = startX; ix < startX+nCellInX; ix ++ ) {
00400 
00401       for( int iy = startY; iy < startY + nCellInY; iy++ ) {
00402 
00403         uint32_t index ;
00404 
00405         if(EEDetId::validDetId(ix,iy,CentralZ)) {
00406           index = EEDetId(ix,iy,CentralZ).rawId();
00407         }
00408         else { continue; }
00409         fillmap[i++] = themap[index];
00410       }
00411    }
00412    uint32_t  centerid = getUnitWithMaxEnergy(themap);
00413 
00414    if ( fillmap[i/2] == themap[centerid] ) 
00415         return true;
00416    else
00417         return false;
00418 }
00419 
00420 
00421 float EcalEndcapSimHitsValidation::eCluster2x2( MapType& themap){
00422   float  E22=0.;
00423   float e012  = themap[0]+themap[1]+themap[2];
00424   float e036  = themap[0]+themap[3]+themap[6];
00425   float e678  = themap[6]+themap[7]+themap[8];
00426   float e258  = themap[2]+themap[5]+themap[8];
00427 
00428   if ( (e012>e678 || e012==e678) && (e036>e258  || e036==e258))
00429     return     E22=themap[0]+themap[1]+themap[3]+themap[4];
00430   else if ( (e012>e678 || e012==e678)  && (e036<e258 || e036==e258) )
00431     return    E22=themap[1]+themap[2]+themap[4]+themap[5];
00432   else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
00433     return     E22=themap[3]+themap[4]+themap[6]+themap[7];
00434   else if ( (e012<e678|| e012==e678)  && (e036<e258|| e036==e258) )
00435     return    E22=themap[4]+themap[5]+themap[7]+themap[8];
00436   else {
00437     return E22;
00438   }
00439 }
00440 
00441 float EcalEndcapSimHitsValidation::eCluster4x4(float e33,  MapType&  themap){
00442   float E44=0.;
00443   float e0_4   = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
00444   float e0_20  = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
00445   float e4_24  = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
00446   float e0_24  = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];
00447   
00448   if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
00449     return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
00450   else if ((e0_4>e0_24 || e0_4==e0_24)  && (e0_20<e4_24 || e0_20==e4_24))
00451     return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
00452   else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
00453     return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
00454   else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
00455     return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
00456   else{
00457     return E44;
00458   }
00459 }
00460 
00461 uint32_t EcalEndcapSimHitsValidation::getUnitWithMaxEnergy(MapType& themap) {
00462   
00463   //look for max
00464   uint32_t unitWithMaxEnergy = 0;
00465   float    maxEnergy = 0.;
00466   
00467   MapType::iterator iter;
00468   for (iter = themap.begin(); iter != themap.end(); iter++) {
00469     
00470     if (maxEnergy < (*iter).second) {
00471       maxEnergy = (*iter).second;       
00472       unitWithMaxEnergy = (*iter).first;
00473     }                           
00474   }
00475   
00476   LogDebug("GeomInfo")
00477     << " max energy of " << maxEnergy 
00478     << " GeV in Unit id " << unitWithMaxEnergy;
00479   return unitWithMaxEnergy;
00480 }
00481 
00482