CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Validation/EcalHits/src/EcalSimHitsValidProducer.cc

Go to the documentation of this file.
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    // the direction of momentum of primary particles
00242     double pInit =0; // etaInit =0, phiInit =0, // UNUSED
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; // UNUSED
00256                   // double theta = acos(std::min(std::max(costheta, -1.),1.)); // UNUSED
00257                   // etaInit = -log(tan(theta/2)); // UNUSED
00258 
00259                   // if ( px != 0 || py != 0) phiInit = atan2(py,px); // UNUSED
00260           }
00261 
00262           thePID = thePrim->GetPDGcode();
00263       }else {
00264           edm::LogWarning("EcalSimHitsValidProducer")
00265           <<" Could not find the primary particle!!";
00266       }
00267   }
00268   // hit map for EB for matrices
00269   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
00270   int EBHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsEB");
00271   int EEHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsEE");
00272   int SEHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsES");
00273  
00274 
00275   CaloG4HitCollection* theEBHC = (CaloG4HitCollection*) allHC->GetHC(EBHCid);
00276   CaloG4HitCollection* theEEHC = (CaloG4HitCollection*) allHC->GetHC(EEHCid);
00277   CaloG4HitCollection* theSEHC = (CaloG4HitCollection*) allHC->GetHC(SEHCid);
00278   
00279   nHitsInEE = theEEHC->entries();
00280   nHitsInEB = theEBHC->entries();
00281   nHitsInES = theSEHC->entries();
00282   totalHits = nHitsInEE + nHitsInEB + nHitsInES;
00283 
00284  //   EB Hit collection start
00285  MapType ebmap;
00286  for (int j=0; j<theEBHC->entries(); j++) {
00287    CaloG4Hit* aHit = (*theEBHC)[j];
00288     totalEInEB += aHit->getEnergyDeposit();
00289     float he     = aHit->getEnergyDeposit();
00290     float htime  = aHit->getTimeSlice();
00291                                                                                                                              
00292     math::XYZPoint hpos = aHit->getEntry();
00293     float htheta    = hpos.theta();
00294     float heta    = -log(tan(htheta * 0.5));
00295     float hphi    = hpos.phi();
00296 
00297     phiOfEBCaloG4Hit.push_back( hphi);
00298     etaOfEBCaloG4Hit.push_back(heta);
00299     tOfEBCaloG4Hit.push_back(htime);
00300     eOfEBCaloG4Hit.push_back(he);
00301     uint32_t  crystid = aHit->getUnitID();
00302     ebmap[crystid] += aHit->getEnergyDeposit();
00303 }
00304 
00305     nCrystalInEB = ebmap.size();
00306 
00307  //   EE Hit collection start
00308  MapType eemap,eezpmap,eezmmap;
00309  
00310  for (int j=0; j<theEEHC->entries(); j++) {
00311     CaloG4Hit* aHit = (*theEEHC)[j];
00312     totalEInEE += aHit->getEnergyDeposit();
00313     float he     = aHit->getEnergyDeposit();
00314     float htime  = aHit->getTimeSlice();
00315 
00316     math::XYZPoint hpos = aHit->getEntry();
00317     float htheta    = hpos.theta();
00318     float heta    = -log(tan(htheta * 0.5));
00319     float hphi    = hpos.phi();
00320     phiOfEECaloG4Hit.push_back( hphi);
00321     etaOfEECaloG4Hit.push_back(heta);
00322     tOfEECaloG4Hit.push_back(htime);
00323     eOfEECaloG4Hit.push_back(he);
00324  
00325     uint32_t  crystid = aHit->getUnitID();
00326     EEDetId myEEid(crystid);
00327     if ( myEEid.zside() == -1 ) {
00328           totalEInEEzm += aHit->getEnergyDeposit();
00329           eOfEEMinusCaloG4Hit.push_back(he);
00330           eezmmap[crystid] += aHit->getEnergyDeposit();
00331      }
00332     if ( myEEid.zside() == 1 ) {
00333           totalEInEEzp += aHit->getEnergyDeposit();
00334           eOfEEPlusCaloG4Hit.push_back(he);
00335           eezpmap[crystid] += aHit->getEnergyDeposit();
00336      }
00337                
00338 
00339     eemap[crystid] += aHit->getEnergyDeposit(); 
00340  }
00341 
00342  nCrystalInEEzm = eezmmap.size();
00343  nCrystalInEEzp = eezpmap.size(); 
00344 
00345  //Hits from ES
00346  for (int j=0; j<theSEHC->entries(); j++) {
00347     CaloG4Hit* aHit = (*theSEHC)[j];
00348     totalEInES += aHit->getEnergyDeposit();
00349     ESDetId esid = ESDetId( aHit->getUnitID());
00350 
00351     if (esid.zside() == -1) {
00352 
00353            totalEInESzm +=  aHit->getEnergyDeposit();
00354 
00355            if (esid.plane() ==  1) {
00356                  nHitsIn1ESzm++;
00357                  eOf1ESzm.push_back(aHit->getEnergyDeposit());
00358             }else if ( esid.plane() ==  2){
00359                   nHitsIn2ESzm++;
00360                   eOf2ESzm.push_back(aHit->getEnergyDeposit());
00361             }
00362     } 
00363     if (esid.zside() == 1) {
00364 
00365            totalEInESzp +=  aHit->getEnergyDeposit();
00366 
00367            if (esid.plane() ==  1) {
00368                  nHitsIn1ESzp++;
00369                  eOf1ESzp.push_back(aHit->getEnergyDeposit());
00370             }else if ( esid.plane() ==  2){
00371                   nHitsIn2ESzp++;
00372                   eOf2ESzp.push_back(aHit->getEnergyDeposit());
00373             }
00374     }
00375 
00376 
00377   
00378   }
00379   
00380    uint32_t  eemaxid = getUnitWithMaxEnergy(eemap);
00381    uint32_t  ebmaxid = getUnitWithMaxEnergy(ebmap); 
00382    if ( eemap[eemaxid] > ebmap[ebmaxid] ) {
00383     uint32_t  centerid = getUnitWithMaxEnergy(eemap);
00384     EEDetId myEEid(centerid);
00385     int ix = myEEid.ix();
00386     int iy = myEEid.iy();
00387     int iz = myEEid.zside();
00388 
00389     ee1 =  energyInEEMatrix(1,1,ix,iy,iz,eemap);
00390     ee9 =  energyInEEMatrix(3,3,ix,iy,iz,eemap);
00391     ee25=  energyInEEMatrix(5,5,ix,iy,iz,eemap);
00392     MapType   neweemap;
00393     if( fillEEMatrix(3,3,ix,iy,iz,neweemap, eemap)) {
00394           ee4 = eCluster2x2(neweemap);
00395     }
00396     if ( fillEEMatrix(5,5,ix,iy,iz,neweemap, eemap)) {
00397           ee16 = eCluster4x4(ee9,neweemap);
00398     }
00399   } else {
00400     uint32_t  ebcenterid = getUnitWithMaxEnergy(ebmap);
00401     EBDetId myEBid(ebcenterid);
00402     int bx = myEBid.ietaAbs();
00403     int by = myEBid.iphi();
00404     int bz = myEBid.zside();
00405     eb1 =  energyInEBMatrix(1,1,bx,by,bz,ebmap);
00406     eb9 =  energyInEBMatrix(3,3,bx,by,bz,ebmap);
00407     eb25=  energyInEBMatrix(5,5,bx,by,bz,ebmap);
00408 
00409     MapType  newebmap;
00410     if( fillEBMatrix(3,3,bx,by,bz,newebmap, ebmap)){
00411           eb4 = eCluster2x2(newebmap);
00412     }
00413     if( fillEBMatrix(5,5,bx,by,bz,newebmap, ebmap)){
00414           eb16 = eCluster4x4(eb9,newebmap);
00415     }
00416   } 
00417 
00418 }
00419 
00420 void
00421 EcalSimHitsValidProducer::update(const G4Step* aStep){
00422  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
00423   G4ThreeVector hitPoint = preStepPoint->GetPosition();
00424   G4VPhysicalVolume* currentPV  = preStepPoint->GetPhysicalVolume();
00425   G4String name         = currentPV->GetName();
00426   std::string crystal;
00427   crystal.assign(name,0,4);
00428 
00429   float Edeposit =  aStep->GetTotalEnergyDeposit();
00430   if ( crystal == "EFRY"&& Edeposit > 0.0){
00431           float z = hitPoint.z();
00432           float detz = fabs(fabs(z)-3200);
00433           int x0 = (int)floor( detz/8.9 );
00434           if ( x0< 26){ 
00435             eEX0[x0] += Edeposit;
00436           }
00437   } 
00438   if(crystal == "EBRY" && Edeposit > 0.0) {
00439             float x = hitPoint.x();
00440             float y = hitPoint.y();
00441             float r = sqrt(x*x +y*y);
00442             float detr = r -1290;
00443             int x0 = (int)floor( detr/8.9);
00444             eBX0[x0] += Edeposit;
00445   }
00446    
00447 }
00448 
00449 bool  EcalSimHitsValidProducer::fillEEMatrix(int nCellInEta, int nCellInPhi,
00450                                    int CentralEta, int CentralPhi,int CentralZ,
00451                                    MapType& fillmap, MapType&  themap)
00452 {
00453    int goBackInEta = nCellInEta/2;
00454    int goBackInPhi = nCellInPhi/2;
00455 
00456    int startEta  =  CentralEta - goBackInEta;
00457    int startPhi  =  CentralPhi - goBackInPhi;
00458 
00459    int i = 0 ;
00460    for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
00461 
00462       for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
00463         uint32_t index;
00464 
00465         if (EEDetId::validDetId(ieta, iphi,CentralZ)) {
00466           index = EEDetId( ieta, iphi,CentralZ).rawId();
00467         } else { continue; }
00468         
00469         fillmap[i++] = themap[index];
00470       }
00471    }
00472    uint32_t  centerid = getUnitWithMaxEnergy(themap);
00473 
00474    if ( fillmap[i/2] == themap[centerid] ) 
00475         return true;
00476    else
00477         return false;
00478 }
00479 
00480 bool  EcalSimHitsValidProducer::fillEBMatrix(int nCellInEta, int nCellInPhi,
00481                                    int CentralEta, int CentralPhi,int CentralZ,
00482                                    MapType& fillmap, MapType&  themap)
00483 {
00484    int goBackInEta = nCellInEta/2;
00485    int goBackInPhi = nCellInPhi/2;
00486 
00487    int startEta  =  CentralZ*CentralEta - goBackInEta;
00488    int startPhi  =  CentralPhi - goBackInPhi;
00489 
00490    int i = 0 ;
00491    for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
00492       for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
00493           uint32_t  index;
00494           if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00495           if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
00496           else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00497           else              { index = EBDetId(ieta,iphi).rawId();     }
00498           fillmap[i++] = themap[index];
00499       }
00500    }
00501 
00502    uint32_t  ebcenterid = getUnitWithMaxEnergy(themap);
00503 
00504    if ( fillmap[i/2] == themap[ebcenterid] )
00505         return true;
00506    else
00507         return false;
00508 }
00509 
00510 
00511 float EcalSimHitsValidProducer::eCluster2x2( MapType& themap){
00512         float  E22=0.;
00513         float e012  = themap[0]+themap[1]+themap[2];
00514         float e036  = themap[0]+themap[3]+themap[6];
00515         float e678  = themap[6]+themap[7]+themap[8];
00516         float e258  = themap[2]+themap[5]+themap[8];
00517 
00518     if ( (e012>e678 || e012==e678) && (e036>e258  || e036==e258))
00519             return     E22=themap[0]+themap[1]+themap[3]+themap[4];
00520     else if ( (e012>e678 || e012==e678)  && (e036<e258 || e036==e258) )
00521             return    E22=themap[1]+themap[2]+themap[4]+themap[5];
00522     else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
00523             return     E22=themap[3]+themap[4]+themap[6]+themap[7];
00524     else if ( (e012<e678|| e012==e678)  && (e036<e258|| e036==e258) )
00525             return    E22=themap[4]+themap[5]+themap[7]+themap[8];
00526     else {
00527             return E22;
00528     }
00529 }
00530 
00531 float EcalSimHitsValidProducer::eCluster4x4(float e33,  MapType&  themap){
00532         float E44=0.;
00533         float e0_4   = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
00534         float e0_20  = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
00535         float e4_24  = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
00536         float e0_24  = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];
00537 
00538    if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
00539            return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
00540    else if ((e0_4>e0_24 || e0_4==e0_24)  && (e0_20<e4_24 || e0_20==e4_24))
00541            return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
00542    else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
00543            return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
00544    else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
00545            return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
00546    else{
00547            return E44;
00548    }
00549 }
00550 
00551 float EcalSimHitsValidProducer::energyInEEMatrix(int nCellInX, int nCellInY,
00552                                         int centralX, int centralY,
00553                                         int centralZ, MapType& themap){
00554 
00555   int   ncristals   = 0;
00556   float totalEnergy = 0.;
00557 
00558   int goBackInX = nCellInX/2;
00559   int goBackInY = nCellInY/2;
00560   int startX    = centralX-goBackInX;
00561   int startY    = centralY-goBackInY;
00562 
00563   for (int ix=startX; ix<startX+nCellInX; ix++) {
00564     for (int iy=startY; iy<startY+nCellInY; iy++) {
00565 
00566       uint32_t index ;
00567       if (EEDetId::validDetId(ix, iy,centralZ)) {
00568         index = EEDetId(ix,iy,centralZ).rawId();
00569       } else { continue; }
00570 
00571       totalEnergy   += themap[index];
00572       ncristals     += 1;
00573 
00574       LogDebug("EcalSimHitsValidProducer")   
00575       << " EnergyInEEMatrix: ix - iy - E = " << ix
00576       << "  " << iy << " "  << themap[index] ;
00577     
00578     }
00579   }
00580 
00581     LogDebug("EcalSimHitsValidProducer")
00582     << " EnergyInEEMatrix: energy in " << nCellInX
00583     << " cells in x times " << nCellInY
00584     << " cells in y matrix = " << totalEnergy
00585     << " for " << ncristals << " crystals";
00586 
00587   return totalEnergy;
00588 
00589 }
00590 
00591 float EcalSimHitsValidProducer::energyInEBMatrix(int nCellInEta, int nCellInPhi,
00592                                         int centralEta, int centralPhi,
00593                                         int centralZ, MapType& themap){
00594 
00595   int   ncristals   = 0;
00596   float totalEnergy = 0.;
00597 
00598   int goBackInEta = nCellInEta/2;
00599   int goBackInPhi = nCellInPhi/2;
00600   int startEta    = centralZ*centralEta-goBackInEta;
00601   int startPhi    = centralPhi-goBackInPhi;
00602 
00603   for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
00604     for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {
00605 
00606       uint32_t index ;
00607       if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
00608       if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
00609       else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
00610       else              { index = EBDetId(ieta,iphi).rawId();     }
00611 
00612       totalEnergy   += themap[index];
00613       ncristals     += 1;
00614 
00615       LogDebug("EcalSimHitsValidProducer")
00616       << " EnergyInEBMatrix: ieta - iphi - E = " << ieta
00617       << "  " << iphi << " "  << themap[index];
00618 
00619     }
00620   }
00621 
00622     LogDebug("EcalSimHitsValidProducer")
00623     << " EnergyInEBMatrix: energy in " << nCellInEta
00624     << " cells in eta times " << nCellInPhi
00625     << " cells in phi matrix = " << totalEnergy
00626     << " for " << ncristals << " crystals";
00627 
00628   return totalEnergy;
00629 
00630 }
00631 
00632 uint32_t EcalSimHitsValidProducer::getUnitWithMaxEnergy(MapType& themap) {
00633 
00634   uint32_t unitWithMaxEnergy = 0;
00635   float    maxEnergy = 0.;
00636 
00637   MapType::iterator iter;
00638   for (iter = themap.begin(); iter != themap.end(); iter++) {
00639 
00640     if (maxEnergy < (*iter).second) {
00641       maxEnergy = (*iter).second;
00642       unitWithMaxEnergy = (*iter).first;
00643     }
00644   }
00645 
00646   
00647    LogDebug("EcalSimHitsValidProducer")
00648    << " Max energy of " << maxEnergy
00649    << " MeV was found in Unit id 0x" << std::hex
00650     << unitWithMaxEnergy << std::dec; 
00651 
00652   return unitWithMaxEnergy;
00653 }
00654