#include <ContainmentCorrectionAnalyzer.h>

Inheritance diagram for ContainmentCorrectionAnalyzer:

Public Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &)

virtual void	beginJob ()

	ContainmentCorrectionAnalyzer (const edm::ParameterSet &)

virtual void	endJob ()

	~ContainmentCorrectionAnalyzer ()

Public Member Functions inherited from edm::EDAnalyzer
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

	EDAnalyzer ()

ModuleDescription const &	moduleDescription () const

std::string	workerType () const

virtual	~EDAnalyzer ()

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

virtual	~EDConsumerBase () noexcept(false)

Private Member Functions
float	ecalEta (float EtaParticle, float Zvertex, float plane_Radius)

void	fillMcTruth (std::vector< SimTrack > &theSimTracks, std::vector< SimVertex > &theSimVertices)

std::vector< EcalSimPhotonMCTruth >	findMcTruth (std::vector< SimTrack > &theSimTracks, std::vector< SimVertex > &theSimVertices)

Private Attributes
edm::EDGetTokenT< reco::SuperClusterCollection >	BarrelSuperClusterCollection_

std::vector< float >	e1

std::vector< float >	e25

std::vector< float >	e9

edm::EDGetTokenT< reco::SuperClusterCollection >	EndcapSuperClusterCollection_

std::map< unsigned, unsigned >	geantToIndex_

TH1F *	h_EB_converted

TH1F *	h_EB_e25EtrueReference

TH1F *	h_EB_e9EtrueReference

TH1F *	h_EB_eRecoEtrueReference

TH1F *	h_EB_eTrue

TH1F *	h_EE_converted

TH1F *	h_EE_e25EtrueReference

TH1F *	h_EE_e9EtrueReference

TH1F *	h_EE_eRecoEtrueReference

TH1F *	h_EE_eTrue

std::vector< int >	iMC

std::vector< int >	isConverted

std::vector< float >	mcEnergy

std::vector< float >	mcEta

std::vector< float >	mcPhi

std::vector< float >	mcPt

int	nMCphotons

int	nRECOphotons

edm::EDGetTokenT< EcalRecHitCollection >	reducedBarrelRecHitCollection_

edm::EDGetTokenT< EcalRecHitCollection >	reducedEndcapRecHitCollection_

std::vector< int >	seedXtal

edm::EDGetTokenT< edm::SimTrackContainer >	SimTrackCollection_

edm::EDGetTokenT< edm::SimVertexContainer >	SimVertexCollection_

std::vector< float >	superClusterEnergy

std::vector< float >	superClusterEt

std::vector< float >	superClusterEta

std::vector< float >	superClusterPhi

std::vector< float >	x_vtx

std::vector< float >	y_vtx

std::vector< float >	z_vtx

Additional Inherited Members
Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &)

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

Definition at line 48 of file ContainmentCorrectionAnalyzer.h.

Constructor & Destructor Documentation

ContainmentCorrectionAnalyzer::ContainmentCorrectionAnalyzer ( const edm::ParameterSet & ps )

explicit

Definition at line 11 of file ContainmentCorrectionAnalyzer.cc.

References edm::ParameterSet::getParameter().

                                                                                      {
 
   BarrelSuperClusterCollection_  = consumes<reco::SuperClusterCollection>(ps.getParameter<InputTag>("BarrelSuperClusterCollection"));
   EndcapSuperClusterCollection_  = consumes<reco::SuperClusterCollection>(ps.getParameter<InputTag>("EndcapSuperClusterCollection"));
   reducedBarrelRecHitCollection_ = consumes<EcalRecHitCollection>(ps.getParameter<InputTag>("reducedBarrelRecHitCollection"));
   reducedEndcapRecHitCollection_ = consumes<EcalRecHitCollection>(ps.getParameter<InputTag>("reducedEndcapRecHitCollection"));
   SimTrackCollection_            = consumes<SimTrackContainer>(ps.getParameter<InputTag>("simTrackCollection"));
   SimVertexCollection_           = consumes<SimVertexContainer>(ps.getParameter<InputTag>("simVertexCollection"));
 }

ContainmentCorrectionAnalyzer::~ContainmentCorrectionAnalyzer ( )

Definition at line 21 of file ContainmentCorrectionAnalyzer.cc.

21 { }

Member Function Documentation

void ContainmentCorrectionAnalyzer::analyze	(	const edm::Event &	evt,
		const edm::EventSetup &	es
	)

virtual

Definition at line 40 of file ContainmentCorrectionAnalyzer.cc.

References deltaR(), MillePedeFileConverter_cfg::e, reco::e1, heepElectronID_HEEPV50_CSA14_25ns_cff::e5x5, photonPostprocessing_cfi::eMax, stringResolutionProvider_cfi::eta, HcalObjRepresent::Fill(), edm::EventSetup::get(), edm::Event::getByToken(), mps_fire::i, edm::EventBase::id(), edm::ESHandleBase::isValid(), edm::HandleBase::isValid(), checklumidiff::l, edm::ProductLabels::module, jetCleaner_cfi::photons, pi, reco::CaloCluster::position(), funct::pow(), reco::SuperCluster::preshowerEnergy(), edm::Handle< T >::product(), edm::ESHandle< T >::product(), reco::SuperCluster::rawEnergy(), reco::SuperCluster::seed(), mathSSE::sqrt(), ecaldqm::topology(), and badGlobalMuonTaggersAOD_cff::vtx.

                                                                                     {
   
   LogInfo("ContainmentCorrectionAnalyzer") << "Analyzing event " << evt.id() << "\n";
   
   // taking the needed collections
   std::vector<SimTrack> theSimTracks;
   Handle<SimTrackContainer> SimTk;
   evt.getByToken(SimTrackCollection_,SimTk);
   Labels l;
   labelsForToken(SimTrackCollection_, l);
 
   if (SimTk.isValid() ) theSimTracks.insert(theSimTracks.end(),SimTk->begin(),SimTk->end());  
   else { LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module; }
 
   std::vector<SimVertex> theSimVertexes;  
   Handle<SimVertexContainer> SimVtx;
   evt.getByToken(SimVertexCollection_,SimVtx);
   labelsForToken(SimVertexCollection_, l);
 
 
   if (SimVtx.isValid()) theSimVertexes.insert(theSimVertexes.end(),SimVtx->begin(),SimVtx->end());
   else {LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module; }
   
   const reco::SuperClusterCollection* BarrelSuperClusters = 0;
   Handle<reco::SuperClusterCollection> pHybridBarrelSuperClusters;
   evt.getByToken(BarrelSuperClusterCollection_, pHybridBarrelSuperClusters);
   labelsForToken(BarrelSuperClusterCollection_, l);
 
   if (pHybridBarrelSuperClusters.isValid()) { BarrelSuperClusters = pHybridBarrelSuperClusters.product(); } 
   else {LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module; }
 
   const reco::SuperClusterCollection* EndcapSuperClusters = 0;
   Handle<reco::SuperClusterCollection> pMulti5x5EndcapSuperClusters;
   evt.getByToken(EndcapSuperClusterCollection_, pMulti5x5EndcapSuperClusters);
   labelsForToken(EndcapSuperClusterCollection_, l);
 
   if (pMulti5x5EndcapSuperClusters.isValid()) EndcapSuperClusters = pMulti5x5EndcapSuperClusters.product();
   else {LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module; }
   
   const EcalRecHitCollection *ebRecHits = 0;
   Handle< EcalRecHitCollection > pEBRecHits;
   evt.getByToken( reducedBarrelRecHitCollection_, pEBRecHits );
   labelsForToken(reducedBarrelRecHitCollection_, l);
 
   if (pEBRecHits.isValid()) ebRecHits = pEBRecHits.product();
   else {LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module; }
   
   const EcalRecHitCollection *eeRecHits = 0;
   Handle< EcalRecHitCollection > pEERecHits;
   evt.getByToken( reducedEndcapRecHitCollection_, pEERecHits );
   labelsForToken(reducedEndcapRecHitCollection_, l);
 
   if (pEERecHits.isValid()) eeRecHits = pEERecHits.product();
   else {LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module; }
   
   const CaloTopology *topology = 0;
   ESHandle<CaloTopology> pTopology;
   es.get<CaloTopologyRecord>().get(pTopology);
   if(pTopology.isValid()) topology = pTopology.product();
   
   std::vector<EcalSimPhotonMCTruth> photons=findMcTruth(theSimTracks,theSimVertexes);
   
   nMCphotons = 0;
   nRECOphotons = 0;
   
   int mc_size = photons.size();
   mcEnergy.resize(mc_size);
   mcEta.resize(mc_size);
   mcPhi.resize(mc_size);
   mcPt.resize(mc_size);
   isConverted.resize(mc_size);
   x_vtx.resize(mc_size);
   y_vtx.resize(mc_size);
   z_vtx.resize(mc_size);
   
   superClusterEnergy.resize(mc_size);
   superClusterEta.resize(mc_size);
   superClusterPhi.resize(mc_size);
   superClusterEt.resize(mc_size);
   e1.resize(mc_size);
   e9.resize(mc_size);
   e25.resize(mc_size);
   seedXtal.resize(mc_size);
   iMC.resize(mc_size);
 
 
   // loop over MC truth photons
   for (unsigned int ipho=0;ipho<photons.size();ipho++) {
 
     math::XYZTLorentzVectorD vtx = photons[ipho].primaryVertex();
     double phiTrue = photons[ipho].fourMomentum().phi();
     double vtxPerp = sqrt(vtx.x()*vtx.x() + vtx.y()*vtx.y()); 
     double etaTrue = ecalEta(photons[ipho].fourMomentum().eta(), vtx.z(), vtxPerp);
     double etTrue  = photons[ipho].fourMomentum().e()/cosh(etaTrue);
     nMCphotons++;
     mcEnergy[nMCphotons-1]=photons[ipho].fourMomentum().e();
     mcEta[nMCphotons-1]=etaTrue;
     mcPhi[nMCphotons-1]=phiTrue;
     mcPt[nMCphotons-1]=etTrue;      
     isConverted[nMCphotons-1]=photons[ipho].isAConversion();
     x_vtx[nMCphotons-1]=vtx.x();
     y_vtx[nMCphotons-1]=vtx.y();
     z_vtx[nMCphotons-1]=vtx.z();
 
     // check histos for MC truth
     if(std::fabs(etaTrue) < 1.479) {
       h_EB_eTrue     ->Fill(photons[ipho].fourMomentum().e());
       h_EB_converted ->Fill(photons[ipho].isAConversion());
     }         
     if(std::fabs(etaTrue) >= 1.6) {
       h_EE_eTrue     ->Fill(photons[ipho].fourMomentum().e());
       h_EE_converted ->Fill(photons[ipho].isAConversion());
     }         
 
     // barrel
     if(std::fabs(etaTrue) < 1.479) {
       double etaCurrent; // , etaFound = 0; // UNUSED
       double phiCurrent;
       // double etCurrent,  etFound  = 0; // UNUSED
       const reco::SuperCluster* nearSC = 0;
       
       double closestParticleDistance = 999;      
       for(reco::SuperClusterCollection::const_iterator aClus = BarrelSuperClusters->begin(); 
       aClus != BarrelSuperClusters->end(); aClus++) {       
     etaCurrent = aClus->position().eta();
     phiCurrent = aClus->position().phi();
     // etCurrent  = aClus->energy()/std::cosh(etaCurrent); // UNUSED
     double deltaR = std::sqrt(std::pow(etaCurrent-etaTrue,2)+std::pow(phiCurrent-phiTrue,2));
     if(deltaR < closestParticleDistance) {
       // etFound  = etCurrent; // UNUSED
       // etaFound = etaCurrent; // UNUSED
       closestParticleDistance = deltaR;
       nearSC=&(*aClus);
     }
       }
 
       if(closestParticleDistance < 0.3) {
     // Is a matched particle dumping informations
     nRECOphotons++;
     superClusterEnergy[nRECOphotons-1]=nearSC->rawEnergy();
     superClusterEta[nRECOphotons-1]=nearSC->position().eta();
     superClusterPhi[nRECOphotons-1]=nearSC->position().phi();
     superClusterEt[nRECOphotons-1]=nearSC->rawEnergy()/std::cosh(nearSC->position().eta());
     iMC[nRECOphotons-1]=nMCphotons-1;
 
     reco::CaloClusterPtr theSeed = nearSC->seed();
     seedXtal[nRECOphotons-1] = EcalClusterTools::getMaximum(*theSeed, ebRecHits).first;
     e1[nRECOphotons-1]       = EcalClusterTools::eMax(*theSeed, ebRecHits);
     e9[nRECOphotons-1]       = EcalClusterTools::e3x3(*theSeed, ebRecHits, topology ); 
     e25[nRECOphotons-1]      = EcalClusterTools::e5x5(*theSeed, ebRecHits, topology ); 
       }
     }
     
     // endcap
     if(std::fabs(etaTrue) >= 1.6) {
       double etaCurrent; // , etaFound = 0; // UNUSED
       double phiCurrent;
       // double etCurrent,  etFound  = 0; // UNUSED
       const reco::SuperCluster* nearSC = 0;
       
       double closestParticleDistance = 999; 
       for(reco::SuperClusterCollection::const_iterator aClus = EndcapSuperClusters->begin(); 
       aClus != EndcapSuperClusters->end(); aClus++) {
     etaCurrent = aClus->position().eta();
     phiCurrent = aClus->position().phi();
     // etCurrent  =  aClus->energy()/std::cosh(etaCurrent);
     double deltaR = std::sqrt(std::pow(etaCurrent-etaTrue,2)+std::pow(phiCurrent-phiTrue,2));
     if(deltaR < closestParticleDistance) {
       // etFound  = etCurrent; // UNUSED
       // etaFound = etaCurrent; // UNUSED
       closestParticleDistance = deltaR;
       nearSC=&(*aClus);
     }
       }
       
       if(closestParticleDistance < 0.3) {
     //Is a matched particle dumping informations
     nRECOphotons++;
     float psEnergy = nearSC->preshowerEnergy();
     superClusterEnergy[nRECOphotons-1]=nearSC->rawEnergy() + psEnergy;
     superClusterEta[nRECOphotons-1]=nearSC->position().eta();
     superClusterPhi[nRECOphotons-1]=nearSC->position().phi();
     superClusterEt[nRECOphotons-1]=(nearSC->rawEnergy() + psEnergy)/std::cosh(nearSC->position().eta());
     iMC[nRECOphotons-1]=nMCphotons-1;
 
     reco::CaloClusterPtr theSeed = nearSC->seed();
     seedXtal[nRECOphotons-1] = EcalClusterTools::getMaximum(*theSeed, eeRecHits).first;
     e1[nRECOphotons-1]       = EcalClusterTools::eMax(*theSeed, eeRecHits) + psEnergy;
     e9[nRECOphotons-1]       = EcalClusterTools::e3x3(*theSeed, eeRecHits, topology ) + psEnergy; 
     e25[nRECOphotons-1]      = EcalClusterTools::e5x5(*theSeed, eeRecHits, topology ) + psEnergy; 
       }
     }  
   }  
 
 
   // containment analysis for unconverted photons in the reference region only
   for (int i=0;i<nRECOphotons;i++) {
 
     // barrel
     if (fabs(superClusterEta[i]) < 1.479 ) {
       if (isConverted[iMC[i]] != 1){
     int ietaAbs = (seedXtal[i]>>9)&0x7F;
     int iphi    = seedXtal[i]&0x1FF;
     if (ietaAbs > 5 && ietaAbs < 21 && ((iphi % 20) > 5) && ((iphi % 20) < 16) ) {
       h_EB_eRecoEtrueReference -> Fill(superClusterEnergy[i]/mcEnergy[iMC[i]]);
       h_EB_e9EtrueReference    -> Fill(e9[i]/mcEnergy[iMC[i]]);
       h_EB_e25EtrueReference   -> Fill(e25[i]/mcEnergy[iMC[i]]);
     }
       }
     }    
   
     // endcap
     if (fabs(superClusterEta[i]) > 1.6 ) {
       if (isConverted[iMC[i]] != 1) {
     if ( fabs(superClusterEta[i]) > 1.7 && fabs(superClusterEta[i] < 2.3)  &&
          ((superClusterPhi[i] > -CLHEP::pi/2. + 0.1 && superClusterPhi[i] < CLHEP::pi/2. - 0.1) ||
           (superClusterPhi[i] > CLHEP::pi/2. + 0.1) ||
           (superClusterPhi[i] < -CLHEP::pi/2. - 0.1) 
           ))
       {
         h_EE_eRecoEtrueReference -> Fill(superClusterEnergy[i]/mcEnergy[iMC[i]]);
         h_EE_e9EtrueReference    -> Fill(e9[i]/mcEnergy[iMC[i]]);
         h_EE_e25EtrueReference   -> Fill(e25[i]/mcEnergy[iMC[i]]);
       }
       }
     }
     
   } // loop over reco photons
 }

void ContainmentCorrectionAnalyzer::beginJob ( void )

virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 23 of file ContainmentCorrectionAnalyzer.cc.

References TFileService::make().

                                              {
 
   Service<TFileService> fs;  
 
   // Define reference histograms
   h_EB_eRecoEtrueReference = fs->make<TH1F>("EB_eRecoEtrueReference","EB_eRecoEtrueReference",440,0.,1.1);
   h_EB_e9EtrueReference    = fs->make<TH1F>("EB_e9EtrueReference",   "EB_e9EtrueReference",   440,0.,1.1);
   h_EB_e25EtrueReference   = fs->make<TH1F>("EB_e25EtrueReference",  "EB_e25EtrueReference",  440,0.,1.1);
   h_EE_eRecoEtrueReference = fs->make<TH1F>("EE_eRecoEtrueReference","EE_eRecoEtrueReference",440,0.,1.1);
   h_EE_e9EtrueReference    = fs->make<TH1F>("EE_e9EtrueReference",   "EE_e9EtrueReference",   440,0.,1.1);
   h_EE_e25EtrueReference   = fs->make<TH1F>("EE_e25EtrueReference",  "EE_e25EtrueReference",  440,0.,1.1);
   h_EB_eTrue               = fs->make<TH1F>("EB_eTrue",              "EB_eTrue",              41,40.,60.);
   h_EE_eTrue               = fs->make<TH1F>("EE_eTrue",              "EE_eTrue",              41,40.,60.);
   h_EB_converted           = fs->make<TH1F>("EB_converted",          "EB_converted",          2,0.,2.);
   h_EE_converted           = fs->make<TH1F>("EE_converted",          "EE_converted",          2,0.,2.);
 }

float ContainmentCorrectionAnalyzer::ecalEta	(	float	EtaParticle,
		float	Zvertex,
		float	plane_Radius
	)

private

Definition at line 272 of file ContainmentCorrectionAnalyzer.cc.

References ETA, etaBarrelEndcap, cmsBatch::log, Geom::pi(), R_ECAL, funct::tan(), and Z_Endcap.

                                                                                                   {  
 
   const float R_ECAL           = 136.5;
   const float Z_Endcap         = 328.0;
   const float etaBarrelEndcap  = 1.479;
 
   if(EtaParticle != 0.) {
     float Theta = 0.0  ;
     float ZEcal = (R_ECAL-plane_Radius)*sinh(EtaParticle)+Zvertex;
     
     if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
     if(Theta<0.0) Theta = Theta+Geom::pi() ;
     
     float ETA = - log(tan(0.5*Theta));
     
     if( fabs(ETA) > etaBarrelEndcap ) {
       float Zend = Z_Endcap ;
       if(EtaParticle<0.0 )  Zend = -Zend ;
       float Zlen = Zend - Zvertex ;
       float RR = Zlen/sinh(EtaParticle);
       Theta = atan((RR+plane_Radius)/Zend);
       if(Theta<0.0) Theta = Theta+Geom::pi() ;
       ETA = - log(tan(0.5*Theta));
     }
     
     return ETA;
   }
   else {
     LogWarning("")  << "[ContainmentCorrectionAnalyzer::ecalEta] Warning: Eta equals to zero, not correcting" ;
     return EtaParticle;
   }
 }

void ContainmentCorrectionAnalyzer::endJob ( void )

virtual

Reimplemented from edm::EDAnalyzer.

Definition at line 270 of file ContainmentCorrectionAnalyzer.cc.

270 { }

void ContainmentCorrectionAnalyzer::fillMcTruth	(	std::vector< SimTrack > &	theSimTracks,
		std::vector< SimVertex > &	theSimVertices
	)

private

Definition at line 438 of file ContainmentCorrectionAnalyzer.cc.

                                                                                                                   {
 
   unsigned nVtx = simVertices.size();
   unsigned nTks = simTracks.size();
   if ( nVtx == 0 ) return;
   // create a map associating geant particle id and position in the event SimTrack vector
   for( unsigned it=0; it<nTks; ++it ) {
     geantToIndex_[ simTracks[it].trackId() ] = it;
   }  
 }

std::vector< EcalSimPhotonMCTruth > ContainmentCorrectionAnalyzer::findMcTruth	(	std::vector< SimTrack > &	theSimTracks,
		std::vector< SimVertex > &	theSimVertices
	)

private

store this electron since it's from a converted photon

Definition at line 306 of file ContainmentCorrectionAnalyzer.cc.

References funct::abs(), TrackValidation_cff::association, reco::if(), objects.autophobj::motherId, SimVertex::parentIndex(), CoreSimVertex::position(), and mps_fire::result.

                                                                                                                                                    {
 
   std::vector<EcalSimPhotonMCTruth> result;
 
   geantToIndex_.clear();
   // int   idTrk1_[10]; // UNUSED
   // int   idTrk2_[10]; // UNUSED
   
   // Local variables  
   // const int SINGLE=1; // UNUSED
   // const int DOUBLE=2; // UNUSED
   // const int PYTHIA=3; // UNUSED
   const int ELECTRON_FLAV=1;
   const int PIZERO_FLAV=2;
   const int PHOTON_FLAV=3;
   
   // int ievtype=0; // UNUSED
   int ievflav=0;
   std::vector<SimTrack*> photonTracks;
   std::vector<SimTrack*> pizeroTracks;
   std::vector<const SimTrack *> trkFromConversion;
   SimVertex primVtx;   
   std::vector<int> convInd;
 
   fillMcTruth(theSimTracks,  theSimVertices);
   int iPV=-1;   
   int partType1=0;
   int partType2=0;
   std::vector<SimTrack>::iterator iFirstSimTk = theSimTracks.begin();
   if (  !(*iFirstSimTk).noVertex() ) {
     iPV =  (*iFirstSimTk).vertIndex();
     int vtxId =   (*iFirstSimTk).vertIndex();
     primVtx = theSimVertices[vtxId];  
     partType1 = (*iFirstSimTk).type();
   }
   
   // Look at a second track
   iFirstSimTk++;
   if ( iFirstSimTk!=  theSimTracks.end() ) {    
     if (  (*iFirstSimTk).vertIndex() == iPV) {
       partType2 = (*iFirstSimTk).type();  
     }
   }
   int npv=0;
   int iPho=0;
   for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
     if (  (*iSimTk).noVertex() ) continue;
     // int vertexId = (*iSimTk).vertIndex(); // UNUSED
     // SimVertex vertex = theSimVertices[vertexId]; // UNUSED
     if ( (*iSimTk).vertIndex() == iPV ) {
       npv++;
       if ( (*iSimTk).type() == 22) {
     convInd.push_back(0);
     photonTracks.push_back( &(*iSimTk) );
     // math::XYZTLorentzVectorD momentum = (*iSimTk).momentum(); // UNUSED
       } 
     }
   } 
   
   if(npv > 4) { // ievtype = PYTHIA; // UNUSED
   } else if(npv == 1) {
     if( abs(partType1) == 11 ) { /* ievtype = SINGLE; ==UNUSED== */ ievflav = ELECTRON_FLAV; } 
     else if(partType1 == 111)  { /* ievtype = SINGLE; ==UNUSED== */ ievflav = PIZERO_FLAV; } 
     else if(partType1 == 22)   { /* ievtype = SINGLE; ==UNUSED== */ ievflav = PHOTON_FLAV; }
   } else if(npv == 2) {
     if (  abs(partType1) == 11 && abs(partType2) == 11 ) { /* ievtype = DOUBLE; ==UNUSED== */ ievflav = ELECTRON_FLAV; } 
     else if(partType1 == 111 && partType2 == 111)        { /* ievtype = DOUBLE; ==UNUSED== */ ievflav = PIZERO_FLAV; } 
     else if(partType1 == 22 && partType2 == 22)          { /* ievtype = DOUBLE; ==UNUSED== */ ievflav = PHOTON_FLAV; }
   }      
   
   //  Look into converted photons  
   int isAconversion=0;   
   if(ievflav == PHOTON_FLAV) {
 
     int nConv=0;
     int iConv=0;
     iPho=0;
     for (std::vector<SimTrack*>::iterator iPhoTk = photonTracks.begin(); iPhoTk != photonTracks.end(); ++iPhoTk){
       trkFromConversion.clear();           
       for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk){
     if (  (*iSimTk).noVertex() )                    continue;
     if ( (*iSimTk).vertIndex() == iPV )             continue; 
     if ( abs((*iSimTk).type()) != 11  )             continue;
     int vertexId = (*iSimTk).vertIndex();
     SimVertex vertex = theSimVertices[vertexId];
     int motherId=-1;
     if ( vertex.parentIndex()  ) {
       unsigned  motherGeantId = vertex.parentIndex(); 
       std::map<unsigned, unsigned >::iterator association = geantToIndex_.find( motherGeantId );
       if(association != geantToIndex_.end() )
         motherId = association->second;
       //int motherType = motherId == -1 ? 0 : theSimTracks[motherId].type();
       
       if ( theSimTracks[motherId].trackId() == (*iPhoTk)->trackId() ) {
         trkFromConversion.push_back(&(*iSimTk ) );
       }
     }    
       } // loop over the SimTracks      
        
       if ( trkFromConversion.size() > 0 ) {
     isAconversion=1;
     nConv++;
     convInd[iPho]=nConv;         
     int convVtxId =  trkFromConversion[0]->vertIndex();
     SimVertex convVtx = theSimVertices[convVtxId];
     math::XYZTLorentzVectorD vtxPosition = convVtx.position();
     // math::XYZTLorentzVectorD momentum = (*iPhoTk)->momentum(); // UNUSED
     if ( nConv <= 10) {         
       if ( trkFromConversion.size() > 1) {
         // idTrk1_[iConv]= trkFromConversion[0]->trackId(); // UNUSED
         // idTrk2_[iConv]= trkFromConversion[1]->trackId(); // UNUSED
       } else {
         // idTrk1_[iConv]=trkFromConversion[0]->trackId(); // UNUSED
         // idTrk2_[iConv]=-1; // UNUSED
       }
     }
     iConv++;
        
     result.push_back( EcalSimPhotonMCTruth(isAconversion, (*iPhoTk)->momentum(),   vtxPosition.pt(),  vtxPosition.z() , vtxPosition,   primVtx.position(), trkFromConversion ));
        } else {
          isAconversion=0;
      math::XYZTLorentzVectorD vtxPosition(0.,0.,0.,0.);
          result.push_back( EcalSimPhotonMCTruth(isAconversion, (*iPhoTk)->momentum(),   vtxPosition.pt(),  vtxPosition.z() , vtxPosition,   primVtx.position(), trkFromConversion ));
        }
       iPho++;   
     } // loop over the primary photons
   }   // Event with one or two photons 
 
   return result;
 }