#include <EcalSimHitsValidProducer.h>

Inheritance diagram for EcalSimHitsValidProducer:

Public Member Functions
	EcalSimHitsValidProducer (const edm::ParameterSet &)
void	produce (edm::Event &, const edm::EventSetup &)
virtual	~EcalSimHitsValidProducer ()
Private Types
typedef std::vector< float >	FloatVector
typedef std::map< uint32_t, float, std::less< uint32_t > >	MapType
Private Member Functions
	EcalSimHitsValidProducer (const EcalSimHitsValidProducer &)
float	eCluster2x2 (MapType &themap)
float	eCluster4x4 (float e33, MapType &themap)
float	energyInEBMatrix (int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType &themap)
float	energyInEEMatrix (int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType &themap)
bool	fillEBMatrix (int nCellInEta, int nCellInPhi, int CentralEta, int CentralPhi, int CentralZ, MapType &fillmap, MapType &themap)
bool	fillEEMatrix (int nCellInEta, int nCellInPhi, int CentralEta, int CentralPhi, int CentralZ, MapType &fillmap, MapType &themap)
void	fillEventInfo (PEcalValidInfo &)
uint32_t	getUnitWithMaxEnergy (MapType &themap)
const EcalSimHitsValidProducer &	operator= (const EcalSimHitsValidProducer &)
void	update (const G4Step *)
	This routine will be called when the appropriate signal arrives.
void	update (const BeginOfEvent *)
	This routine will be called when the appropriate signal arrives.
void	update (const EndOfEvent *)
	This routine will be called when the appropriate signal arrives.
Private Attributes
float	eb1
float	eb16
float	eb25
float	eb4
float	eb9
float	eBX0 [26]
float	ee1
float	ee16
float	ee25
float	ee4
float	ee9
float	eEX0 [26]
FloatVector	eOf1ES
FloatVector	eOf1ESzm
FloatVector	eOf1ESzp
FloatVector	eOf2ES
FloatVector	eOf2ESzm
FloatVector	eOf2ESzp
FloatVector	eOfEBCaloG4Hit
FloatVector	eOfEECaloG4Hit
FloatVector	eOfEEMinusCaloG4Hit
FloatVector	eOfEEPlusCaloG4Hit
FloatVector	eOfESCaloG4Hit
FloatVector	etaOfEBCaloG4Hit
FloatVector	etaOfEECaloG4Hit
FloatVector	etaOfESCaloG4Hit
std::string	label
int	nCrystalInEB
int	nCrystalInEEzm
int	nCrystalInEEzp
int	nHitsIn1ES
int	nHitsIn1ESzm
int	nHitsIn1ESzp
int	nHitsIn2ES
int	nHitsIn2ESzm
int	nHitsIn2ESzp
int	nHitsInEB
int	nHitsInEE
int	nHitsInES
FloatVector	phiOfEBCaloG4Hit
FloatVector	phiOfEECaloG4Hit
FloatVector	phiOfESCaloG4Hit
math::XYZTLorentzVector	theMomentum
int	thePID
math::XYZTLorentzVector	theVertex
FloatVector	tOfEBCaloG4Hit
FloatVector	tOfEECaloG4Hit
FloatVector	tOfESCaloG4Hit
float	totalEInEB
float	totalEInEE
float	totalEInEEzm
float	totalEInEEzp
float	totalEInES
float	totalEInESzm
float	totalEInESzp
int	totalHits
FloatVector	zOfES

Detailed Description

Definition at line 21 of file EcalSimHitsValidProducer.h.

Member Typedef Documentation

typedef std::vector<float> EcalSimHitsValidProducer::FloatVector [private]

Definition at line 26 of file EcalSimHitsValidProducer.h.

typedef std::map<uint32_t,float,std::less<uint32_t> > EcalSimHitsValidProducer::MapType [private]

Definition at line 27 of file EcalSimHitsValidProducer.h.

Constructor & Destructor Documentation

EcalSimHitsValidProducer::EcalSimHitsValidProducer ( const edm::ParameterSet & iPSet )

Definition at line 20 of file EcalSimHitsValidProducer.cc.

References eBX0, eEX0, i, and label.

                                                                               :
     ee1(0.0),ee4(0.0),ee9(0.0),ee16(0.0),ee25(0.0),
     eb1(0.0),eb4(0.0),eb9(0.0),eb16(0.0),eb25(0.0),
     totalEInEE(0.0),totalEInEB(0),totalEInES(0.0),totalEInEEzp(0.0),totalEInEEzm(0.0),
     totalEInESzp(0.0),totalEInESzm(0.0), 
     totalHits(0),nHitsInEE(0),nHitsInEB(0),nHitsInES(0),nHitsIn1ES(0),nHitsIn2ES(0),
     nCrystalInEB(0),nCrystalInEEzp(0),nCrystalInEEzm(0),
     nHitsIn1ESzp(0),nHitsIn1ESzm(0),nHitsIn2ESzp(0),nHitsIn2ESzm(0),
     thePID(0),
     label(iPSet.getUntrackedParameter<std::string>("instanceLabel","EcalValidInfo") ) 
{
   produces<PEcalValidInfo>(label);
  
  for ( int i = 0; i<26; i++ ) {
     eBX0[i] = 0.0;
     eEX0[i] = 0.0;
  }
  
}

EcalSimHitsValidProducer::~EcalSimHitsValidProducer ( ) [virtual]

Definition at line 41 of file EcalSimHitsValidProducer.cc.

                                                   {

}

EcalSimHitsValidProducer::EcalSimHitsValidProducer ( const EcalSimHitsValidProducer & ) [private]

Member Function Documentation

float EcalSimHitsValidProducer::eCluster2x2 ( MapType & themap ) [private]

Definition at line 512 of file EcalSimHitsValidProducer.cc.

Referenced by update().

                                                           {
        float  E22=0.;
        float e012  = themap[0]+themap[1]+themap[2];
        float e036  = themap[0]+themap[3]+themap[6];
        float e678  = themap[6]+themap[7]+themap[8];
        float e258  = themap[2]+themap[5]+themap[8];

    if ( (e012>e678 || e012==e678) && (e036>e258  || e036==e258))
            return     E22=themap[0]+themap[1]+themap[3]+themap[4];
    else if ( (e012>e678 || e012==e678)  && (e036<e258 || e036==e258) )
            return    E22=themap[1]+themap[2]+themap[4]+themap[5];
    else if ( (e012<e678 || e012==e678) && (e036>e258 || e036==e258))
            return     E22=themap[3]+themap[4]+themap[6]+themap[7];
    else if ( (e012<e678|| e012==e678)  && (e036<e258|| e036==e258) )
            return    E22=themap[4]+themap[5]+themap[7]+themap[8];
    else {
            return E22;
    }
}

float EcalSimHitsValidProducer::eCluster4x4	(	float	e33,
		MapType &	themap
	)		`[private]`

Definition at line 532 of file EcalSimHitsValidProducer.cc.

Referenced by update().

                                                                       {
        float E44=0.;
        float e0_4   = themap[0]+themap[1]+themap[2]+themap[3]+themap[4];
        float e0_20  = themap[0]+themap[5]+themap[10]+themap[15]+themap[20];
        float e4_24  = themap[4]+themap[9]+themap[14]+themap[19]+themap[24];
        float e0_24  = themap[20]+themap[21]+themap[22]+themap[23]+themap[24];

   if ((e0_4>e0_24 || e0_4==e0_24) && (e0_20>e4_24|| e0_20==e4_24))
           return E44=e33+themap[0]+themap[1]+themap[2]+themap[3]+themap[5]+themap[10]+themap[15];
   else if ((e0_4>e0_24 || e0_4==e0_24)  && (e0_20<e4_24 || e0_20==e4_24))
           return E44=e33+themap[1]+themap[2]+themap[3]+themap[4]+themap[9]+themap[14]+themap[19];
   else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20>e4_24 || e0_20==e4_24))
           return E44=e33+themap[5]+themap[10]+themap[15]+themap[20]+themap[21]+themap[22]+themap[23];
   else if ((e0_4<e0_24|| e0_4==e0_24) && (e0_20<e4_24 || e0_20==e4_24))
           return E44=e33+themap[21]+themap[22]+themap[23]+themap[24]+themap[9]+themap[14]+themap[19];
   else{
           return E44;
   }
}

float EcalSimHitsValidProducer::energyInEBMatrix	(	int	nCellInX,
		int	nCellInY,
		int	centralX,
		int	centralY,
		int	centralZ,
		MapType &	themap
	)		`[private]`

Definition at line 592 of file EcalSimHitsValidProducer.cc.

References abs, getHLTprescales::index, LogDebug, and DetId::rawId().

Referenced by update().

                                                                      {

  int   ncristals   = 0;
  float totalEnergy = 0.;

  int goBackInEta = nCellInEta/2;
  int goBackInPhi = nCellInPhi/2;
  int startEta    = centralZ*centralEta-goBackInEta;
  int startPhi    = centralPhi-goBackInPhi;

  for (int ieta=startEta; ieta<startEta+nCellInEta; ieta++) {
    for (int iphi=startPhi; iphi<startPhi+nCellInPhi; iphi++) {

      uint32_t index ;
      if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
      if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
      else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
      else              { index = EBDetId(ieta,iphi).rawId();     }

      totalEnergy   += themap[index];
      ncristals     += 1;

      LogDebug("EcalSimHitsValidProducer")
      << " EnergyInEBMatrix: ieta - iphi - E = " << ieta
      << "  " << iphi << " "  << themap[index];

    }
  }

    LogDebug("EcalSimHitsValidProducer")
    << " EnergyInEBMatrix: energy in " << nCellInEta
    << " cells in eta times " << nCellInPhi
    << " cells in phi matrix = " << totalEnergy
    << " for " << ncristals << " crystals";

  return totalEnergy;

}

float EcalSimHitsValidProducer::energyInEEMatrix	(	int	nCellInX,
		int	nCellInY,
		int	centralX,
		int	centralY,
		int	centralZ,
		MapType &	themap
	)		`[private]`

Definition at line 552 of file EcalSimHitsValidProducer.cc.

References getHLTprescales::index, LogDebug, DetId::rawId(), and EEDetId::validDetId().

Referenced by update().

                                                                      {

  int   ncristals   = 0;
  float totalEnergy = 0.;

  int goBackInX = nCellInX/2;
  int goBackInY = nCellInY/2;
  int startX    = centralX-goBackInX;
  int startY    = centralY-goBackInY;

  for (int ix=startX; ix<startX+nCellInX; ix++) {
    for (int iy=startY; iy<startY+nCellInY; iy++) {

      uint32_t index ;
      if (EEDetId::validDetId(ix, iy,centralZ)) {
        index = EEDetId(ix,iy,centralZ).rawId();
      } else { continue; }

      totalEnergy   += themap[index];
      ncristals     += 1;

      LogDebug("EcalSimHitsValidProducer")   
      << " EnergyInEEMatrix: ix - iy - E = " << ix
      << "  " << iy << " "  << themap[index] ;
    
    }
  }

    LogDebug("EcalSimHitsValidProducer")
    << " EnergyInEEMatrix: energy in " << nCellInX
    << " cells in x times " << nCellInY
    << " cells in y matrix = " << totalEnergy
    << " for " << ncristals << " crystals";

  return totalEnergy;

}

bool EcalSimHitsValidProducer::fillEBMatrix	(	int	nCellInEta,
		int	nCellInPhi,
		int	CentralEta,
		int	CentralPhi,
		int	CentralZ,
		MapType &	fillmap,
		MapType &	themap
	)		`[private]`

Definition at line 481 of file EcalSimHitsValidProducer.cc.

References abs, getUnitWithMaxEnergy(), i, getHLTprescales::index, and DetId::rawId().

Referenced by update().

{
   int goBackInEta = nCellInEta/2;
   int goBackInPhi = nCellInPhi/2;

   int startEta  =  CentralZ*CentralEta - goBackInEta;
   int startPhi  =  CentralPhi - goBackInPhi;

   int i = 0 ;
   for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {
      for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
          uint32_t  index;
          if (abs(ieta) > 85 || abs(ieta)<1 ) { continue; }
          if (iphi< 1)      { index = EBDetId(ieta,iphi+360).rawId(); }
          else if(iphi>360) { index = EBDetId(ieta,iphi-360).rawId(); }
          else              { index = EBDetId(ieta,iphi).rawId();     }
          fillmap[i++] = themap[index];
      }
   }

   uint32_t  ebcenterid = getUnitWithMaxEnergy(themap);

   if ( fillmap[i/2] == themap[ebcenterid] )
        return true;
   else
        return false;
}

bool EcalSimHitsValidProducer::fillEEMatrix	(	int	nCellInEta,
		int	nCellInPhi,
		int	CentralEta,
		int	CentralPhi,
		int	CentralZ,
		MapType &	fillmap,
		MapType &	themap
	)		`[private]`

Definition at line 450 of file EcalSimHitsValidProducer.cc.

References getUnitWithMaxEnergy(), i, getHLTprescales::index, DetId::rawId(), and EEDetId::validDetId().

Referenced by update().

{
   int goBackInEta = nCellInEta/2;
   int goBackInPhi = nCellInPhi/2;

   int startEta  =  CentralEta - goBackInEta;
   int startPhi  =  CentralPhi - goBackInPhi;

   int i = 0 ;
   for ( int ieta = startEta; ieta < startEta+nCellInEta; ieta ++ ) {

      for( int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++ ) {
        uint32_t index;

        if (EEDetId::validDetId(ieta, iphi,CentralZ)) {
          index = EEDetId( ieta, iphi,CentralZ).rawId();
        } else { continue; }
        
        fillmap[i++] = themap[index];
      }
   }
   uint32_t  centerid = getUnitWithMaxEnergy(themap);

   if ( fillmap[i/2] == themap[centerid] ) 
        return true;
   else
        return false;
}

void EcalSimHitsValidProducer::fillEventInfo ( PEcalValidInfo & product ) [private]

Definition at line 54 of file EcalSimHitsValidProducer.cc.

Referenced by produce().

{
 if( ee1 != 0 ){
   product.ee1 = ee1;
   product.ee4 = ee4;
   product.ee9 = ee9;
   product.ee16 = ee16;
   product.ee25 = ee25;
   for ( int i = 0; i<26; i++ ) {
      product.eEX0.push_back( eEX0[i]);
   }
 } 

 if( eb1 != 0 ){
   product.eb1 = eb1;
   product.eb4 = eb4;
   product.eb9 = eb9;
   product.eb16 = eb16;
   product.eb25 = eb25;
   for ( int i = 0; i<26; i++ ) {
      product.eBX0.push_back( eBX0[i]);
   }
 }

   product.totalEInEE  = totalEInEE;
   product.totalEInEB  = totalEInEB;
   product.totalEInES  = totalEInES;

   product.totalEInEEzp  = totalEInEEzp;
   product.totalEInEEzm  = totalEInEEzm;

   product.totalEInESzp  = totalEInESzp;
   product.totalEInESzm  = totalEInESzm;

   product.totalHits  = totalHits;
   product.nHitsInEE  = nHitsInEE;
   product.nHitsInEB  = nHitsInEB;
   product.nHitsInES  = nHitsInES;
   product.nHitsIn1ES = nHitsIn1ES;
   product.nHitsIn2ES = nHitsIn2ES;
   product.nCrystalInEB   = nCrystalInEB;
   product.nCrystalInEEzp = nCrystalInEEzp;
   product.nCrystalInEEzm = nCrystalInEEzm; 

   product.nHitsIn1ESzp = nHitsIn1ESzp;
   product.nHitsIn1ESzm = nHitsIn1ESzm;
   product.nHitsIn2ESzp = nHitsIn2ESzp;
   product.nHitsIn2ESzm = nHitsIn2ESzm;


   product.eOf1ES = eOf1ES;
   product.eOf2ES = eOf2ES;
   product.zOfES  = zOfES;

   product.eOf1ESzp = eOf1ESzp;
   product.eOf1ESzm = eOf1ESzm;
   product.eOf2ESzp = eOf2ESzp;
   product.eOf2ESzm = eOf2ESzm;


   product.phiOfEECaloG4Hit = phiOfEECaloG4Hit;
   product.etaOfEECaloG4Hit = etaOfEECaloG4Hit;
   product.eOfEECaloG4Hit   = eOfEECaloG4Hit;
   product.eOfEEPlusCaloG4Hit    = eOfEEPlusCaloG4Hit;
   product.eOfEEMinusCaloG4Hit   = eOfEEMinusCaloG4Hit;
   product.tOfEECaloG4Hit        = tOfEECaloG4Hit;

   product.phiOfESCaloG4Hit = phiOfESCaloG4Hit;
   product.etaOfESCaloG4Hit = etaOfESCaloG4Hit;
   product.eOfESCaloG4Hit   = eOfESCaloG4Hit;
   product.tOfESCaloG4Hit   = tOfESCaloG4Hit;

   product.phiOfEBCaloG4Hit = phiOfEBCaloG4Hit;
   product.etaOfEBCaloG4Hit = etaOfEBCaloG4Hit;
   product.eOfEBCaloG4Hit   = eOfEBCaloG4Hit;
   product.tOfEBCaloG4Hit   = tOfEBCaloG4Hit;


   product.theMomentum = theMomentum;
   product.theVertex   = theVertex;
   product.thePID      = thePID;
}

uint32_t EcalSimHitsValidProducer::getUnitWithMaxEnergy ( MapType & themap ) [private]

Definition at line 633 of file EcalSimHitsValidProducer.cc.

References LogDebug.

Referenced by fillEBMatrix(), fillEEMatrix(), and update().

                                                                       {

  uint32_t unitWithMaxEnergy = 0;
  float    maxEnergy = 0.;

  MapType::iterator iter;
  for (iter = themap.begin(); iter != themap.end(); iter++) {

    if (maxEnergy < (*iter).second) {
      maxEnergy = (*iter).second;
      unitWithMaxEnergy = (*iter).first;
    }
  }

  
   LogDebug("EcalSimHitsValidProducer")
   << " Max energy of " << maxEnergy
   << " MeV was found in Unit id 0x" << std::hex
    << unitWithMaxEnergy << std::dec; 

  return unitWithMaxEnergy;
}

const EcalSimHitsValidProducer& EcalSimHitsValidProducer::operator= ( const EcalSimHitsValidProducer & ) [private]

void EcalSimHitsValidProducer::produce	(	edm::Event &	e,
		const edm::EventSetup &
	)		`[virtual]`

Implements SimProducer.

Definition at line 46 of file EcalSimHitsValidProducer.cc.

References fillEventInfo(), label, and edm::Event::put().

{
   std::auto_ptr<PEcalValidInfo> product(new PEcalValidInfo );
   fillEventInfo(*product);
   e.put(product,label);
}

void EcalSimHitsValidProducer::update ( const G4Step * ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const G4Step * >.

Definition at line 422 of file EcalSimHitsValidProducer.cc.

References eBX0, eEX0, mergeVDriftHistosByStation::name, csvReporter::r, mathSSE::sqrt(), x, detailsBasic3DVector::y, and z.

                                                   {
 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
  G4ThreeVector hitPoint = preStepPoint->GetPosition();
  G4VPhysicalVolume* currentPV  = preStepPoint->GetPhysicalVolume();
  G4String name         = currentPV->GetName();
  std::string crystal;
  crystal.assign(name,0,4);

  float Edeposit =  aStep->GetTotalEnergyDeposit();
  if ( crystal == "EFRY"&& Edeposit > 0.0){
          float z = hitPoint.z();
          float detz = fabs(fabs(z)-3200);
          int x0 = (int)floor( detz/8.9 );
          if ( x0< 26){ 
            eEX0[x0] += Edeposit;
          }
  } 
  if(crystal == "EBRY" && Edeposit > 0.0) {
            float x = hitPoint.x();
            float y = hitPoint.y();
            float r = sqrt(x*x +y*y);
            float detr = r -1290;
            int x0 = (int)floor( detr/8.9);
            eBX0[x0] += Edeposit;
  }
   
}

void EcalSimHitsValidProducer::update ( const BeginOfEvent * ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfEvent * >.

Definition at line 138 of file EcalSimHitsValidProducer.cc.

References eb1, eb16, eb25, eb4, eb9, eBX0, ee1, ee16, ee25, ee4, ee9, eEX0, eOf1ES, eOf1ESzm, eOf1ESzp, eOf2ES, eOf2ESzm, eOf2ESzp, eOfEBCaloG4Hit, eOfEECaloG4Hit, eOfEEMinusCaloG4Hit, eOfEEPlusCaloG4Hit, eOfESCaloG4Hit, etaOfEBCaloG4Hit, etaOfEECaloG4Hit, etaOfESCaloG4Hit, i, nCrystalInEB, nCrystalInEEzm, nCrystalInEEzp, nHitsIn1ES, nHitsIn1ESzm, nHitsIn1ESzp, nHitsIn2ES, nHitsIn2ESzm, nHitsIn2ESzp, nHitsInEB, nHitsInEE, nHitsInES, phiOfEBCaloG4Hit, phiOfEECaloG4Hit, phiOfESCaloG4Hit, tOfEBCaloG4Hit, tOfEECaloG4Hit, tOfESCaloG4Hit, totalEInEB, totalEInEE, totalEInEEzm, totalEInEEzp, totalEInES, totalEInESzm, totalEInESzp, totalHits, and zOfES.

                                                   {
  ee1  = 0.0;
  ee4  = 0.0;
  ee9  = 0.0;
  ee16 = 0.0;
  ee25 = 0.0;

  eb1  = 0.0;
  eb4  = 0.0;
  eb9  = 0.0;
  eb16 = 0.0;
  eb25 = 0.0;

  totalEInEE = 0.0 ;
  totalEInEB = 0.0 ;
  totalEInES = 0.0 ;
  
  totalEInEEzp = 0.0 ;
  totalEInEEzm = 0.0 ;
  totalEInESzp = 0.0 ;
  totalEInESzm = 0.0 ;

  totalHits  = 0 ;
  nHitsInEE  = 0 ;
  nHitsInEB  = 0 ;
  nHitsInES  = 0 ;
  nHitsIn1ES = 0 ;
  nHitsIn2ES = 0 ;
  nCrystalInEB   = 0;
  nCrystalInEEzp = 0;
  nCrystalInEEzm = 0;
 

  nHitsIn1ESzp = 0 ;
  nHitsIn1ESzm = 0 ;
  nHitsIn2ESzp = 0 ;
  nHitsIn2ESzm = 0 ;
 
  for ( int i = 0; i<26; i++ ) {
     eBX0[i] = 0.0;
     eEX0[i] = 0.0;
  }

  eOf1ES.clear();
  eOf2ES.clear();
  zOfES.clear();

  eOf1ESzp.clear();
  eOf1ESzm.clear();
  eOf2ESzp.clear();
  eOf2ESzm.clear();

  phiOfEECaloG4Hit.clear();
  etaOfEECaloG4Hit.clear();
  tOfEECaloG4Hit.clear();
  eOfEECaloG4Hit.clear();
  eOfEEPlusCaloG4Hit.clear();
  eOfEEMinusCaloG4Hit.clear();

  phiOfESCaloG4Hit.clear();
  etaOfESCaloG4Hit.clear();
  tOfESCaloG4Hit.clear();
  eOfESCaloG4Hit.clear();

  phiOfEBCaloG4Hit.clear();
  etaOfEBCaloG4Hit.clear();
  tOfEBCaloG4Hit.clear();
  eOfEBCaloG4Hit.clear();
}

void EcalSimHitsValidProducer::update ( const EndOfEvent * ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfEvent * >.

Definition at line 209 of file EcalSimHitsValidProducer.cc.

References eb1, eb16, eb25, eb4, eb9, eCluster2x2(), eCluster4x4(), ee1, ee16, ee25, ee4, ee9, energyInEBMatrix(), energyInEEMatrix(), eOf1ESzm, eOf1ESzp, eOf2ESzm, eOf2ESzp, eOfEBCaloG4Hit, eOfEECaloG4Hit, eOfEEMinusCaloG4Hit, eOfEEPlusCaloG4Hit, etaOfEBCaloG4Hit, etaOfEECaloG4Hit, fillEBMatrix(), fillEEMatrix(), CaloG4Hit::getEnergyDeposit(), CaloG4Hit::getEntry(), CaloG4Hit::getTimeSlice(), CaloG4Hit::getUnitID(), getUnitWithMaxEnergy(), i, EBDetId::ietaAbs(), EBDetId::iphi(), EEDetId::ix(), EEDetId::iy(), j, funct::log(), max(), min, nCrystalInEB, nCrystalInEEzm, nCrystalInEEzp, nHitsIn1ESzm, nHitsIn1ESzp, nHitsIn2ESzm, nHitsIn2ESzp, nHitsInEB, nHitsInEE, nHitsInES, npart, phiOfEBCaloG4Hit, phiOfEECaloG4Hit, ESDetId::plane(), funct::pow(), mathSSE::sqrt(), funct::tan(), theMomentum, thePID, theta(), theVertex, tOfEBCaloG4Hit, tOfEECaloG4Hit, totalEInEB, totalEInEE, totalEInEEzm, totalEInEEzp, totalEInES, totalEInESzm, totalEInESzp, totalHits, ESDetId::zside(), EEDetId::zside(), and EBDetId::zside().

                                                     {

  int trackID = 0;
  G4PrimaryParticle * thePrim = 0;
  int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
  if ( nvertex <= 0) {
          edm::LogWarning("EcalSimHitsValidProducer")
                <<" No Vertex in this Event!";
  }else {
    for ( int i = 0; i< nvertex; i++){
          G4PrimaryVertex * avertex =(*evt)()->GetPrimaryVertex(i);
          if ( avertex == 0 )
                  edm::LogWarning("EcalSimHitsValidProducer")
                  <<" Pointer to vertex is NULL!";
          else {
             float x0 = avertex->GetX0();
             float y0 = avertex->GetY0();
             float z0 = avertex->GetZ0();
             float t0 = avertex->GetT0();
             theVertex.SetCoordinates(x0,y0,z0,t0);
    
             int npart = avertex->GetNumberOfParticle();
             if ( npart == 0)
                edm::LogWarning("EcalSimHitsValidProducer")
                <<" No primary particle in this event";
             else {
                 if ( thePrim == 0)
                   thePrim = avertex->GetPrimary(trackID);
             }
          }
     }

   // the direction of momentum of primary particles
    double etaInit =0, phiInit =0, pInit =0;
     if ( thePrim != 0){
          double  px = thePrim -> GetPx();
          double  py = thePrim -> GetPy();
          double  pz = thePrim -> GetPz();
          theMomentum.SetCoordinates(px,py,pz,0.);

          pInit =sqrt( pow(px,2.) + pow(py,2.) + pow(pz,2.));
          if ( pInit == 0)
                  edm::LogWarning("EcalSimHitsValidProducer") 
                  <<" Primary has p = 0 ; ";
          else {
                  theMomentum.SetE(pInit);
                  double costheta  = pz/pInit;
                  double theta = acos(std::min(std::max(costheta, -1.),1.));
                  etaInit = -log(tan(theta/2));

                  if ( px != 0 || py != 0)
                          phiInit = atan2(py,px);
          }

          thePID = thePrim->GetPDGcode();
      }else {
          edm::LogWarning("EcalSimHitsValidProducer")
          <<" Could not find the primary particle!!";
      }
  }
  // hit map for EB for matrices
  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
  int EBHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsEB");
  int EEHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsEE");
  int SEHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsES");
 

  CaloG4HitCollection* theEBHC = (CaloG4HitCollection*) allHC->GetHC(EBHCid);
  CaloG4HitCollection* theEEHC = (CaloG4HitCollection*) allHC->GetHC(EEHCid);
  CaloG4HitCollection* theSEHC = (CaloG4HitCollection*) allHC->GetHC(SEHCid);
  
  nHitsInEE = theEEHC->entries();
  nHitsInEB = theEBHC->entries();
  nHitsInES = theSEHC->entries();
  totalHits = nHitsInEE + nHitsInEB + nHitsInES;

 //   EB Hit collection start
 MapType ebmap;
 for (int j=0; j<theEBHC->entries(); j++) {
   CaloG4Hit* aHit = (*theEBHC)[j];
    totalEInEB += aHit->getEnergyDeposit();
    float he     = aHit->getEnergyDeposit();
    float htime  = aHit->getTimeSlice();
                                                                                                                             
    math::XYZPoint hpos = aHit->getEntry();
    float htheta    = hpos.theta();
    float heta    = -log(tan(htheta * 0.5));
    float hphi    = hpos.phi();

    phiOfEBCaloG4Hit.push_back( hphi);
    etaOfEBCaloG4Hit.push_back(heta);
    tOfEBCaloG4Hit.push_back(htime);
    eOfEBCaloG4Hit.push_back(he);
    uint32_t  crystid = aHit->getUnitID();
    ebmap[crystid] += aHit->getEnergyDeposit();
}

    nCrystalInEB = ebmap.size();

 //   EE Hit collection start
 MapType eemap,eezpmap,eezmmap;
 
 for (int j=0; j<theEEHC->entries(); j++) {
    CaloG4Hit* aHit = (*theEEHC)[j];
    totalEInEE += aHit->getEnergyDeposit();
    float he     = aHit->getEnergyDeposit();
    float htime  = aHit->getTimeSlice();

    math::XYZPoint hpos = aHit->getEntry();
    float htheta    = hpos.theta();
    float heta    = -log(tan(htheta * 0.5));
    float hphi    = hpos.phi();
    phiOfEECaloG4Hit.push_back( hphi);
    etaOfEECaloG4Hit.push_back(heta);
    tOfEECaloG4Hit.push_back(htime);
    eOfEECaloG4Hit.push_back(he);
 
    uint32_t  crystid = aHit->getUnitID();
    EEDetId myEEid(crystid);
    if ( myEEid.zside() == -1 ) {
          totalEInEEzm += aHit->getEnergyDeposit();
          eOfEEMinusCaloG4Hit.push_back(he);
          eezmmap[crystid] += aHit->getEnergyDeposit();
     }
    if ( myEEid.zside() == 1 ) {
          totalEInEEzp += aHit->getEnergyDeposit();
          eOfEEPlusCaloG4Hit.push_back(he);
          eezpmap[crystid] += aHit->getEnergyDeposit();
     }
               

    eemap[crystid] += aHit->getEnergyDeposit(); 
 }

 nCrystalInEEzm = eezmmap.size();
 nCrystalInEEzp = eezpmap.size(); 

 //Hits from ES
 for (int j=0; j<theSEHC->entries(); j++) {
    CaloG4Hit* aHit = (*theSEHC)[j];
    totalEInES += aHit->getEnergyDeposit();
    ESDetId esid = ESDetId( aHit->getUnitID());

    if (esid.zside() == -1) {

           totalEInESzm +=  aHit->getEnergyDeposit();

           if (esid.plane() ==  1) {
                 nHitsIn1ESzm++;
                 eOf1ESzm.push_back(aHit->getEnergyDeposit());
            }else if ( esid.plane() ==  2){
                  nHitsIn2ESzm++;
                  eOf2ESzm.push_back(aHit->getEnergyDeposit());
            }
    } 
    if (esid.zside() == 1) {

           totalEInESzp +=  aHit->getEnergyDeposit();

           if (esid.plane() ==  1) {
                 nHitsIn1ESzp++;
                 eOf1ESzp.push_back(aHit->getEnergyDeposit());
            }else if ( esid.plane() ==  2){
                  nHitsIn2ESzp++;
                  eOf2ESzp.push_back(aHit->getEnergyDeposit());
            }
    }


  
  }
  
   uint32_t  eemaxid = getUnitWithMaxEnergy(eemap);
   uint32_t  ebmaxid = getUnitWithMaxEnergy(ebmap); 
   if ( eemap[eemaxid] > ebmap[ebmaxid] ) {
    uint32_t  centerid = getUnitWithMaxEnergy(eemap);
    EEDetId myEEid(centerid);
    int ix = myEEid.ix();
    int iy = myEEid.iy();
    int iz = myEEid.zside();

    ee1 =  energyInEEMatrix(1,1,ix,iy,iz,eemap);
    ee9 =  energyInEEMatrix(3,3,ix,iy,iz,eemap);
    ee25=  energyInEEMatrix(5,5,ix,iy,iz,eemap);
    MapType   neweemap;
    if( fillEEMatrix(3,3,ix,iy,iz,neweemap, eemap)) {
          ee4 = eCluster2x2(neweemap);
    }
    if ( fillEEMatrix(5,5,ix,iy,iz,neweemap, eemap)) {
          ee16 = eCluster4x4(ee9,neweemap);
    }
  } else {
    uint32_t  ebcenterid = getUnitWithMaxEnergy(ebmap);
    EBDetId myEBid(ebcenterid);
    int bx = myEBid.ietaAbs();
    int by = myEBid.iphi();
    int bz = myEBid.zside();
    eb1 =  energyInEBMatrix(1,1,bx,by,bz,ebmap);
    eb9 =  energyInEBMatrix(3,3,bx,by,bz,ebmap);
    eb25=  energyInEBMatrix(5,5,bx,by,bz,ebmap);

    MapType  newebmap;
    if( fillEBMatrix(3,3,bx,by,bz,newebmap, ebmap)){
          eb4 = eCluster2x2(newebmap);
    }
    if( fillEBMatrix(5,5,bx,by,bz,newebmap, ebmap)){
          eb16 = eCluster4x4(eb9,newebmap);
    }
  } 

}