00001 #include "Validation/EcalHits/interface/EcalSimHitsValidProducer.h"
00002 #include "FWCore/Framework/interface/Event.h"
00003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00004
00005 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
00006 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
00007 #include "SimG4Core/Notification/interface/EndOfEvent.h"
00008 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
00009 #include "SimDataFormats/ValidationFormats/interface/PValidationFormats.h"
00010 #include "DataFormats/Math/interface/Point3D.h"
00011
00012 #include "G4Step.hh"
00013 #include "G4SDManager.hh"
00014 #include "G4PrimaryParticle.hh"
00015 #include "G4PrimaryVertex.hh"
00016
00017 #include <iostream>
00018
00019
00020 EcalSimHitsValidProducer::EcalSimHitsValidProducer(const edm::ParameterSet& iPSet) :
00021 ee1(0.0),ee4(0.0),ee9(0.0),ee16(0.0),ee25(0.0),
00022 eb1(0.0),eb4(0.0),eb9(0.0),eb16(0.0),eb25(0.0),
00023 totalEInEE(0.0),totalEInEB(0),totalEInES(0.0),totalEInEEzp(0.0),totalEInEEzm(0.0),
00024 totalEInESzp(0.0),totalEInESzm(0.0),
00025 totalHits(0),nHitsInEE(0),nHitsInEB(0),nHitsInES(0),nHitsIn1ES(0),nHitsIn2ES(0),
00026 nCrystalInEB(0),nCrystalInEEzp(0),nCrystalInEEzm(0),
00027 nHitsIn1ESzp(0),nHitsIn1ESzm(0),nHitsIn2ESzp(0),nHitsIn2ESzm(0),
00028 thePID(0),
00029 label(iPSet.getUntrackedParameter<std::string>("instanceLabel","EcalValidInfo") )
00030 {
00031 produces<PEcalValidInfo>(label);
00032
00033 for ( int i = 0; i<26; i++ ) {
00034 eBX0[i] = 0.0;
00035 eEX0[i] = 0.0;
00036 }
00037
00038 }
00039
00040
00041 EcalSimHitsValidProducer::~EcalSimHitsValidProducer(){
00042
00043 }
00044
00045 void
00046 EcalSimHitsValidProducer::produce(edm::Event& e, const edm::EventSetup&)
00047 {
00048 std::auto_ptr<PEcalValidInfo> product(new PEcalValidInfo );
00049 fillEventInfo(*product);
00050 e.put(product,label);
00051 }
00052
00053 void
00054 EcalSimHitsValidProducer::fillEventInfo(PEcalValidInfo& product)
00055 {
00056 if( ee1 != 0 ){
00057 product.ee1 = ee1;
00058 product.ee4 = ee4;
00059 product.ee9 = ee9;
00060 product.ee16 = ee16;
00061 product.ee25 = ee25;
00062 for ( int i = 0; i<26; i++ ) {
00063 product.eEX0.push_back( eEX0[i]);
00064 }
00065 }
00066
00067 if( eb1 != 0 ){
00068 product.eb1 = eb1;
00069 product.eb4 = eb4;
00070 product.eb9 = eb9;
00071 product.eb16 = eb16;
00072 product.eb25 = eb25;
00073 for ( int i = 0; i<26; i++ ) {
00074 product.eBX0.push_back( eBX0[i]);
00075 }
00076 }
00077
00078 product.totalEInEE = totalEInEE;
00079 product.totalEInEB = totalEInEB;
00080 product.totalEInES = totalEInES;
00081
00082 product.totalEInEEzp = totalEInEEzp;
00083 product.totalEInEEzm = totalEInEEzm;
00084
00085 product.totalEInESzp = totalEInESzp;
00086 product.totalEInESzm = totalEInESzm;
00087
00088 product.totalHits = totalHits;
00089 product.nHitsInEE = nHitsInEE;
00090 product.nHitsInEB = nHitsInEB;
00091 product.nHitsInES = nHitsInES;
00092 product.nHitsIn1ES = nHitsIn1ES;
00093 product.nHitsIn2ES = nHitsIn2ES;
00094 product.nCrystalInEB = nCrystalInEB;
00095 product.nCrystalInEEzp = nCrystalInEEzp;
00096 product.nCrystalInEEzm = nCrystalInEEzm;
00097
00098 product.nHitsIn1ESzp = nHitsIn1ESzp;
00099 product.nHitsIn1ESzm = nHitsIn1ESzm;
00100 product.nHitsIn2ESzp = nHitsIn2ESzp;
00101 product.nHitsIn2ESzm = nHitsIn2ESzm;
00102
00103
00104 product.eOf1ES = eOf1ES;
00105 product.eOf2ES = eOf2ES;
00106 product.zOfES = zOfES;
00107
00108 product.eOf1ESzp = eOf1ESzp;
00109 product.eOf1ESzm = eOf1ESzm;
00110 product.eOf2ESzp = eOf2ESzp;
00111 product.eOf2ESzm = eOf2ESzm;
00112
00113
00114 product.phiOfEECaloG4Hit = phiOfEECaloG4Hit;
00115 product.etaOfEECaloG4Hit = etaOfEECaloG4Hit;
00116 product.eOfEECaloG4Hit = eOfEECaloG4Hit;
00117 product.eOfEEPlusCaloG4Hit = eOfEEPlusCaloG4Hit;
00118 product.eOfEEMinusCaloG4Hit = eOfEEMinusCaloG4Hit;
00119 product.tOfEECaloG4Hit = tOfEECaloG4Hit;
00120
00121 product.phiOfESCaloG4Hit = phiOfESCaloG4Hit;
00122 product.etaOfESCaloG4Hit = etaOfESCaloG4Hit;
00123 product.eOfESCaloG4Hit = eOfESCaloG4Hit;
00124 product.tOfESCaloG4Hit = tOfESCaloG4Hit;
00125
00126 product.phiOfEBCaloG4Hit = phiOfEBCaloG4Hit;
00127 product.etaOfEBCaloG4Hit = etaOfEBCaloG4Hit;
00128 product.eOfEBCaloG4Hit = eOfEBCaloG4Hit;
00129 product.tOfEBCaloG4Hit = tOfEBCaloG4Hit;
00130
00131
00132 product.theMomentum = theMomentum;
00133 product.theVertex = theVertex;
00134 product.thePID = thePID;
00135 }
00136
00137 void
00138 EcalSimHitsValidProducer::update(const BeginOfEvent*){
00139 ee1 = 0.0;
00140 ee4 = 0.0;
00141 ee9 = 0.0;
00142 ee16 = 0.0;
00143 ee25 = 0.0;
00144
00145 eb1 = 0.0;
00146 eb4 = 0.0;
00147 eb9 = 0.0;
00148 eb16 = 0.0;
00149 eb25 = 0.0;
00150
00151 totalEInEE = 0.0 ;
00152 totalEInEB = 0.0 ;
00153 totalEInES = 0.0 ;
00154
00155 totalEInEEzp = 0.0 ;
00156 totalEInEEzm = 0.0 ;
00157 totalEInESzp = 0.0 ;
00158 totalEInESzm = 0.0 ;
00159
00160 totalHits = 0 ;
00161 nHitsInEE = 0 ;
00162 nHitsInEB = 0 ;
00163 nHitsInES = 0 ;
00164 nHitsIn1ES = 0 ;
00165 nHitsIn2ES = 0 ;
00166 nCrystalInEB = 0;
00167 nCrystalInEEzp = 0;
00168 nCrystalInEEzm = 0;
00169
00170
00171 nHitsIn1ESzp = 0 ;
00172 nHitsIn1ESzm = 0 ;
00173 nHitsIn2ESzp = 0 ;
00174 nHitsIn2ESzm = 0 ;
00175
00176 for ( int i = 0; i<26; i++ ) {
00177 eBX0[i] = 0.0;
00178 eEX0[i] = 0.0;
00179 }
00180
00181 eOf1ES.clear();
00182 eOf2ES.clear();
00183 zOfES.clear();
00184
00185 eOf1ESzp.clear();
00186 eOf1ESzm.clear();
00187 eOf2ESzp.clear();
00188 eOf2ESzm.clear();
00189
00190 phiOfEECaloG4Hit.clear();
00191 etaOfEECaloG4Hit.clear();
00192 tOfEECaloG4Hit.clear();
00193 eOfEECaloG4Hit.clear();
00194 eOfEEPlusCaloG4Hit.clear();
00195 eOfEEMinusCaloG4Hit.clear();
00196
00197 phiOfESCaloG4Hit.clear();
00198 etaOfESCaloG4Hit.clear();
00199 tOfESCaloG4Hit.clear();
00200 eOfESCaloG4Hit.clear();
00201
00202 phiOfEBCaloG4Hit.clear();
00203 etaOfEBCaloG4Hit.clear();
00204 tOfEBCaloG4Hit.clear();
00205 eOfEBCaloG4Hit.clear();
00206 }
00207
00208 void
00209 EcalSimHitsValidProducer::update(const EndOfEvent* evt){
00210
00211 int trackID = 0;
00212 G4PrimaryParticle * thePrim = 0;
00213 int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
00214 if ( nvertex <= 0) {
00215 edm::LogWarning("EcalSimHitsValidProducer")
00216 <<" No Vertex in this Event!";
00217 }else {
00218 for ( int i = 0; i< nvertex; i++){
00219 G4PrimaryVertex * avertex =(*evt)()->GetPrimaryVertex(i);
00220 if ( avertex == 0 )
00221 edm::LogWarning("EcalSimHitsValidProducer")
00222 <<" Pointer to vertex is NULL!";
00223 else {
00224 float x0 = avertex->GetX0();
00225 float y0 = avertex->GetY0();
00226 float z0 = avertex->GetZ0();
00227 float t0 = avertex->GetT0();
00228 theVertex.SetCoordinates(x0,y0,z0,t0);
00229
00230 int npart = avertex->GetNumberOfParticle();
00231 if ( npart == 0)
00232 edm::LogWarning("EcalSimHitsValidProducer")
00233 <<" No primary particle in this event";
00234 else {
00235 if ( thePrim == 0)
00236 thePrim = avertex->GetPrimary(trackID);
00237 }
00238 }
00239 }
00240
00241
00242 double etaInit =0, phiInit =0, pInit =0;
00243 if ( thePrim != 0){
00244 double px = thePrim -> GetPx();
00245 double py = thePrim -> GetPy();
00246 double pz = thePrim -> GetPz();
00247 theMomentum.SetCoordinates(px,py,pz,0.);
00248
00249 pInit =sqrt( pow(px,2.) + pow(py,2.) + pow(pz,2.));
00250 if ( pInit == 0)
00251 edm::LogWarning("EcalSimHitsValidProducer")
00252 <<" Primary has p = 0 ; ";
00253 else {
00254 theMomentum.SetE(pInit);
00255 double costheta = pz/pInit;
00256 double theta = acos(std::min(std::max(costheta, -1.),1.));
00257 etaInit = -log(tan(theta/2));
00258
00259 if ( px != 0 || py != 0)
00260 phiInit = atan2(py,px);
00261 }
00262
00263 thePID = thePrim->GetPDGcode();
00264 }else {
00265 edm::LogWarning("EcalSimHitsValidProducer")
00266 <<" Could not find the primary particle!!";
00267 }
00268 }
00269
00270 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
00271 int EBHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsEB");
00272 int EEHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsEE");
00273 int SEHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsES");
00274
00275
00276 CaloG4HitCollection* theEBHC = (CaloG4HitCollection*) allHC->GetHC(EBHCid);
00277 CaloG4HitCollection* theEEHC = (CaloG4HitCollection*) allHC->GetHC(EEHCid);
00278 CaloG4HitCollection* theSEHC = (CaloG4HitCollection*) allHC->GetHC(SEHCid);
00279
00280 nHitsInEE = theEEHC->entries();
00281 nHitsInEB = theEBHC->entries();
00282 nHitsInES = theSEHC->entries();
00283 totalHits = nHitsInEE + nHitsInEB + nHitsInES;
00284
00285
00286 MapType ebmap;
00287 for (int j=0; j<theEBHC->entries(); j++) {
00288 CaloG4Hit* aHit = (*theEBHC)[j];
00289 totalEInEB += aHit->getEnergyDeposit();
00290 float he = aHit->getEnergyDeposit();
00291 float htime = aHit->getTimeSlice();
00292
00293 math::XYZPoint hpos = aHit->getEntry();
00294 float htheta = hpos.theta();
00295 float heta = -log(tan(htheta * 0.5));
00296 float hphi = hpos.phi();
00297
00298 phiOfEBCaloG4Hit.push_back( hphi);
00299 etaOfEBCaloG4Hit.push_back(heta);
00300 tOfEBCaloG4Hit.push_back(htime);
00301 eOfEBCaloG4Hit.push_back(he);
00302 uint32_t crystid = aHit->getUnitID();
00303 ebmap[crystid] += aHit->getEnergyDeposit();
00304 }
00305
00306 nCrystalInEB = ebmap.size();
00307
00308
00309 MapType eemap,eezpmap,eezmmap;
00310
00311 for (int j=0; j<theEEHC->entries(); j++) {
00312 CaloG4Hit* aHit = (*theEEHC)[j];
00313 totalEInEE += aHit->getEnergyDeposit();
00314 float he = aHit->getEnergyDeposit();
00315 float htime = aHit->getTimeSlice();
00316
00317 math::XYZPoint hpos = aHit->getEntry();
00318 float htheta = hpos.theta();
00319 float heta = -log(tan(htheta * 0.5));
00320 float hphi = hpos.phi();
00321 phiOfEECaloG4Hit.push_back( hphi);
00322 etaOfEECaloG4Hit.push_back(heta);
00323 tOfEECaloG4Hit.push_back(htime);
00324 eOfEECaloG4Hit.push_back(he);
00325
00326 uint32_t crystid = aHit->getUnitID();
00327 EEDetId myEEid(crystid);
00328 if ( myEEid.zside() == -1 ) {
00329 totalEInEEzm += aHit->getEnergyDeposit();
00330 eOfEEMinusCaloG4Hit.push_back(he);
00331 eezmmap[crystid] += aHit->getEnergyDeposit();
00332 }
00333 if ( myEEid.zside() == 1 ) {
00334 totalEInEEzp += aHit->getEnergyDeposit();
00335 eOfEEPlusCaloG4Hit.push_back(he);
00336 eezpmap[crystid] += aHit->getEnergyDeposit();
00337 }
00338
00339
00340 eemap[crystid] += aHit->getEnergyDeposit();
00341 }
00342
00343 nCrystalInEEzm = eezmmap.size();
00344 nCrystalInEEzp = eezpmap.size();
00345
00346
00347 for (int j=0; j<theSEHC->entries(); j++) {
00348 CaloG4Hit* aHit = (*theSEHC)[j];
00349 totalEInES += aHit->getEnergyDeposit();
00350 ESDetId esid = ESDetId( aHit->getUnitID());
00351
00352 if (esid.zside() == -1) {
00353
00354 totalEInESzm += aHit->getEnergyDeposit();
00355
00356 if (esid.plane() == 1) {
00357 nHitsIn1ESzm++;
00358 eOf1ESzm.push_back(aHit->getEnergyDeposit());
00359 }else if ( esid.plane() == 2){
00360 nHitsIn2ESzm++;
00361 eOf2ESzm.push_back(aHit->getEnergyDeposit());
00362 }
00363 }
00364 if (esid.zside() == 1) {
00365
00366 totalEInESzp += aHit->getEnergyDeposit();
00367
00368 if (esid.plane() == 1) {
00369 nHitsIn1ESzp++;
00370 eOf1ESzp.push_back(aHit->getEnergyDeposit());
00371 }else if ( esid.plane() == 2){
00372 nHitsIn2ESzp++;
00373 eOf2ESzp.push_back(aHit->getEnergyDeposit());
00374 }
00375 }
00376
00377
00378
00379 }
00380
00381 uint32_t eemaxid = getUnitWithMaxEnergy(eemap);
00382 uint32_t ebmaxid = getUnitWithMaxEnergy(ebmap);
00383 if ( eemap[eemaxid] > ebmap[ebmaxid] ) {
00384 uint32_t centerid = getUnitWithMaxEnergy(eemap);
00385 EEDetId myEEid(centerid);
00386 int ix = myEEid.ix();
00387 int iy = myEEid.iy();
00388 int iz = myEEid.zside();
00389
00390 ee1 = energyInEEMatrix(1,1,ix,iy,iz,eemap);
00391 ee9 = energyInEEMatrix(3,3,ix,iy,iz,eemap);
00392 ee25= energyInEEMatrix(5,5,ix,iy,iz,eemap);
00393 MapType neweemap;
00394 if( fillEEMatrix(3,3,ix,iy,iz,neweemap, eemap)) {
00395 ee4 = eCluster2x2(neweemap);
00396 }
00397 if ( fillEEMatrix(5,5,ix,iy,iz,neweemap, eemap)) {
00398 ee16 = eCluster4x4(ee9,neweemap);
00399 }
00400 } else {
00401 uint32_t ebcenterid = getUnitWithMaxEnergy(ebmap);
00402 EBDetId myEBid(ebcenterid);
00403 int bx = myEBid.ietaAbs();
00404 int by = myEBid.iphi();
00405 int bz = myEBid.zside();
00406 eb1 = energyInEBMatrix(1,1,bx,by,bz,ebmap);
00407 eb9 = energyInEBMatrix(3,3,bx,by,bz,ebmap);
00408 eb25= energyInEBMatrix(5,5,bx,by,bz,ebmap);
00409
00410 MapType newebmap;
00411 if( fillEBMatrix(3,3,bx,by,bz,newebmap, ebmap)){
00412 eb4 = eCluster2x2(newebmap);
00413 }
00414 if( fillEBMatrix(5,5,bx,by,bz,newebmap, ebmap)){
00415 eb16 = eCluster4x4(eb9,newebmap);
00416 }
00417 }
00418
00419 }
00420
00421 void
00422 EcalSimHitsValidProducer::update(const G4Step* aStep){
00423 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
00424 G4ThreeVector hitPoint = preStepPoint->GetPosition();
00425 G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
00426 G4String name = currentPV->GetName();
00427 std::string crystal;
00428 crystal.assign(name,0,4);
00429
00430 float Edeposit = aStep->GetTotalEnergyDeposit();
00431 if ( crystal == "EFRY"&& Edeposit > 0.0){
00432 float z = hitPoint.z();
00433 float detz = fabs(fabs(z)-3200);
00434 int x0 = (int)floor( detz/8.9 );
00435 if ( x0< 26){
00436 eEX0[x0] += Edeposit;
00437 }
00438 }
00439 if(crystal == "EBRY" && Edeposit > 0.0) {
00440 float x = hitPoint.x();
00441 float y = hitPoint.y();
00442 float r = sqrt(x*x +y*y);
00443 float detr = r -1290;
00444 int x0 = (int)floor( detr/8.9);
00445 eBX0[x0] += Edeposit;
00446 }
00447
00448 }
00449
00450 bool EcalSimHitsValidProducer::fillEEMatrix(int nCellInEta, int nCellInPhi,
00451 int CentralEta, int CentralPhi,int CentralZ,
00452 MapType& fillmap, MapType& themap)
00453 {
00454 int goBackInEta = nCellInEta/2;
00455 int goBackInPhi = nCellInPhi/2;
00456
00457 int startEta = CentralEta - goBackInEta;
00458 int startPhi = CentralPhi - goBackInPhi;
00459
00460 int i = 0 ;
00461 for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
00462
00463 for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
00464 uint32_t index;
00465
00466 if (EEDetId::validDetId(ieta, iphi,CentralZ)) {
00467 index = EEDetId( ieta, iphi,CentralZ).rawId();
00468 } else { continue; }
00469
00470 fillmap[i++] = themap[index];
00471 }
00472 }
00473 uint32_t centerid = getUnitWithMaxEnergy(themap);
00474
00475 if ( fillmap[i/2] == themap[centerid] )
00476 return true;
00477 else
00478 return false;
00479 }
00480
00481 bool EcalSimHitsValidProducer::fillEBMatrix(int nCellInEta, int nCellInPhi,
00482 int CentralEta, int CentralPhi,int CentralZ,
00483 MapType& fillmap, MapType& themap)
00484 {
00485 int goBackInEta = nCellInEta/2;
00486 int goBackInPhi = nCellInPhi/2;
00487
00488 int startEta = CentralZ*CentralEta - goBackInEta;
00489 int startPhi = CentralPhi - goBackInPhi;
00490
00491 int i = 0 ;
00492 for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
00493 for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
00494 uint32_t index;
00495 if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00496 if (iphi< 1) { index = EBDetId(ieta,iphi+360).rawId(); }
00497 else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00498 else { index = EBDetId(ieta,iphi).rawId(); }
00499 fillmap[i++] = themap[index];
00500 }
00501 }
00502
00503 uint32_t ebcenterid = getUnitWithMaxEnergy(themap);
00504
00505 if ( fillmap[i/2] == themap[ebcenterid] )
00506 return true;
00507 else
00508 return false;
00509 }
00510
00511
00512 float EcalSimHitsValidProducer::eCluster2x2( MapType& themap){
00513 float E22=0.;
00514 float e012 = themap[0]+themap[1]+themap[2];
00515 float e036 = themap[0]+themap[3]+themap[6];
00516 float e678 = themap[6]+themap[7]+themap[8];
00517 float e258 = themap[2]+themap[5]+themap[8];
00518
00519 if ( (e012>e678 || e012==e678) && (e036>e258 || e036==e258))
00520 return E22=themap[0]+themap[1]+themap[3]+themap[4];
00521 else if ( (e012>e678 || e012==e678) && (e036<e258 || e036==e258) )
00522 return E22=themap[1]+themap[2]+themap[4]+themap[5];
00523 else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
00524 return E22=themap[3]+themap[4]+themap[6]+themap[7];
00525 else if ( (e012<e678|| e012==e678) && (e036<e258|| e036==e258) )
00526 return E22=themap[4]+themap[5]+themap[7]+themap[8];
00527 else {
00528 return E22;
00529 }
00530 }
00531
00532 float EcalSimHitsValidProducer::eCluster4x4(float e33, MapType& themap){
00533 float E44=0.;
00534 float e0_4 = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
00535 float e0_20 = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
00536 float e4_24 = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
00537 float e0_24 = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];
00538
00539 if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
00540 return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
00541 else if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
00542 return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
00543 else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
00544 return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
00545 else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
00546 return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
00547 else{
00548 return E44;
00549 }
00550 }
00551
00552 float EcalSimHitsValidProducer::energyInEEMatrix(int nCellInX, int nCellInY,
00553 int centralX, int centralY,
00554 int centralZ, MapType& themap){
00555
00556 int ncristals = 0;
00557 float totalEnergy = 0.;
00558
00559 int goBackInX = nCellInX/2;
00560 int goBackInY = nCellInY/2;
00561 int startX = centralX-goBackInX;
00562 int startY = centralY-goBackInY;
00563
00564 for (int ix=startX; ix<startX+nCellInX; ix++) {
00565 for (int iy=startY; iy<startY+nCellInY; iy++) {
00566
00567 uint32_t index ;
00568 if (EEDetId::validDetId(ix, iy,centralZ)) {
00569 index = EEDetId(ix,iy,centralZ).rawId();
00570 } else { continue; }
00571
00572 totalEnergy += themap[index];
00573 ncristals += 1;
00574
00575 LogDebug("EcalSimHitsValidProducer")
00576 << " EnergyInEEMatrix: ix - iy - E = " << ix
00577 << " " << iy << " " << themap[index] ;
00578
00579 }
00580 }
00581
00582 LogDebug("EcalSimHitsValidProducer")
00583 << " EnergyInEEMatrix: energy in " << nCellInX
00584 << " cells in x times " << nCellInY
00585 << " cells in y matrix = " << totalEnergy
00586 << " for " << ncristals << " crystals";
00587
00588 return totalEnergy;
00589
00590 }
00591
00592 float EcalSimHitsValidProducer::energyInEBMatrix(int nCellInEta, int nCellInPhi,
00593 int centralEta, int centralPhi,
00594 int centralZ, MapType& themap){
00595
00596 int ncristals = 0;
00597 float totalEnergy = 0.;
00598
00599 int goBackInEta = nCellInEta/2;
00600 int goBackInPhi = nCellInPhi/2;
00601 int startEta = centralZ*centralEta-goBackInEta;
00602 int startPhi = centralPhi-goBackInPhi;
00603
00604 for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
00605 for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
00606
00607 uint32_t index ;
00608 if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00609 if (iphi< 1) { index = EBDetId(ieta,iphi+360).rawId(); }
00610 else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00611 else { index = EBDetId(ieta,iphi).rawId(); }
00612
00613 totalEnergy += themap[index];
00614 ncristals += 1;
00615
00616 LogDebug("EcalSimHitsValidProducer")
00617 << " EnergyInEBMatrix: ieta - iphi - E = " << ieta
00618 << " " << iphi << " " << themap[index];
00619
00620 }
00621 }
00622
00623 LogDebug("EcalSimHitsValidProducer")
00624 << " EnergyInEBMatrix: energy in " << nCellInEta
00625 << " cells in eta times " << nCellInPhi
00626 << " cells in phi matrix = " << totalEnergy
00627 << " for " << ncristals << " crystals";
00628
00629 return totalEnergy;
00630
00631 }
00632
00633 uint32_t EcalSimHitsValidProducer::getUnitWithMaxEnergy(MapType& themap) {
00634
00635 uint32_t unitWithMaxEnergy = 0;
00636 float maxEnergy = 0.;
00637
00638 MapType::iterator iter;
00639 for (iter = themap.begin(); iter != themap.end(); iter++) {
00640
00641 if (maxEnergy < (*iter).second) {
00642 maxEnergy = (*iter).second;
00643 unitWithMaxEnergy = (*iter).first;
00644 }
00645 }
00646
00647
00648 LogDebug("EcalSimHitsValidProducer")
00649 << " Max energy of " << maxEnergy
00650 << " MeV was found in Unit id 0x" << std::hex
00651 << unitWithMaxEnergy << std::dec;
00652
00653 return unitWithMaxEnergy;
00654 }
00655
1.7.3