Reference Manual

Go to the documentation of this file.
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(){
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( isim->energy() > 0 ) {
00241       if( meEEhitLog10Energy_ ) meEEhitLog10Energy_->Fill(log10(isim->energy()));
00242       int log10i = int( ( log10(isim->energy()) + 10. ) * 10. );
00243       if( log10i >=0 && log10i < 140 ) econtr[log10i] += isim->energy();
00244     }
00245     if (meEEHitEnergy2_) meEEHitEnergy2_->Fill(isim->energy());
00246     eemap[crystid] += isim->energy();
00247 
00248 
00249   }
00250   
00251   if (meEEzpCrystals_) meEEzpCrystals_->Fill(eemapzp.size());
00252   if (meEEzmCrystals_) meEEzmCrystals_->Fill(eemapzm.size());
00253   
00254   if (meEEcrystalEnergy_) {
00255     for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = eemap.begin(); it != eemap.end(); ++it ) meEEcrystalEnergy_->Fill((*it).second);
00256   }
00257   if (meEEcrystalEnergy2_) {
00258     for (std::map<uint32_t,float,std::less<uint32_t> >::iterator it = eemap.begin(); it != eemap.end(); ++it ) meEEcrystalEnergy2_->Fill((*it).second);
00259   }
00260     
00261   if (meEEzpHits_)    meEEzpHits_->Fill(nEEzpHits);
00262   if (meEEzmHits_)    meEEzmHits_->Fill(nEEzmHits);
00263     
00264   
00265   int nEEHits = nEEzmHits + nEEzpHits;
00266   if (nEEHits > 0) {
00267     
00268     uint32_t  eecenterid = getUnitWithMaxEnergy(eemap);
00269     EEDetId myEEid(eecenterid);
00270     int bx = myEEid.ix();
00271     int by = myEEid.iy();
00272     int bz = myEEid.zside();
00273     ee1 =  energyInMatrixEE(1,1,bx,by,bz,eemap);
00274     if (meEEe1_) meEEe1_->Fill(ee1);
00275     ee9 =  energyInMatrixEE(3,3,bx,by,bz,eemap);
00276     if (meEEe9_) meEEe9_->Fill(ee9);
00277     ee25=  energyInMatrixEE(5,5,bx,by,bz,eemap);
00278     if (meEEe25_) meEEe25_->Fill(ee25);
00279     
00280     std::vector<uint32_t> ids25; ids25 = getIdsAroundMax(5,5,bx,by,bz,eemap);
00281 
00282     for( unsigned i=0; i<25; i++ ) {
00283       for( unsigned int j=0; j<CaloHitMap[ids25[i]].size(); j++ ) {
00284         if( CaloHitMap[ids25[i]][j]->energy() > 0 ) {
00285           int log10i = int( ( log10( CaloHitMap[ids25[i]][j]->energy()) + 10. ) * 10. );
00286           if( log10i >=0 && log10i < 140 ) econtr25[log10i] += CaloHitMap[ids25[i]][j]->energy();
00287         }
00288       }
00289     }
00290 
00291     MapType  neweemap;
00292     if( fillEEMatrix(3,3,bx,by,bz,neweemap, eemap)){
00293       ee4 = eCluster2x2(neweemap);
00294       if (meEEe4_) meEEe4_->Fill(ee4);
00295     }
00296     if( fillEEMatrix(5,5,bx,by,bz,neweemap, eemap)){
00297       ee16 = eCluster4x4(ee9,neweemap); 
00298       if (meEEe16_) meEEe16_->Fill(ee16);
00299     }
00300     
00301     if (meEEe1oe4_   && ee4  > 0.1 ) meEEe1oe4_  ->Fill(ee1/ee4);
00302     if (meEEe1oe9_   && ee9  > 0.1 ) meEEe1oe9_  ->Fill(ee1/ee9);
00303     if (meEEe4oe9_   && ee9  > 0.1 ) meEEe4oe9_  ->Fill(ee4/ee9);
00304     if (meEEe9oe16_  && ee16 > 0.1 ) meEEe9oe16_ ->Fill(ee9/ee16);
00305     if (meEEe16oe25_ && ee25 > 0.1 ) meEEe16oe25_->Fill(ee16/ee25);
00306     if (meEEe1oe25_  && ee25 > 0.1 ) meEEe1oe25_ ->Fill(ee1/ee25);
00307     if (meEEe9oe25_  && ee25 > 0.1 ) meEEe9oe25_ ->Fill(ee9/ee25);
00308 
00309     if( meEEhitLog10EnergyNorm_ && (EEetzp_+EEetzm_) != 0 ) {
00310       for( int i=0; i<140; i++ ) {
00311         meEEhitLog10EnergyNorm_->Fill( -10.+(float(i)+0.5)/10., econtr[i]/(EEetzp_+EEetzm_) );
00312       }
00313     }
00314 
00315     if( meEEhitLog10Energy25Norm_ && ee25 != 0 ) {
00316       for( int i=0; i<140; i++ ) {
00317         meEEhitLog10Energy25Norm_->Fill( -10.+(float(i)+0.5)/10., econtr25[i]/ee25 );
00318       }
00319     }
00320     
00321   }
00322     
00323   if( MyPEcalValidInfo.isValid() ) {
00324     if ( MyPEcalValidInfo->ee1x1() > 0. ) {
00325       std::vector<float>  EX0 = MyPEcalValidInfo->eX0();
00326       if (meEELongitudinalShower_) meEELongitudinalShower_->Reset();
00327       for (int myStep=0; myStep< 26; myStep++ ) { 
00328         eRLength[myStep] += EX0[myStep]; 
00329         if (meEELongitudinalShower_) meEELongitudinalShower_->Fill(float(myStep), eRLength[myStep]/myEntries);
00330       }
00331     }
00332   }
00333 }
00334 
00335 float EcalEndcapSimHitsValidation::energyInMatrixEE(int nCellInX, int nCellInY,
00336                                         int centralX, int centralY,
00337                                         int centralZ, MapType& themap){
00338   
00339   int   ncristals   = 0;
00340   float totalEnergy = 0.;
00341         
00342   int goBackInX = nCellInX/2;
00343   int goBackInY = nCellInY/2;
00344   int startX    = centralX-goBackInX;
00345   int startY    = centralY-goBackInY;
00346   
00347   for (int ix=startX; ix<startX+nCellInX; ix++) {
00348     for (int iy=startY; iy<startY+nCellInY; iy++) {
00349       uint32_t index ;
00350       
00351       if(EEDetId::validDetId(ix,iy,centralZ)) {
00352         index = EEDetId(ix,iy,centralZ).rawId();
00353       }
00354       else { continue; }
00355 
00356       totalEnergy   += themap[index];
00357       ncristals     += 1;
00358     }
00359   }
00360   
00361   LogDebug("GeomInfo")
00362     << nCellInX << " x " << nCellInY 
00363     << " EE matrix energy = " << totalEnergy
00364     << " for " << ncristals << " crystals";
00365   return totalEnergy;
00366   
00367 }
00368 
00369 std::vector<uint32_t> EcalEndcapSimHitsValidation::getIdsAroundMax( int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType& themap){
00370   
00371   int   ncristals   = 0;
00372   std::vector<uint32_t> ids( nCellInX*nCellInY );
00373         
00374   int goBackInX = nCellInX/2;
00375   int goBackInY = nCellInY/2;
00376   int startX    = centralX-goBackInX;
00377   int startY    = centralY-goBackInY;
00378   
00379   for (int ix=startX; ix<startX+nCellInX; ix++) {
00380     for (int iy=startY; iy<startY+nCellInY; iy++) {
00381       uint32_t index ;
00382       
00383       if(EEDetId::validDetId(ix,iy,centralZ)) {
00384         index = EEDetId(ix,iy,centralZ).rawId();
00385       }
00386       else { continue; }
00387 
00388       ids[ncristals] = index;
00389       ncristals     += 1;
00390     }
00391   }
00392   
00393   return ids;
00394 }
00395 
00396 bool  EcalEndcapSimHitsValidation::fillEEMatrix(int nCellInX, int nCellInY,
00397                                     int CentralX, int CentralY,int CentralZ,
00398                                     MapType& fillmap, MapType&  themap)
00399 {
00400    int goBackInX = nCellInX/2;
00401    int goBackInY = nCellInY/2;
00402 
00403    int startX  =  CentralX - goBackInX;
00404    int startY  =  CentralY - goBackInY;
00405 
00406    int i = 0 ;
00407    for ( int ix = startX; ix < startX+nCellInX; ix ++ ) {
00408 
00409       for( int iy = startY; iy < startY + nCellInY; iy++ ) {
00410 
00411         uint32_t index ;
00412 
00413         if(EEDetId::validDetId(ix,iy,CentralZ)) {
00414           index = EEDetId(ix,iy,CentralZ).rawId();
00415         }
00416         else { continue; }
00417         fillmap[i++] = themap[index];
00418       }
00419    }
00420    uint32_t  centerid = getUnitWithMaxEnergy(themap);
00421 
00422    if ( fillmap[i/2] == themap[centerid] ) 
00423         return true;
00424    else
00425         return false;
00426 }
00427 
00428 
00429 float EcalEndcapSimHitsValidation::eCluster2x2( MapType& themap){
00430   float  E22=0.;
00431   float e012  = themap[0]+themap[1]+themap[2];
00432   float e036  = themap[0]+themap[3]+themap[6];
00433   float e678  = themap[6]+themap[7]+themap[8];
00434   float e258  = themap[2]+themap[5]+themap[8];
00435 
00436   if ( (e012>e678 || e012==e678) && (e036>e258  || e036==e258))
00437     return     E22=themap[0]+themap[1]+themap[3]+themap[4];
00438   else if ( (e012>e678 || e012==e678)  && (e036<e258 || e036==e258) )
00439     return    E22=themap[1]+themap[2]+themap[4]+themap[5];
00440   else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
00441     return     E22=themap[3]+themap[4]+themap[6]+themap[7];
00442   else if ( (e012<e678|| e012==e678)  && (e036<e258|| e036==e258) )
00443     return    E22=themap[4]+themap[5]+themap[7]+themap[8];
00444   else {
00445     return E22;
00446   }
00447 }
00448 
00449 float EcalEndcapSimHitsValidation::eCluster4x4(float e33,  MapType&  themap){
00450   float E44=0.;
00451   float e0_4   = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
00452   float e0_20  = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
00453   float e4_24  = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
00454   float e0_24  = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];
00455   
00456   if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
00457     return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
00458   else if ((e0_4>e0_24 || e0_4==e0_24)  && (e0_20<e4_24 || e0_20==e4_24))
00459     return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
00460   else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
00461     return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
00462   else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
00463     return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
00464   else{
00465     return E44;
00466   }
00467 }
00468 
00469 uint32_t EcalEndcapSimHitsValidation::getUnitWithMaxEnergy(MapType& themap) {
00470   
00471   //look for max
00472   uint32_t unitWithMaxEnergy = 0;
00473   float    maxEnergy = 0.;
00474   
00475   MapType::iterator iter;
00476   for (iter = themap.begin(); iter != themap.end(); iter++) {
00477     
00478     if (maxEnergy < (*iter).second) {
00479       maxEnergy = (*iter).second;       
00480       unitWithMaxEnergy = (*iter).first;
00481     }                           
00482   }
00483   
00484   LogDebug("GeomInfo")
00485     << " max energy of " << maxEnergy 
00486     << " GeV in Unit id " << unitWithMaxEnergy;
00487   return unitWithMaxEnergy;
00488 }
00489 
00490
/data/refman/pasoursint/CMSSW_5_3_10/src/Validation/EcalHits/src/EcalEndcapSimHitsValidation.cc
