00001
00002
00003
00004
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
00023 verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
00024
00025
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
00166
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
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
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
1.7.3