ContainmentCorrectionAnalyzer Class Reference

#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 ()
ProductHolderIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const
void	itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
void	itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const
void	labelsForToken (EDGetToken iToken, Labels &oLabels) const
void	modulesDependentUpon (std::string const &iProcessName, std::string const &iModuleLabel, bool iPrint, std::vector< char const * > &oModuleLabels) const
void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const
bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const
bool	registeredToConsumeMany (TypeID const &, BranchType) const
void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)
virtual	~EDConsumerBase ()

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.

{ }

Member Function Documentation

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

virtual

Implements edm::EDAnalyzer.

Definition at line 40 of file ContainmentCorrectionAnalyzer.cc.

References deltaR(), alignCSCRings::e, reco::e1, eta, HcalObjRepresent::Fill(), edm::EventSetup::get(), edm::Event::getByToken(), i, edm::EventBase::id(), edm::ESHandleBase::isValid(), edm::HandleBase::isValid(), cmsLHEtoEOSManager::l, edm::ProductLabels::module, interactiveExample::photons, pi, reco::CaloCluster::position(), funct::pow(), reco::SuperCluster::preshowerEnergy(), edm::Handle< T >::product(), edm::ESHandle< class >::product(), reco::SuperCluster::rawEnergy(), reco::SuperCluster::seed(), mathSSE::sqrt(), and ecaldqm::topology().

                                                                                    {
  
  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, dqm-mbProfile::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.

{ }

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(), reco::if(), SimVertex::parentIndex(), CoreSimVertex::position(), and query::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;
}

Member Data Documentation

edm::EDGetTokenT<reco::SuperClusterCollection> ContainmentCorrectionAnalyzer::BarrelSuperClusterCollection_

private

Definition at line 65 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::e1

private

Definition at line 76 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::e25

private

Definition at line 76 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::e9

private

Definition at line 76 of file ContainmentCorrectionAnalyzer.h.

edm::EDGetTokenT<reco::SuperClusterCollection> ContainmentCorrectionAnalyzer::EndcapSuperClusterCollection_

private

Definition at line 66 of file ContainmentCorrectionAnalyzer.h.

std::map<unsigned, unsigned> ContainmentCorrectionAnalyzer::geantToIndex_

private

Definition at line 80 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EB_converted

private

Definition at line 90 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EB_e25EtrueReference

private

Definition at line 84 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EB_e9EtrueReference

private

Definition at line 83 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EB_eRecoEtrueReference

private

Definition at line 82 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EB_eTrue

private

Definition at line 88 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EE_converted

private

Definition at line 91 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EE_e25EtrueReference

private

Definition at line 87 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EE_e9EtrueReference

private

Definition at line 86 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EE_eRecoEtrueReference

private

Definition at line 85 of file ContainmentCorrectionAnalyzer.h.

TH1F* ContainmentCorrectionAnalyzer::h_EE_eTrue

private

Definition at line 89 of file ContainmentCorrectionAnalyzer.h.

std::vector<int> ContainmentCorrectionAnalyzer::iMC

private

Definition at line 78 of file ContainmentCorrectionAnalyzer.h.

std::vector<int> ContainmentCorrectionAnalyzer::isConverted

private

Definition at line 72 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::mcEnergy

private

Definition at line 71 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::mcEta

private

Definition at line 71 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::mcPhi

private

Definition at line 71 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::mcPt

private

Definition at line 71 of file ContainmentCorrectionAnalyzer.h.

int ContainmentCorrectionAnalyzer::nMCphotons

private

Definition at line 70 of file ContainmentCorrectionAnalyzer.h.

int ContainmentCorrectionAnalyzer::nRECOphotons

private

Definition at line 74 of file ContainmentCorrectionAnalyzer.h.

edm::EDGetTokenT<EcalRecHitCollection> ContainmentCorrectionAnalyzer::reducedBarrelRecHitCollection_

private

Definition at line 67 of file ContainmentCorrectionAnalyzer.h.

edm::EDGetTokenT<EcalRecHitCollection> ContainmentCorrectionAnalyzer::reducedEndcapRecHitCollection_

private

Definition at line 68 of file ContainmentCorrectionAnalyzer.h.

std::vector<int> ContainmentCorrectionAnalyzer::seedXtal

private

Definition at line 77 of file ContainmentCorrectionAnalyzer.h.

edm::EDGetTokenT<edm::SimTrackContainer> ContainmentCorrectionAnalyzer::SimTrackCollection_

private

Definition at line 63 of file ContainmentCorrectionAnalyzer.h.

edm::EDGetTokenT<edm::SimVertexContainer> ContainmentCorrectionAnalyzer::SimVertexCollection_

private

Definition at line 64 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::superClusterEnergy

private

Definition at line 75 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::superClusterEt

private

Definition at line 75 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::superClusterEta

private

Definition at line 75 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::superClusterPhi

private

Definition at line 75 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::x_vtx

private

Definition at line 73 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::y_vtx

private

Definition at line 73 of file ContainmentCorrectionAnalyzer.h.

std::vector<float> ContainmentCorrectionAnalyzer::z_vtx

private

Definition at line 73 of file ContainmentCorrectionAnalyzer.h.