#include <SimpleConvertedPhotonAnalyzer.h>

Inheritance diagram for SimpleConvertedPhotonAnalyzer:

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

void	beginJob () override

void	endJob () override

	SimpleConvertedPhotonAnalyzer (const edm::ParameterSet &)

	~SimpleConvertedPhotonAnalyzer () override

Public Member Functions inherited from edm::one::EDAnalyzer<>
	EDAnalyzer ()=default

SerialTaskQueue *	globalLuminosityBlocksQueue () final

SerialTaskQueue *	globalRunsQueue () final

bool	wantsGlobalLuminosityBlocks () const final

bool	wantsGlobalRuns () const final

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

	EDAnalyzerBase ()

ModuleDescription const &	moduleDescription () const

bool	wantsStreamLuminosityBlocks () const

bool	wantsStreamRuns () const

	~EDAnalyzerBase () override

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

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	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

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

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

virtual	~EDConsumerBase () noexcept(false)

Private Member Functions
float	etaTransformation (float a, float b)

Private Attributes
std::string	convertedPhotonCollection_

std::string	convertedPhotonCollectionProducer_

TFile *	fOutputFile_

std::string	fOutputFileName_

TH2F *	h2_tk_inPtVsR_

TH2F *	h2_tk_nHitsVsR_

TH1F *	h_deltaEta_

TH1F *	h_deltaPhi_

TH1F *	h_ErecoEMC_

TH1F *	h_MCConvE_

TH1F *	h_MCConvEta_

TH1F *	h_MCConvPt_

TH1F *	h_MCphoE_

TH1F *	h_MCphoEta_

TH1F *	h_MCphoPhi_

TH1F *	h_phoE_

TH1F *	h_phoEta_

TH1F *	h_phoPhi_

TH1F *	h_scE_

TH1F *	h_scEta_

TH1F *	h_scPhi_

std::string	HepMCLabel

int	nEvt_

int	nMatched_

int	nMCPho_

std::string	SimHitLabel

std::string	SimTkLabel

std::string	SimVtxLabel

PhotonMCTruthFinder *	thePhotonMCTruthFinder_

Additional Inherited Members
Public Types inherited from edm::one::EDAnalyzerBase
typedef EDAnalyzerBase	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

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

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

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 ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

template<typename ESProduct , Transition Tr = Transition::Event>
auto	esConsumes (eventsetup::EventSetupRecordKey const &, ESInputTag const &tag)

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 22 of file SimpleConvertedPhotonAnalyzer.h.

Constructor & Destructor Documentation

SimpleConvertedPhotonAnalyzer::SimpleConvertedPhotonAnalyzer ( const edm::ParameterSet & pset )

explicit

Definition at line 47 of file SimpleConvertedPhotonAnalyzer.cc.

References convertedPhotonCollection_, convertedPhotonCollectionProducer_, edm::ParameterSet::getParameter(), and AlCaHLTBitMon_QueryRunRegistry::string.

    : fOutputFileName_( pset.getUntrackedParameter<string>("HistOutFile",std::string("TestConversions.root")) ),
      fOutputFile_(nullptr)
 {
 
   convertedPhotonCollectionProducer_ = pset.getParameter<std::string>("phoProducer");
   convertedPhotonCollection_ = pset.getParameter<std::string>("convertedPhotonCollection");
   //
 
 
 }

SimpleConvertedPhotonAnalyzer::~SimpleConvertedPhotonAnalyzer ( )

override

Definition at line 61 of file SimpleConvertedPhotonAnalyzer.cc.

References thePhotonMCTruthFinder_.

                                                               {
 
   delete thePhotonMCTruthFinder_;
 
 }

Member Function Documentation

void SimpleConvertedPhotonAnalyzer::analyze	(	const edm::Event &	e,
		const edm::EventSetup &
	)

override

Loop over recontructed photons

End loop over Reco particles

End loop over MC particles

Definition at line 156 of file SimpleConvertedPhotonAnalyzer.cc.

References convertedPhotonCollection_, convertedPhotonCollectionProducer_, gather_cfg::cout, delta, spr::deltaEta, hiPixelPairStep_cff::deltaPhi, etaTransformation(), HcalObjRepresent::Fill(), PhotonMCTruthFinder::find(), edm::Event::getByLabel(), h2_tk_inPtVsR_, h2_tk_nHitsVsR_, h_deltaEta_, h_deltaPhi_, h_ErecoEMC_, h_MCConvEta_, h_MCphoE_, h_MCphoEta_, h_MCphoPhi_, h_scE_, h_scEta_, h_scPhi_, mps_fire::i, edm::EventBase::id(), nEvt_, nMatched_, nMCPho_, Geom::pi(), funct::pow(), edm::Handle< T >::product(), mathSSE::sqrt(), thePhotonMCTruthFinder_, and Geom::twoPi().

 {
 
   using namespace edm;
   const float etaPhiDistance=0.01;
   // Fiducial region
   //UNUSED const float TRK_BARL =0.9;
   const float BARL = 1.4442; // DAQ TDR p.290
   const float END_LO = 1.566;
   const float END_HI = 2.5;
  // Electron mass
   //UNUSED const Float_t mElec= 0.000511;
 
   nEvt_++;
   LogInfo("ConvertedPhotonAnalyzer") << "ConvertedPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
   //  LogDebug("ConvertedPhotonAnalyzer") << "ConvertedPhotonAnalyzer Analyzing event number: "  << e.id() << " Global Counter " << nEvt_ <<"\n";
   std::cout << "ConvertedPhotonAnalyzer Analyzing event number: "  << e.id() << " Global Counter " << nEvt_ <<"\n";
 
   Handle<reco::ConversionCollection> convertedPhotonHandle;
   e.getByLabel(convertedPhotonCollectionProducer_, convertedPhotonCollection_ , convertedPhotonHandle);
   const reco::ConversionCollection phoCollection = *(convertedPhotonHandle.product());
   std::cout  << "ConvertedPhotonAnalyzer  Converted photon collection size " << phoCollection.size() << "\n";
 
   std::cout  << " ConvertedPhotonAnalyzer Looking for MC truth " << "\n";
 
   //get simtrack info
   std::vector<SimTrack> theSimTracks;
   std::vector<SimVertex> theSimVertices;
 
   edm::Handle<SimTrackContainer> SimTk;
   edm::Handle<SimVertexContainer> SimVtx;
   e.getByLabel("g4SimHits",SimTk);
   e.getByLabel("g4SimHits",SimVtx);
 
   theSimTracks.insert(theSimTracks.end(),SimTk->begin(),SimTk->end());
   theSimVertices.insert(theSimVertices.end(),SimVtx->begin(),SimVtx->end());
 
   std::vector<PhotonMCTruth> mcPhotons=thePhotonMCTruthFinder_->find (theSimTracks,  theSimVertices);
   std::cout << " ConvertedPhotonAnalyzer mcPhotons size " <<  mcPhotons.size() << std::endl;
 
   // Loop over simulated photons
   //UNUSED int iDet=0;
   //UNUSED int iRadius=-1;
   //UNUSED int indPho=0;
 
   for ( std::vector<PhotonMCTruth>::const_iterator mcPho=mcPhotons.begin(); mcPho !=mcPhotons.end(); mcPho++) {
     float mcPhi= (*mcPho).fourMomentum().phi();
     float mcEta= (*mcPho).fourMomentum().pseudoRapidity();
     mcEta = etaTransformation(mcEta, (*mcPho).primaryVertex().z() );
 
     if ( ! (  fabs(mcEta) <= BARL || ( fabs(mcEta) >= END_LO && fabs(mcEta) <=END_HI ) ) ) {
           continue;
     } // all ecal fiducial region
 
     std::cout << " ConvertedPhotonAnalyzer MC Photons before matching  " <<  std::endl;
     std::cout << " ConvertedPhotonAnalyzer Photons isAconversion " << (*mcPho).isAConversion() << " mcMatchingPhoton energy " << (*mcPho).fourMomentum().e()  << " conversion vertex R " << (*mcPho).vertex().perp() << " Z " <<  (*mcPho).vertex().z() <<  " x " << (*mcPho).vertex().x() << " y " <<(*mcPho).vertex().y()  << " z " << (*mcPho).vertex().z() <<  std::endl;
     std::cout << " ConvertedPhotonAnalyzer mcEta " << mcEta<< " mcPhi " << mcPhi <<  std::endl;
 
     h_MCphoE_->Fill(  (*mcPho).fourMomentum().e());
     h_MCphoEta_->Fill( (*mcPho).fourMomentum().eta());
     h_MCphoPhi_->Fill(  (*mcPho).fourMomentum().phi());
 
     // keep only visible conversions
     if (  (*mcPho).isAConversion() == 0 ) continue;
 
     nMCPho_++;
 
     h_MCConvEta_ ->Fill ( fabs( (*mcPho).fourMomentum().pseudoRapidity()) - 0.001);
 
     bool REJECTED;
 
     std::cout   << " ConvertedPhotonAnalyzer  Starting loop over photon candidates " << "\n";
     for( reco::ConversionCollection::const_iterator  iPho = phoCollection.begin(); iPho != phoCollection.end(); iPho++) {
       REJECTED=false;
 
       std::cout  << " ConvertedPhotonAnalyzer Reco SC energy " << (*iPho).caloCluster()[0]->energy() <<  "\n";
 
       float phiClu=(*iPho).caloCluster()[0]->phi();
       float etaClu=(*iPho).caloCluster()[0]->eta();
       float deltaPhi = phiClu-mcPhi;
       float deltaEta = etaClu-mcEta;
 
       if ( deltaPhi > Geom::pi()  ) deltaPhi -= Geom::twoPi();
       if ( deltaPhi < -Geom::pi() ) deltaPhi += Geom::twoPi();
       deltaPhi=std::pow(deltaPhi,2);
       deltaEta=std::pow(deltaEta,2);
       float delta =  deltaPhi+deltaEta ;
       if (  delta >= etaPhiDistance  )  REJECTED=true;
 
       //      if ( ! (  fabs(etaClu) <= BARL || ( fabs(etaClu) >= END_LO && fabs(etaClu) <=END_HI ) ) ) REJECTED=true;
 
       if ( REJECTED ) continue;
       std::cout << " MATCHED " << std::endl;
       nMatched_++;
 
       std::cout << " ConvertedPhotonAnalyzer Matching candidate " << std::endl;
 
       std::cout << " ConvertedPhotonAnalyzer Photons isAconversion " << (*mcPho).isAConversion() << " mcMatchingPhoton energy " <<  (*mcPho).fourMomentum().e()  << " ConvertedPhotonAnalyzer conversion vertex R " <<  (*mcPho).vertex().perp() << " Z " <<  (*mcPho).vertex().z() <<  std::endl;
 
       h_ErecoEMC_->Fill(   (*iPho).caloCluster()[0]->energy()/(*mcPho).fourMomentum().e());
       h_deltaPhi_-> Fill ( (*iPho).caloCluster()[0]->position().phi()- mcPhi);
       h_deltaEta_-> Fill ( (*iPho).caloCluster()[0]->position().eta()- mcEta);
 
       h_scE_->Fill( (*iPho).caloCluster()[0]->energy() );
       h_scEta_->Fill( (*iPho).caloCluster()[0]->position().eta() );
       h_scPhi_->Fill( (*iPho).caloCluster()[0]->position().phi() );
 
       for (unsigned int i=0; i<(*iPho).tracks().size(); i++) {
     std::cout  << " ConvertedPhotonAnalyzer Reco Track charge " <<  (*iPho).tracks()[i]->charge() << "  Num of RecHits " << (*iPho).tracks()[i]->recHitsSize() << " inner momentum " <<  sqrt (  (*iPho).tracks()[i]->innerMomentum().Mag2() )  <<  "\n";
 
         h2_tk_nHitsVsR_ -> Fill (  (*mcPho).vertex().perp(), (*iPho).tracks()[i]->recHitsSize()   );
     h2_tk_inPtVsR_->Fill  (  (*mcPho).vertex().perp(),  sqrt( (*iPho).tracks()[i]->innerMomentum().Mag2() ) );
 
       }
 
     }   
 
   }
 
 }

void SimpleConvertedPhotonAnalyzer::beginJob ( void )

overridevirtual

Reco - MC

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 69 of file SimpleConvertedPhotonAnalyzer.cc.

References fOutputFile_, fOutputFileName_, h2_tk_inPtVsR_, h2_tk_nHitsVsR_, h_deltaEta_, h_deltaPhi_, h_ErecoEMC_, h_MCConvE_, h_MCConvEta_, h_MCConvPt_, h_MCphoE_, h_MCphoEta_, h_MCphoPhi_, h_phoE_, h_phoEta_, h_phoPhi_, h_scE_, h_scEta_, h_scPhi_, nEvt_, nMatched_, nMCPho_, mathSSE::return(), and thePhotonMCTruthFinder_.

 {
 
 
   nEvt_=0;
   nMCPho_=0;
   nMatched_=0;
 
 
   thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
 
   fOutputFile_   = new TFile( fOutputFileName_.c_str(), "RECREATE" ) ;
 
   h_ErecoEMC_  = new TH1F("deltaE","    delta(reco-mc) energy",100,0.,2.);
   h_deltaPhi_ = new TH1F("deltaPhi","  delta(reco-mc) phi",100,-0.1, 0.1);
   h_deltaEta_ = new TH1F("deltaEta","  delta(reco-mc) eta",100,-0.05, 0.05);
 
   h_MCphoE_ = new TH1F("MCphoE","MC photon energy",100,0.,100.);
   h_MCphoPhi_ = new TH1F("MCphoPhi","MC photon phi",40,-3.14, 3.14);
   h_MCphoEta_ = new TH1F("MCphoEta","MC photon eta",40,-3., 3.);
 
 
   h_MCConvE_ = new TH1F("MCConvE","MC converted photon energy",100,0.,100.);
   h_MCConvPt_ = new TH1F("MCConvPt","MC converted photon pt",100,0.,100.);
   h_MCConvEta_ = new TH1F("MCConvEta","MC converted photon eta",50, 0., 2.5);
 
   h_scE_ = new TH1F("scE","Uncorrected converted photons : SC Energy ",100,0., 200.);
   h_scEta_ = new TH1F("scEta","Uncorrected converted photons:  SC Eta ",40,-3., 3.);
   h_scPhi_ = new TH1F("scPhi","Uncorrected converted photons: SC Phi ",40, -3.14, 3.14);
   //
   h_phoE_ = new TH1F("phoE","Uncorrected converted photons :  Energy ",100,0., 200.);
   h_phoEta_ = new TH1F("phoEta","Uncorrected converted photons:   Eta ",40,-3., 3.);
   h_phoPhi_ = new TH1F("phoPhi","Uncorrected converted photons:  Phi ",40,  -3.14, 3.14);
 
   // Recontructed tracks from converted photon candidates
   h2_tk_nHitsVsR_ = new TH2F("tknHitsVsR","Tracks Hits vs R  ", 12,0.,120.,20,0.5, 20.5);
   h2_tk_inPtVsR_ = new TH2F("tkInPtvsR","Tracks inner Pt vs R  ", 12,0.,120.,100,0., 100.);
 
 
   return ;
 }

void SimpleConvertedPhotonAnalyzer::endJob ( void )

overridevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 314 of file SimpleConvertedPhotonAnalyzer.cc.

References gather_cfg::cout, fOutputFile_, nEvt_, and mathSSE::return().

 {
 
    fOutputFile_->Write() ;
    fOutputFile_->Close() ;
 
    edm::LogInfo("ConvertedPhotonAnalyzer") << "Analyzed " << nEvt_  << "\n";
    // std::cout  << "::endJob Analyzed " << nEvt_ << " events " << " with total " << nPho_ << " Photons " << "\n";
    std::cout  << "ConvertedPhotonAnalyzer::endJob Analyzed " << nEvt_ << " events " << "\n";
 
    return ;
 }

float SimpleConvertedPhotonAnalyzer::etaTransformation	(	float	a,
		float	b
	)

private

Definition at line 116 of file SimpleConvertedPhotonAnalyzer.cc.

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

Referenced by analyze().

                                                                                            {
 
 //---Definitions
     const float PI    = 3.1415927;
     //UNUSED const float TWOPI = 2.0*PI;
 
 //---Definitions for ECAL
     const float R_ECAL           = 136.5;
     const float Z_Endcap         = 328.0;
     const float etaBarrelEndcap  = 1.479;
 
 //---ETA correction
 
     float Theta = 0.0  ;
         float ZEcal = R_ECAL*sinh(EtaParticle)+Zvertex;
 
     if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
     if(Theta<0.0) Theta = Theta+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/Zend);
        if(Theta<0.0) Theta = Theta+PI ;
        ETA = - log(tan(0.5*Theta));
       }
 //---Return the result
         return ETA;
 //---end
 }