#include <PreshowerPhiClusterProducer.h>

Inheritance diagram for PreshowerPhiClusterProducer:

Public Types
typedef math::XYZPoint	Point

Public Types inherited from edm::stream::EDProducer<>
typedef CacheContexts< T... >	CacheTypes

typedef CacheTypes::GlobalCache	GlobalCache

typedef AbilityChecker< T... >	HasAbility

typedef CacheTypes::LuminosityBlockCache	LuminosityBlockCache

typedef LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >	LuminosityBlockContext

typedef CacheTypes::LuminosityBlockSummaryCache	LuminosityBlockSummaryCache

typedef CacheTypes::RunCache	RunCache

typedef RunContextT< RunCache, GlobalCache >	RunContext

typedef CacheTypes::RunSummaryCache	RunSummaryCache

Public Member Functions
	PreshowerPhiClusterProducer (const edm::ParameterSet &ps)

void	produce (edm::Event &evt, const edm::EventSetup &es) override

void	set (const edm::EventSetup &es)

	~PreshowerPhiClusterProducer () override

Public Member Functions inherited from edm::stream::EDProducer<>
	EDProducer ()=default

bool	hasAbilityToProduceInLumis () const final

bool	hasAbilityToProduceInRuns () const final

Private Attributes
double	aEta_ [4]

double	alpha0_

double	alpha1_

double	alpha2_

double	alpha3_

std::string	assocSClusterCollection_

double	bEta_ [4]

edm::EDGetTokenT< reco::SuperClusterCollection >	endcapSClusterToken_

edm::ESHandle< ESChannelStatus >	esChannelStatus_

edm::ESHandle< ESEEIntercalibConstants >	esEEInterCalib_

edm::ESHandle< ESGain >	esgain_

edm::ESHandle< ESMIPToGeVConstant >	esMIPToGeV_

edm::ESHandle< ESMissingEnergyCalibration >	esMissingECalib_

float	esPhiClusterDeltaEta_

float	esPhiClusterDeltaPhi_

double	etThresh_

double	gamma0_

double	gamma1_

double	gamma2_

double	gamma3_

double	mip_

int	nEvt_

PreshowerPhiClusterAlgo *	presh_algo

std::string	preshClusterCollectionX_

std::string	preshClusterCollectionY_

edm::EDGetTokenT< EcalRecHitCollection >	preshHitToken_

Detailed Description

Definition at line 23 of file PreshowerPhiClusterProducer.h.

Member Typedef Documentation

typedef math::XYZPoint PreshowerPhiClusterProducer::Point

Definition at line 27 of file PreshowerPhiClusterProducer.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 44 of file PreshowerPhiClusterProducer.cc.

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

                                                                                   {
 
   // use configuration file to setup input/output collection names 
   preshHitToken_   =  
       consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("preshRecHitProducer"));
   
   // Name of a SuperClusterCollection to make associations:
   endcapSClusterToken_   = 
       consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("endcapSClusterProducer"));
   
   // Output collections:
   preshClusterCollectionX_ = ps.getParameter<std::string>("preshClusterCollectionX");
   preshClusterCollectionY_ = ps.getParameter<std::string>("preshClusterCollectionY");
   
   assocSClusterCollection_ = ps.getParameter<std::string>("assocSClusterCollection");
   
   produces< reco::PreshowerClusterCollection >(preshClusterCollectionX_);
   produces< reco::PreshowerClusterCollection >(preshClusterCollectionY_);
   produces< reco::SuperClusterCollection >(assocSClusterCollection_);
   
   float esStripECut          = ps.getParameter<double>("esStripEnergyCut");
   esPhiClusterDeltaEta_      = ps.getParameter<double>("esPhiClusterDeltaEta");
   esPhiClusterDeltaPhi_      = ps.getParameter<double>("esPhiClusterDeltaPhi");
 
   etThresh_                  = ps.getParameter<double>("etThresh");
 
   presh_algo = new PreshowerPhiClusterAlgo(esStripECut);
 }

PreshowerPhiClusterProducer::~PreshowerPhiClusterProducer ( )

override

Definition at line 73 of file PreshowerPhiClusterProducer.cc.

                                                           {
   delete presh_algo;
 }

Member Function Documentation

void PreshowerPhiClusterProducer::produce	(	edm::Event &	evt,
		const edm::EventSetup &	es
	)

override

Definition at line 77 of file PreshowerPhiClusterProducer.cc.

References edm::SortedCollection< T, SORT >::begin(), GetRecoTauVFromDQM_MC_cff::cl2, reco::deltaPhi(), MillePedeFileConverter_cfg::e, reco::e1, reco::e2, DetId::Ecal, EcalPreshower, edm::SortedCollection< T, SORT >::end(), SimCluster::energy(), reco::CaloCluster::energy(), ESCondObjectContainer< T >::find(), edm::EventSetup::get(), edm::Event::getByToken(), ESCondObjectContainer< T >::getMap(), CaloGeometry::getSubdetectorGeometry(), LogTrace, fastPrimaryVertexProducer_cfi::maxDeltaPhi, fftjetdijetfilter_cfi::minDeltaPhi, eostools::move(), point, edm::Handle< T >::product(), edm::PtrVector< T >::push_back(), edm::Event::put(), TrackInfoProducer_cfi::rechits, SimDataFormats::CaloAnalysis::sc, reco::PreshowerCluster::setBCRef(), funct::sin(), edm::SortedCollection< T, SORT >::size(), SimCluster::theta(), X, DOFs::Y, and DOFs::Z.

                                                                                 {
   
   edm::Handle< EcalRecHitCollection >   pRecHits;
   edm::Handle< reco::SuperClusterCollection > pSuperClusters;
   
   // get the ECAL geometry:
   edm::ESHandle<CaloGeometry> geoHandle;
   es.get<CaloGeometryRecord>().get(geoHandle);
   
   // retrieve ES-EE intercalibration constants and channel status
   set(es);
   const ESChannelStatus *channelStatus = esChannelStatus_.product();
   
   const CaloSubdetectorGeometry *geometry_p = (geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower));
   
   // create unique_ptr to a PreshowerClusterCollection
   auto clusters_p1 = std::make_unique<reco::PreshowerClusterCollection>();
   auto clusters_p2 = std::make_unique<reco::PreshowerClusterCollection>();
   // create new collection of corrected super clusters
   auto superclusters_p = std::make_unique<reco::SuperClusterCollection>();
   
   CaloSubdetectorTopology * topology_p=nullptr;
   if (geometry_p)
     topology_p  = new EcalPreshowerTopology(geoHandle);
   
   // fetch the product (pSuperClusters)
   evt.getByToken(endcapSClusterToken_, pSuperClusters);   
   const reco::SuperClusterCollection* SClusts = pSuperClusters.product();
   
   // fetch the product (RecHits)
   evt.getByToken( preshHitToken_, pRecHits);
   // pointer to the object in the product
   const EcalRecHitCollection* rechits = pRecHits.product(); // EcalRecHitCollection hit_collection = *rhcHandle;
   
   LogTrace("EcalClusters") << "PreshowerPhiClusterProducerInfo: ### Total # of preshower RecHits: "<< rechits->size();
   
   // make the map of rechits:
   std::map<DetId, EcalRecHit> rechits_map;
   EcalRecHitCollection::const_iterator it;
   for (it = rechits->begin(); it != rechits->end(); it++) {
     // remove bad ES rechits
     if (it->recoFlag()==1 || it->recoFlag()==14 || (it->recoFlag()<=10 && it->recoFlag()>=5)) continue;
     //Make the map of DetID, EcalRecHit pairs
     rechits_map.insert(std::make_pair(it->id(), *it));   
   }
   // The set of used DetID's for a given event:
   std::set<DetId> used_strips;
   used_strips.clear();
   LogTrace("EcalClusters") << "PreshowerPhiClusterProducerInfo: ### rechits_map of size " << rechits_map.size() <<" was created!";
   
   reco::PreshowerClusterCollection clusters1, clusters2;   // output collection of corrected PCs
   reco::SuperClusterCollection new_SC; // output collection of corrected SCs
   
   //make cycle over super clusters
   reco::SuperClusterCollection::const_iterator it_super;
   int isc  = 0;
   for (it_super=SClusts->begin();  it_super!=SClusts->end(); ++it_super) {     
 
     float e1     = 0;
     float e2     = 0;
     float deltaE = 0;
 
     reco::CaloClusterPtrVector new_BC; 
     ++isc;
     LogTrace("EcalClusters")<< " superE = " << it_super->energy() << " superETA = " << it_super->eta() << " superPHI = " << it_super->phi() ;
     
     //cout<<"=== new SC ==="<<endl;
     //cout<<"superE = "<<it_super->energy()<<" superETA = "<<it_super->eta()<<" superPHI = "<<it_super->phi()<<endl;
 
     int nBC = 0;
     int condP1 = 1; // 0: dead channel; 1: active channel
     int condP2 = 1;
     float maxDeltaPhi = 0;
     float minDeltaPhi = 0;
     float refPhi = 0; 
 
     reco::CaloCluster_iterator bc_iter = it_super->clustersBegin();
     for ( ; bc_iter !=it_super->clustersEnd(); ++bc_iter) {
       if (nBC == 0) {
     refPhi = (*bc_iter)->phi();
       } else {
     if (reco::deltaPhi((*bc_iter)->phi(), refPhi) > 0 && reco::deltaPhi((*bc_iter)->phi(), refPhi) > maxDeltaPhi)
       maxDeltaPhi = reco::deltaPhi((*bc_iter)->phi(), refPhi);
     if (reco::deltaPhi((*bc_iter)->phi(), refPhi) < 0 && reco::deltaPhi((*bc_iter)->phi(), refPhi) < minDeltaPhi)
       minDeltaPhi = reco::deltaPhi((*bc_iter)->phi(), refPhi);
     //cout<<"delta phi : "<<reco::deltaPhi((*bc_iter)->phi(), refPhi)<<endl;
       }
       //cout<<"BC : "<<nBC<<" "<<(*bc_iter)->energy()<<" "<<(*bc_iter)->eta()<<" "<<(*bc_iter)->phi()<<endl; 
       nBC++;
     }
     maxDeltaPhi += esPhiClusterDeltaPhi_;
     minDeltaPhi -= esPhiClusterDeltaPhi_;
 
     nBC = 0;
     for (bc_iter = it_super->clustersBegin() ; bc_iter !=it_super->clustersEnd(); ++bc_iter) {  
       if (geometry_p) {
     
     // Get strip position at intersection point of the line EE - Vertex:
     double X = (*bc_iter)->x();
     double Y = (*bc_iter)->y();
     double Z = (*bc_iter)->z();        
     const GlobalPoint point(X,Y,Z);    
     
     DetId tmp1 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 1);
     DetId tmp2 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 2);
     ESDetId strip1 = (tmp1 == DetId(0)) ? ESDetId(0) : ESDetId(tmp1);
     ESDetId strip2 = (tmp2 == DetId(0)) ? ESDetId(0) : ESDetId(tmp2);     
 
     if (nBC == 0) {
       if (strip1 != ESDetId(0) && strip2 != ESDetId(0)) {
         ESChannelStatusMap::const_iterator status_p1 = channelStatus->getMap().find(strip1);
         ESChannelStatusMap::const_iterator status_p2 = channelStatus->getMap().find(strip2);
         if (status_p1->getStatusCode() == 1) condP1 = 0;
         if (status_p2->getStatusCode() == 1) condP2 = 0;
       } else if (strip1 == ESDetId(0)) {
         condP1 = 0;
       } else if (strip2 == ESDetId(0)) {
         condP2 = 0;
       }
       
       //cout<<"starting cluster : "<<maxDeltaPhi<<" "<<minDeltaPhi<<" "<<esPhiClusterDeltaEta_<<endl;
       //cout<<"do plane 1 === "<<strip1.zside()<<" "<<strip1.plane()<<" "<<strip1.six()<<" "<<strip1.siy()<<" "<<strip1.strip()<<endl;
       // Get a vector of ES clusters (found by the PreshSeeded algorithm) associated with a given EE basic cluster.           
       reco::PreshowerCluster cl1 = presh_algo->makeOneCluster(strip1,&used_strips,&rechits_map,geometry_p,topology_p,esPhiClusterDeltaEta_,minDeltaPhi,maxDeltaPhi);   
       cl1.setBCRef(*bc_iter);
       clusters1.push_back(cl1);
       e1 += cl1.energy();       
       
       //cout<<"do plane 2 === "<<strip2.zside()<<" "<<strip2.plane()<<" "<<strip2.six()<<" "<<strip2.siy()<<" "<<strip2.strip()<<endl;
       reco::PreshowerCluster cl2 = presh_algo->makeOneCluster(strip2,&used_strips,&rechits_map,geometry_p,topology_p,esPhiClusterDeltaEta_,minDeltaPhi,maxDeltaPhi); 
       cl2.setBCRef(*bc_iter);
       clusters2.push_back(cl2);
       e2 += cl2.energy();
     }
       }
 
       new_BC.push_back(*bc_iter);
       nBC++;
     }  // end of cycle over BCs
 
     LogTrace("EcalClusters") << " For SC #" << isc-1 << ", containing " 
                  << it_super->clustersSize()                  
                  << " basic clusters, PreshowerPhiClusterAlgo made " 
                  << clusters1.size() << " in X plane and " 
                  << clusters2.size()       
                  << " in Y plane " << " preshower clusters " ;
     
     // update energy of the SuperCluster    
     if(e1+e2 > 1.0e-10) {
 
       e1 = e1 / mip_; // GeV to #MIPs
       e2 = e2 / mip_;
       
       if (condP1 == 1 && condP2 == 1) {
     deltaE = gamma0_*(e1 + alpha0_*e2);       
       } else if (condP1 == 1 && condP2 == 0) {
     deltaE = gamma1_*(e1 + alpha1_*e2);       
       } else if (condP1 == 0 && condP2 == 1) {
     deltaE = gamma2_*(e1 + alpha2_*e2);       
       } else if (condP1 == 0 && condP2 == 0) {
     deltaE = gamma3_*(e1 + alpha3_*e2);       
       }
     }
     
     //corrected Energy
     float E = it_super->energy() + deltaE;
     
     LogTrace("EcalClusters") << " Creating corrected SC ";
     
     reco::SuperCluster sc(E, it_super->position(), it_super->seed(), new_BC, deltaE);
     sc.setPreshowerEnergyPlane1(e1*mip_);
     sc.setPreshowerEnergyPlane2(e2*mip_);
     if (condP1 == 1 && condP2 == 1) sc.setPreshowerPlanesStatus(0);
     else if (condP1 == 1 && condP2 == 0) sc.setPreshowerPlanesStatus(1);
     else if (condP1 == 0 && condP2 == 1) sc.setPreshowerPlanesStatus(2);
     else if (condP1 == 0 && condP2 == 0) sc.setPreshowerPlanesStatus(3);
 
     if (etThresh_>0){ // calling postion().theta can be expensive
         if (sc.energy()*sin(sc.position().theta())>etThresh_ ) 
             new_SC.push_back(sc);
     } else {
         new_SC.push_back(sc);
     }
     
   } // end of cycle over SCs
   
   // copy the preshower clusters into collections and put in the Event:
   clusters_p1->assign(clusters1.begin(), clusters1.end());
   clusters_p2->assign(clusters2.begin(), clusters2.end());
   // put collection of preshower clusters to the event
   evt.put(std::move(clusters_p1), preshClusterCollectionX_ );
   evt.put(std::move(clusters_p2), preshClusterCollectionY_ );
   LogTrace("EcalClusters") << "Preshower clusters added to the event" ;
   
   // put collection of corrected super clusters to the event
   superclusters_p->assign(new_SC.begin(), new_SC.end());
   evt.put(std::move(superclusters_p), assocSClusterCollection_);
   LogTrace("EcalClusters") << "Corrected SClusters added to the event" ;
   
   if (topology_p) delete topology_p;
 }

void PreshowerPhiClusterProducer::set ( const edm::EventSetup & es )

Definition at line 279 of file PreshowerPhiClusterProducer.cc.

                                                              {
   
   es.get<ESGainRcd>().get(esgain_);
   const ESGain *gain = esgain_.product();
   
   double ESGain = gain->getESGain();
 
   es.get<ESMIPToGeVConstantRcd>().get(esMIPToGeV_);
   const ESMIPToGeVConstant *mipToGeV = esMIPToGeV_.product();
 
   mip_ = (ESGain == 1) ? mipToGeV->getESValueLow() : mipToGeV->getESValueHigh(); 
 
   es.get<ESChannelStatusRcd>().get(esChannelStatus_);
 
   es.get<ESEEIntercalibConstantsRcd>().get(esEEInterCalib_);
   const ESEEIntercalibConstants *esEEInterCalib = esEEInterCalib_.product();
 
   // both planes work
   gamma0_ = (ESGain == 1) ? 0.02 : esEEInterCalib->getGammaHigh0();
   alpha0_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow0() : esEEInterCalib->getAlphaHigh0();
 
   // only first plane works 
   gamma1_ = (ESGain == 1) ? (0.02 * esEEInterCalib->getGammaLow1()) : esEEInterCalib->getGammaHigh1();
   alpha1_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow1() : esEEInterCalib->getAlphaHigh1();
 
   // only second plane works
   gamma2_ = (ESGain == 1) ? (0.02 * esEEInterCalib->getGammaLow2()) : esEEInterCalib->getGammaHigh2();
   alpha2_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow2() : esEEInterCalib->getAlphaHigh2();
 
   // both planes do not work
   gamma3_ = (ESGain == 1) ? 0.02 : esEEInterCalib->getGammaHigh3();
   alpha3_ = (ESGain == 1) ? esEEInterCalib->getAlphaLow3() : esEEInterCalib->getAlphaHigh3();
 
   es.get<ESMissingEnergyCalibrationRcd>().get(esMissingECalib_);
   const ESMissingEnergyCalibration * esMissingECalib = esMissingECalib_.product();
 
   // |eta| < 1.9
   aEta_[0] = esMissingECalib->getConstAEta0();
   bEta_[0] = esMissingECalib->getConstBEta0();
 
   // 1.9 < |eta| < 2.1
   aEta_[1] = esMissingECalib->getConstAEta1();
   bEta_[1] = esMissingECalib->getConstBEta1();
 
   // 2.1 < |eta| < 2.3
   aEta_[2] = esMissingECalib->getConstAEta2();
   bEta_[2] = esMissingECalib->getConstBEta2();
 
   // 2.3 < |eta| < 2.5
   aEta_[3] = esMissingECalib->getConstAEta3();
   bEta_[3] = esMissingECalib->getConstBEta3();
 
 }