#include <Pi0FixedMassWindowCalibration.h>

Inheritance diagram for Pi0FixedMassWindowCalibration:

Public Member Functions
virtual void	beginOfJob ()
	Called at beginning of job. More...

virtual Status	duringLoop (const edm::Event &, const edm::EventSetup &)
	Called at each event. More...

virtual void	endOfJob ()
	Called at end of job. More...

virtual Status	endOfLoop (const edm::EventSetup &, unsigned int iLoop)
	Called at end of loop. More...

	Pi0FixedMassWindowCalibration (const edm::ParameterSet &iConfig)
	Constructor. More...

virtual void	produce (edm::Event &, const edm::EventSetup &)
	Dummy implementation (job done in duringLoop) More...

virtual void	startingNewLoop (unsigned int iLoop)
	Called at beginning of loop. More...

	~Pi0FixedMassWindowCalibration ()
	Destructor. More...

Public Member Functions inherited from edm::ESProducerLooper
	ESProducerLooper ()

virtual std::set < eventsetup::EventSetupRecordKey >	modifyingRecords () const

Public Member Functions inherited from edm::ESProducer
	ESProducer ()

virtual	~ESProducer ()

Public Member Functions inherited from edm::ESProxyFactoryProducer
	ESProxyFactoryProducer ()

virtual void	newInterval (const eventsetup::EventSetupRecordKey &iRecordType, const ValidityInterval &iInterval)
	overrides DataProxyProvider method More...

virtual	~ESProxyFactoryProducer ()

Public Member Functions inherited from edm::eventsetup::DataProxyProvider
	DataProxyProvider ()

const ComponentDescription &	description () const

bool	isUsingRecord (const EventSetupRecordKey &) const

const KeyedProxies &	keyedProxies (const EventSetupRecordKey &iRecordKey) const

void	resetProxies (const EventSetupRecordKey &iRecordType)

void	resetProxiesIfTransient (const EventSetupRecordKey &iRecordType)

void	setAppendToDataLabel (const edm::ParameterSet &)

void	setDescription (const ComponentDescription &iDescription)

std::set< EventSetupRecordKey >	usingRecords () const

virtual	~DataProxyProvider ()

Public Member Functions inherited from edm::EventSetupRecordIntervalFinder
const eventsetup::ComponentDescription &	descriptionForFinder () const

	EventSetupRecordIntervalFinder ()

std::set < eventsetup::EventSetupRecordKey >	findingForRecords () const

const ValidityInterval &	findIntervalFor (const eventsetup::EventSetupRecordKey &, const IOVSyncValue &)

void	setDescriptionForFinder (const eventsetup::ComponentDescription &iDescription)

virtual	~EventSetupRecordIntervalFinder ()

Public Member Functions inherited from edm::EDLooper
	EDLooper ()

	EDLooper (EDLooper const &)=delete

EDLooper &	operator= (EDLooper const &)=delete

virtual	~EDLooper ()

Public Member Functions inherited from edm::EDLooperBase
virtual void	attachTo (ActivityRegistry &)
	Override this method if you need to monitor the state of the processing. More...

virtual void	beginOfJob (EventSetup const &)

void	copyInfo (ScheduleInfo const &)

void	doBeginLuminosityBlock (LuminosityBlockPrincipal &, EventSetup const &)

void	doBeginRun (RunPrincipal &, EventSetup const &)

Status	doDuringLoop (EventPrincipal &eventPrincipal, EventSetup const &es, ProcessingController &)

void	doEndLuminosityBlock (LuminosityBlockPrincipal &, EventSetup const &)

Status	doEndOfLoop (EventSetup const &es)

void	doEndRun (RunPrincipal &, EventSetup const &)

void	doStartingNewLoop ()

	EDLooperBase ()

	EDLooperBase (EDLooperBase const &)=delete

EDLooperBase &	operator= (EDLooperBase const &)=delete

void	prepareForNextLoop (eventsetup::EventSetupProvider *esp)

void	setActionTable (ActionTable const *actionTable)

void	setModuleChanger (ModuleChanger const *)

virtual	~EDLooperBase ()

Private Attributes
std::vector< DetId >	barrelCells

std::string	barrelClusterCollection_

std::string	barrelClusterShapeAssociation_

std::string	barrelHitCollection_

std::string	barrelHitProducer_

std::string	barrelHits_

std::string	clustershapecollectionEB_

std::string	ecalHitsProducer_

const EcalRecHitCollection *	ecalRecHitBarrelCollection

bool	isfirstcall_

IslandClusterAlgo *	island_p

double	mwxtals [85][360][2]

int	nevent

double	newCalibs_barl [85][360][2]

double	oldCalibs_barl [85][360][2]

PositionCalc	posCalculator_

const EcalRecHitCollection *	recalibEcalRecHitCollection

std::map< DetId, EcalRecHit > *	recHitsEB_map

double	selePi0DetaBelt_

double	selePi0DRBelt_

double	selePi0EtBeltIsoRatioMax_

double	selePi0MinvMeanFixed_

double	selePi0MinvSigmaFixed_

double	selePi0PtGammaOneMin_

double	selePi0PtGammaTwoMin_

double	selePi0PtPi0Min_

double	selePi0S4S9GammaOneMin_

double	selePi0S4S9GammaTwoMin_

double	selePi0S9S25GammaOneMin_

double	selePi0S9S25GammaTwoMin_

ClusterShapeAlgo	shapeAlgo_

TFile *	theFile

unsigned int	theMaxLoops

edm::ParameterSet	theParameterSet

IslandClusterAlgo::VerbosityLevel	verbosity

double	wxtals [85][360][2]

Additional Inherited Members
Public Types inherited from edm::eventsetup::DataProxyProvider
typedef std::vector< std::pair < DataKey, boost::shared_ptr < DataProxy > > >	KeyedProxies

typedef std::vector < EventSetupRecordKey >	Keys

typedef std::map < EventSetupRecordKey, KeyedProxies >	RecordProxies

Public Types inherited from edm::EDLooperBase
enum	Status { kContinue, kStop }

Static Public Member Functions inherited from edm::eventsetup::DataProxyProvider
static void	prevalidate (ConfigurationDescriptions &)

Protected Member Functions inherited from edm::ESProducerLooper
virtual void	registerFactoryWithKey (const eventsetup::EventSetupRecordKey &iRecord, std::auto_ptr< eventsetup::ProxyFactoryBase > &iFactory, const std::string &iLabel=std::string())

void	setIntervalFor (const eventsetup::EventSetupRecordKey &iKey, const IOVSyncValue &iTime, ValidityInterval &oInterval)

Protected Member Functions inherited from edm::ESProducer
template<typename T >
void	setWhatProduced (T *iThis, const es::Label &iLabel=es::Label())

template<typename T >
void	setWhatProduced (T iThis, const char iLabel)

template<typename T >
void	setWhatProduced (T *iThis, const std::string &iLabel)

template<typename T , typename TDecorator >
void	setWhatProduced (T *iThis, const TDecorator &iDec, const es::Label &iLabel=es::Label())

template<typename T , typename TReturn , typename TRecord >
void	setWhatProduced (T iThis, TReturn(T::iMethod)(const TRecord &), const es::Label &iLabel=es::Label())

template<typename T , typename TReturn , typename TRecord , typename TArg >
void	setWhatProduced (T iThis, TReturn(T::iMethod)(const TRecord &), const TArg &iDec, const es::Label &iLabel=es::Label())

Protected Member Functions inherited from edm::ESProxyFactoryProducer
template<class TFactory >
void	registerFactory (std::auto_ptr< TFactory > iFactory, const std::string &iLabel=std::string())

virtual void	registerProxies (const eventsetup::EventSetupRecordKey &iRecord, KeyedProxies &aProxyList)
	override DataProxyProvider method More...

Protected Member Functions inherited from edm::eventsetup::DataProxyProvider
void	eraseAll (const EventSetupRecordKey &iRecordKey)
	deletes all the Proxies in aStream More...

void	invalidateProxies (const EventSetupRecordKey &iRecordKey)

template<class T >
void	usingRecord ()

void	usingRecordWithKey (const EventSetupRecordKey &)

Protected Member Functions inherited from edm::EventSetupRecordIntervalFinder
template<class T >
void	findingRecord ()

void	findingRecordWithKey (const eventsetup::EventSetupRecordKey &)

Protected Member Functions inherited from edm::EDLooperBase
ModuleChanger const *	moduleChanger () const
	This only returns a non-zero value during the call to endOfLoop. More...

ScheduleInfo const *	scheduleInfo () const
	This returns a non-zero value after the constructor has been called. More...

Detailed Description

Definition at line 58 of file Pi0FixedMassWindowCalibration.h.

Constructor & Destructor Documentation

Pi0FixedMassWindowCalibration::Pi0FixedMassWindowCalibration ( const edm::ParameterSet & iConfig )

Constructor.

Definition at line 34 of file Pi0FixedMassWindowCalibration.cc.

References barrelClusterCollection_, barrelClusterShapeAssociation_, clustershapecollectionEB_, gather_cfg::cout, edm::ParameterSet::getParameter(), island_p, IslandClusterAlgo::pDEBUG, IslandClusterAlgo::pERROR, IslandClusterAlgo::pINFO, posCalculator_, IslandClusterAlgo::pWARNING, selePi0DetaBelt_, selePi0DRBelt_, selePi0EtBeltIsoRatioMax_, selePi0MinvMeanFixed_, selePi0MinvSigmaFixed_, selePi0PtGammaOneMin_, selePi0PtGammaTwoMin_, selePi0PtPi0Min_, selePi0S4S9GammaOneMin_, selePi0S4S9GammaTwoMin_, selePi0S9S25GammaOneMin_, selePi0S9S25GammaTwoMin_, shapeAlgo_, AlCaHLTBitMon_QueryRunRegistry::string, theParameterSet, and verbosity.

                                                                                            :
   theMaxLoops( iConfig.getUntrackedParameter<unsigned int>("maxLoops",0) ),
   ecalHitsProducer_( iConfig.getParameter< std::string > ("ecalRecHitsProducer") ),
   barrelHits_( iConfig.getParameter< std::string > ("barrelHitCollection") )
 {
 
   std::cout << "[Pi0FixedMassWindowCalibration] Constructor "<<std::endl;
   // The verbosity level
   std::string verbosityString = iConfig.getParameter<std::string>("VerbosityLevel");
   if      (verbosityString == "DEBUG")   verbosity = IslandClusterAlgo::pDEBUG;
   else if (verbosityString == "WARNING") verbosity = IslandClusterAlgo::pWARNING;
   else if (verbosityString == "INFO")    verbosity = IslandClusterAlgo::pINFO;
   else                                   verbosity = IslandClusterAlgo::pERROR;
 
   // The names of the produced cluster collections
   barrelClusterCollection_  = iConfig.getParameter<std::string>("barrelClusterCollection");
 
   // Island algorithm parameters
   double barrelSeedThreshold = iConfig.getParameter<double>("IslandBarrelSeedThr");
   double endcapSeedThreshold = iConfig.getParameter<double>("IslandEndcapSeedThr");
 
 
   // Selection algorithm parameters
   selePi0PtGammaOneMin_ = iConfig.getParameter<double>("selePi0PtGammaOneMin");
   selePi0PtGammaTwoMin_ = iConfig.getParameter<double>("selePi0PtGammaTwoMin");
 
   selePi0DRBelt_ = iConfig.getParameter<double>("selePi0DRBelt");
   selePi0DetaBelt_ = iConfig.getParameter<double>("selePi0DetaBelt");
 
   selePi0PtPi0Min_ = iConfig.getParameter<double>("selePi0PtPi0Min");
 
   selePi0S4S9GammaOneMin_ = iConfig.getParameter<double>("selePi0S4S9GammaOneMin");
   selePi0S4S9GammaTwoMin_ = iConfig.getParameter<double>("selePi0S4S9GammaTwoMin");
   selePi0S9S25GammaOneMin_ = iConfig.getParameter<double>("selePi0S9S25GammaOneMin");
   selePi0S9S25GammaTwoMin_ = iConfig.getParameter<double>("selePi0S9S25GammaTwoMin");
  
   selePi0EtBeltIsoRatioMax_ = iConfig.getParameter<double>("selePi0EtBeltIsoRatioMax");
 
   selePi0MinvMeanFixed_ = iConfig.getParameter<double>("selePi0MinvMeanFixed");
   selePi0MinvSigmaFixed_ = iConfig.getParameter<double>("selePi0MinvSigmaFixed");
 
 
 
 
   // Parameters for the position calculation:
   edm::ParameterSet posCalcParameters = 
     iConfig.getParameter<edm::ParameterSet>("posCalcParameters");
   posCalculator_ = PositionCalc(posCalcParameters); 
   shapeAlgo_ = ClusterShapeAlgo(posCalcParameters);
   clustershapecollectionEB_ = iConfig.getParameter<std::string>("clustershapecollectionEB");
 
   //AssociationMap
   barrelClusterShapeAssociation_ = iConfig.getParameter<std::string>("barrelShapeAssociation");
 
   island_p = new IslandClusterAlgo(barrelSeedThreshold, endcapSeedThreshold, posCalculator_,verbosity);
 
   theParameterSet=iConfig;
 
 
 }

Pi0FixedMassWindowCalibration::~Pi0FixedMassWindowCalibration ( )

Destructor.

Definition at line 99 of file Pi0FixedMassWindowCalibration.cc.

 {
 
 }

Member Function Documentation

void Pi0FixedMassWindowCalibration::beginOfJob ( )

virtual

Called at beginning of job.

Reimplemented from edm::EDLooperBase.

Definition at line 107 of file Pi0FixedMassWindowCalibration.cc.

References isfirstcall_.

 {
 
   //std::cout << "[Pi0FixedMassWindowCalibration] beginOfJob "<<std::endl;
 
   isfirstcall_=true;
   
 }

edm::EDLooper::Status Pi0FixedMassWindowCalibration::duringLoop	(	const edm::Event &	event,
		const edm::EventSetup &	setup
	)

virtual

Called at each event.

Implements edm::EDLooper.

Definition at line 218 of file Pi0FixedMassWindowCalibration.cc.

References abs, IslandClusterAlgo::barrel, barrelCells, barrelHits_, edm::SortedCollection< T, SORT >::begin(), ClusterShapeAlgo::Calculate(), calib, dtNoiseDBValidation_cfg::cerr, funct::cos(), gather_cfg::cout, DetId::Ecal, EcalBarrel, ecalHitsProducer_, EcalPreshower, ecalRecHitBarrelCollection, EcalCondObjectContainer< T >::end(), edm::SortedCollection< T, SORT >::end(), cppFunctionSkipper::exception, create_public_lumi_plots::exp, EcalCondObjectContainer< T >::find(), geometry, edm::EventSetup::get(), CaloGeometry::getValidDetIds(), i, EBDetId::ieta(), EBDetId::iphi(), isfirstcall_, island_p, j, gen::k, edm::EDLooperBase::kContinue, create_public_lumi_plots::log, IslandClusterAlgo::makeClusters(), mwxtals, nevent, newCalibs_barl, oldCalibs_barl, edm::ESHandle< class >::product(), DetId::rawId(), recalibEcalRecHitCollection, recHitsEB_map, selePi0DetaBelt_, selePi0DRBelt_, selePi0EtBeltIsoRatioMax_, selePi0MinvMeanFixed_, selePi0MinvSigmaFixed_, selePi0PtGammaOneMin_, selePi0PtPi0Min_, selePi0S4S9GammaOneMin_, selePi0S4S9GammaTwoMin_, selePi0S9S25GammaOneMin_, selePi0S9S25GammaTwoMin_, shapeAlgo_, funct::sin(), EcalCondObjectContainer< T >::size(), edm::SortedCollection< T, SORT >::size(), mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, funct::tan(), wxtals, and EBDetId::zside().

 {
   using namespace edm;
   using namespace std;
 
   // this chunk used to belong to beginJob(isetup). Moved here
   // with the beginJob without arguments migration
   
   if (isfirstcall_){
    // initialize arrays
 
     for (int sign=0; sign<2; sign++) {
       for (int ieta=0; ieta<85; ieta++) {
     for (int iphi=0; iphi<360; iphi++) {
       oldCalibs_barl[ieta][iphi][sign]=0.;
       newCalibs_barl[ieta][iphi][sign]=0.;
       wxtals[ieta][iphi][sign]=0.;
       mwxtals[ieta][iphi][sign]=0.;
     }
       }
     }
     
     // get initial constants out of DB
     
     edm::ESHandle<EcalIntercalibConstants> pIcal;
     EcalIntercalibConstantMap imap;
     
     try {
       setup.get<EcalIntercalibConstantsRcd>().get(pIcal);
       std::cout << "Taken EcalIntercalibConstants" << std::endl;
       imap = pIcal.product()->getMap();
       std::cout << "imap.size() = " << imap.size() << std::endl;
     } catch ( std::exception& ex ) {     
       std::cerr << "Error! can't get EcalIntercalibConstants " << std::endl;
     } 
     
     // get the ecal geometry:
     edm::ESHandle<CaloGeometry> geoHandle;
     setup.get<CaloGeometryRecord>().get(geoHandle);
     const CaloGeometry& geometry = *geoHandle;
     //const CaloSubdetectorGeometry *barrelGeometry = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
     
     // loop over all barrel crystals
     barrelCells = geometry.getValidDetIds(DetId::Ecal, EcalBarrel);
     std::vector<DetId>::const_iterator barrelIt;
     for (barrelIt=barrelCells.begin(); barrelIt!=barrelCells.end(); barrelIt++) {
       EBDetId eb(*barrelIt);
       
       // get the initial calibration constants
       EcalIntercalibConstantMap::const_iterator itcalib = imap.find(eb.rawId());
       if ( itcalib == imap.end() ) {
     // FIXME -- throw error
       }
       EcalIntercalibConstant calib = (*itcalib);
       int sign = eb.zside()>0 ? 1 : 0;
       oldCalibs_barl[abs(eb.ieta())-1][eb.iphi()-1][sign] = calib;
       if (eb.iphi()==1) std::cout << "Read old constant for crystal "
                   << " (" << eb.ieta() << "," << eb.iphi()
                   << ") : " << calib << std::endl;
       
       
       
     }
     isfirstcall_=false;
   }
   
 
 
   nevent++;
 
   if ((nevent<100 && nevent%10==0) 
       ||(nevent<1000 && nevent%100==0) 
       ||(nevent<10000 && nevent%100==0) 
       ||(nevent<100000 && nevent%1000==0) 
       ||(nevent<10000000 && nevent%1000==0))
     std::cout << "[Pi0FixedMassWindowCalibration] Events processed: "<<nevent<<std::endl;
 
   recHitsEB_map = new std::map<DetId, EcalRecHit>();
 
   EcalRecHitCollection* recalibEcalRecHitCollection( new EcalRecHitCollection );
 
   int nRecHitsEB=0;
   Handle<EcalRecHitCollection> pEcalRecHitBarrelCollection;
   event.getByLabel(ecalHitsProducer_, barrelHits_, pEcalRecHitBarrelCollection);
   const EcalRecHitCollection* ecalRecHitBarrelCollection = pEcalRecHitBarrelCollection.product();
   cout << " ECAL Barrel RecHits # "<< ecalRecHitBarrelCollection->size() <<endl;
   for(EcalRecHitCollection::const_iterator aRecHitEB = ecalRecHitBarrelCollection->begin(); aRecHitEB != ecalRecHitBarrelCollection->end(); aRecHitEB++) {
     //cout << " ECAL Barrel RecHit #,E,time,det,subdetid: "<<nRecHitsEB<<" "<<aRecHitEB->energy()<<" "<<aRecHitEB->time()<<" "<<aRecHitEB->detid().det()<<" "<<aRecHitEB->detid().subdetId()<<endl;
 
 
     EBDetId ebrhdetid = aRecHitEB->detid();
     //cout << " EBDETID: z,ieta,iphi "<<ebrhdetid.zside()<<" "<<ebrhdetid.ieta()<<" "<<ebrhdetid.iphi()<<endl;
     //cout << " EBDETID: tower_ieta,tower_iphi "<<ebrhdetid.tower_ieta()<<" "<<ebrhdetid.tower_iphi()<<endl;
     //cout << " EBDETID: iSM, ic "<<ebrhdetid.ism()<<" "<<ebrhdetid.ic()<<endl;
 
     int sign = ebrhdetid.zside()>0 ? 1 : 0;
     EcalRecHit aHit(aRecHitEB->id(),aRecHitEB->energy()*oldCalibs_barl[abs(ebrhdetid.ieta())-1][ebrhdetid.iphi()-1][sign],aRecHitEB->time());
     recalibEcalRecHitCollection->push_back(aHit);
 
     nRecHitsEB++;
 
   }
 
   //  cout<<" Recalib size: "<<recalibEcalRecHitCollection->size()<<endl;
   int irecalib=0;
   for(EcalRecHitCollection::const_iterator aRecHitEB = recalibEcalRecHitCollection->begin(); aRecHitEB != recalibEcalRecHitCollection->end(); aRecHitEB++) {
     //cout << " [recalibrated] ECAL Barrel RecHit #,E,time,det,subdetid: "<<irecalib<<" "<<aRecHitEB->energy()<<" "<<aRecHitEB->time()<<" "<<aRecHitEB->detid().det()<<" "<<aRecHitEB->detid().subdetId()<<endl;
 
     //    EBDetId ebrhdetid = aRecHitEB->detid();
     //cout << " [recalibrated] EBDETID: z,ieta,iphi "<<ebrhdetid.zside()<<" "<<ebrhdetid.ieta()<<" "<<ebrhdetid.iphi()<<endl;
     //cout << " [recalibrated] EBDETID: tower_ieta,tower_iphi "<<ebrhdetid.tower_ieta()<<" "<<ebrhdetid.tower_iphi()<<endl;
     //cout << " [recalibrated] EBDETID: iSM, ic "<<ebrhdetid.ism()<<" "<<ebrhdetid.ic()<<endl;
 
     std::pair<DetId, EcalRecHit> map_entry(aRecHitEB->id(), *aRecHitEB);
     recHitsEB_map->insert(map_entry);
 
 
     irecalib++;
 
   }
 
   
   // get the geometry and topology from the event setup:
   edm::ESHandle<CaloGeometry> geoHandle;
   setup.get<CaloGeometryRecord>().get(geoHandle);
 
   const CaloSubdetectorGeometry *geometry_p;
   CaloSubdetectorTopology *topology_p;
 
   std::string clustershapetag;
   geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
   topology_p = new EcalBarrelTopology(geoHandle);
 
   const CaloSubdetectorGeometry *geometryES_p;
   geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
 
   /*
   reco::BasicClusterCollection clusters;
   clusters = island_p->makeClusters(ecalRecHitBarrelCollection,geometry_p,topology_p,geometryES_p,IslandClusterAlgo::barrel);
   
   //Create associated ClusterShape objects.
   std::vector <reco::ClusterShape> ClusVec;
   for (int erg=0;erg<int(clusters.size());++erg){
     reco::ClusterShape TestShape = shapeAlgo_.Calculate(clusters[erg],ecalRecHitBarrelCollection,geometry_p,topology_p);
     ClusVec.push_back(TestShape);
   }
 
   //Put clustershapes in event, but retain a Handle on them.
   std::auto_ptr<reco::ClusterShapeCollection> clustersshapes_p(new reco::ClusterShapeCollection);
   clustersshapes_p->assign(ClusVec.begin(), ClusVec.end());
 
   cout<<"[Pi0Calibration] Basic Cluster collection size: "<<clusters.size()<<endl;
   cout<<"[Pi0Calibration] Basic Cluster Shape Collection size: "<<clustersshapes_p->size()<<endl;
 
   int iClus=0;
   for(reco::BasicClusterCollection::const_iterator aClus = clusters.begin(); aClus != clusters.end(); aClus++) {
     cout<<" CLUSTER : #,NHits,e,et,eta,phi,e2x2,e3x3,e5x5: "<<iClus<<" "<<aClus->getHitsByDetId().size()<<" "<<aClus->energy()<<" "<<aClus->energy()*sin(aClus->position().theta())<<" "<<aClus->position().eta()<<" "<<aClus->position().phi()<<" "<<(*clustersshapes_p)[iClus].e2x2()<<" "<<(*clustersshapes_p)[iClus].e3x3()<<" "<<(*clustersshapes_p)[iClus].e5x5()<<endl; 
     iClus++;
   }
   */
 
   // recalibrated clusters
   reco::BasicClusterCollection clusters_recalib;
   clusters_recalib = island_p->makeClusters(recalibEcalRecHitCollection,geometry_p,topology_p,geometryES_p,IslandClusterAlgo::barrel);
   
   //Create associated ClusterShape objects.
   std::vector <reco::ClusterShape> ClusVec_recalib;
   for (int erg=0;erg<int(clusters_recalib.size());++erg){
     reco::ClusterShape TestShape_recalib = shapeAlgo_.Calculate(clusters_recalib[erg],recalibEcalRecHitCollection,geometry_p,topology_p);
     ClusVec_recalib.push_back(TestShape_recalib);
   }
 
   //Put clustershapes in event, but retain a Handle on them.
   std::auto_ptr<reco::ClusterShapeCollection> clustersshapes_p_recalib(new reco::ClusterShapeCollection);
   clustersshapes_p_recalib->assign(ClusVec_recalib.begin(), ClusVec_recalib.end());
 
   cout<<"[Pi0FixedMassWindowCalibration][recalibration] Basic Cluster collection size: "<<clusters_recalib.size()<<endl;
   cout<<"[Pi0FixedMassWindowCalibration][recalibration] Basic Cluster Shape Collection size: "<<clustersshapes_p_recalib->size()<<endl;
 
 
   // pizero selection
 
   // Get ECAL Barrel Island Basic Clusters collection
   // ECAL Barrel Island Basic Clusters 
   static const int MAXBCEB = 200;
   static const int MAXBCEBRH = 200;
   int nIslandBCEB;
   float eIslandBCEB[MAXBCEB];
   float etIslandBCEB[MAXBCEB];
   float etaIslandBCEB[MAXBCEB];
   float phiIslandBCEB[MAXBCEB];
   float e2x2IslandBCEB[MAXBCEB];
   float e3x3IslandBCEB[MAXBCEB];
   float e5x5IslandBCEB[MAXBCEB];
   // indexes to the RecHits assiciated with 
   // ECAL Barrel Island Basic Clusters
   int nIslandBCEBRecHits[MAXBCEB];
   //  int indexIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
   int ietaIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
   int iphiIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
   int zsideIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
   float eIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
 
   nIslandBCEB=0;
   for(int i=0; i<MAXBCEB; i++){
     eIslandBCEB[i] = 0;
     etIslandBCEB[i] = 0;
     etaIslandBCEB[i] = 0;
     phiIslandBCEB[i] = 0;
     e2x2IslandBCEB[i] = 0;
     e3x3IslandBCEB[i] = 0;
     e5x5IslandBCEB[i] = 0;
      nIslandBCEBRecHits[i] = 0;
      for(int j=0;j<MAXBCEBRH;j++){
        //       indexIslandBCEBRecHits[i][j] = 0;
        ietaIslandBCEBRecHits[i][j] = 0;
        iphiIslandBCEBRecHits[i][j] = 0;
        zsideIslandBCEBRecHits[i][j] = 0;
        eIslandBCEBRecHits[i][j] = 0;
      }
   }
 
 
   int iClus_recalib=0;
   for(reco::BasicClusterCollection::const_iterator aClus = clusters_recalib.begin(); aClus != clusters_recalib.end(); aClus++) {
     cout<<" CLUSTER [recalibration] : #,NHits,e,et,eta,phi,e2x2,e3x3,e5x5: "<<iClus_recalib<<" "<<aClus->size()<<" "<<aClus->energy()<<" "<<aClus->energy()*sin(aClus->position().theta())<<" "<<aClus->position().eta()<<" "<<aClus->position().phi()<<" "<<(*clustersshapes_p_recalib)[iClus_recalib].e2x2()<<" "<<(*clustersshapes_p_recalib)[iClus_recalib].e3x3()<<" "<<(*clustersshapes_p_recalib)[iClus_recalib].e5x5()<<endl; 
 
     eIslandBCEB[nIslandBCEB] = aClus->energy();
     etIslandBCEB[nIslandBCEB] = aClus->energy()*sin(aClus->position().theta());
     etaIslandBCEB[nIslandBCEB] = aClus->position().eta();
     phiIslandBCEB[nIslandBCEB] = aClus->position().phi();
     
     e2x2IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e2x2();    
     e3x3IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e3x3();    
     e5x5IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e5x5();    
 
     nIslandBCEBRecHits[nIslandBCEB] = aClus->size();
 
     std::vector<std::pair< DetId,float> > hits = aClus->hitsAndFractions();
     std::vector<std::pair< DetId,float> >::iterator hit;
     std::map<DetId, EcalRecHit>::iterator aHit;
 
     int irhcount=0;
     for(hit = hits.begin(); hit != hits.end(); hit++)
       {
         // need to get hit by DetID in order to get energy
         aHit = recHitsEB_map->find((*hit).first);
         //cout << "       RecHit #: "<<irhcount<<" from Basic Cluster with Energy: "<<aHit->second.energy()<<endl; 
 
         EBDetId sel_rh = aHit->second.detid();
         //cout << "       RecHit: z,ieta,iphi "<<sel_rh.zside()<<" "<<sel_rh.ieta()<<" "<<sel_rh.iphi()<<endl;
         //cout << "       RecHit: tower_ieta,tower_iphi "<<sel_rh.tower_ieta()<<" "<<sel_rh.tower_iphi()<<endl;
         //cout << "       RecHit: iSM, ic "<<sel_rh.ism()<<" "<<sel_rh.ic()<<endl;
 
     ietaIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.ieta();
     iphiIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.iphi();
     zsideIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.zside();
     eIslandBCEBRecHits[nIslandBCEB][irhcount] = aHit->second.energy();
 
         irhcount++;
       }
     nIslandBCEB++;
     iClus_recalib++;
 
 
   }
 
   // Selection, based on ECAL Barrel Basic Clusters
 
   if (nIslandBCEB > 1) 
     {
       for(int i=0 ; i<nIslandBCEB ; i++)
         {
           for(int j=i+1 ; j<nIslandBCEB ; j++)
             {
 
               if( etIslandBCEB[i]>selePi0PtGammaOneMin_ && etIslandBCEB[j]>selePi0PtGammaOneMin_) 
         {
                 
           float theta_0 = 2. * atan(exp(-etaIslandBCEB[i]));
           float theta_1 = 2. * atan(exp(-etaIslandBCEB[j]));
                 
           float p0x = eIslandBCEB[i] * sin(theta_0) * cos(phiIslandBCEB[i]);
           float p1x = eIslandBCEB[j] * sin(theta_1) * cos(phiIslandBCEB[j]);
                 
           float p0y = eIslandBCEB[i] * sin(theta_0) * sin(phiIslandBCEB[i]);
           float p1y = eIslandBCEB[j] * sin(theta_1) * sin(phiIslandBCEB[j]);
                 
           float p0z = eIslandBCEB[i] * cos(theta_0);
           float p1z = eIslandBCEB[j] * cos(theta_1);
 
           float pi0_px = p0x + p1x;
           float pi0_py = p0y + p1y;
           float pi0_pz = p0z + p1z;
                 
           float pi0_ptot = sqrt (pi0_px*pi0_px + pi0_py*pi0_py + pi0_pz*pi0_pz); 
 
           float pi0_theta = acos(pi0_pz/pi0_ptot);
           float pi0_eta = -log(tan(pi0_theta/2));
           float pi0_phi = atan(pi0_py/pi0_px);
           //cout << " pi0_theta, pi0_eta, pi0_phi "<<pi0_theta<<" "<<pi0_eta<<" "<<pi0_phi<<endl;
 
           // belt isolation
                 
           float et_belt=0;
           for(Int_t k=0 ; k<nIslandBCEB ; k++)
             {
               if( (k != i) && (k != j) ) 
             {
               float dr_pi0_k = sqrt( (etaIslandBCEB[k]-pi0_eta)*(etaIslandBCEB[k]-pi0_eta) + (phiIslandBCEB[k]-pi0_phi)*(phiIslandBCEB[k]-pi0_phi) );
               float deta_pi0_k = fabs(etaIslandBCEB[k]-pi0_eta);
               if ( (dr_pi0_k<selePi0DRBelt_) && (deta_pi0_k<selePi0DetaBelt_) ) et_belt = et_belt + etIslandBCEB[k];
             }
             }
 
                 
           float pt_pi0 = sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
           //float dr_pi0 = sqrt ( (etaIslandBCEB[i]-etaIslandBCEB[j])*(etaIslandBCEB[i]-etaIslandBCEB[j]) + (phiIslandBCEB[i]-phiIslandBCEB[j])*(phiIslandBCEB[i]-phiIslandBCEB[j]) );
 
           //cout <<" pi0 pt,dr:  "<<pt_pi0<<" "<<dr_pi0<<endl; 
           if (pt_pi0 > selePi0PtPi0Min_) 
             {
               float m_inv = sqrt ( (eIslandBCEB[i] + eIslandBCEB[j])*(eIslandBCEB[i] + eIslandBCEB[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );  
               cout <<" pi0 (pt>2.5 GeV) m_inv = "<<m_inv<<endl;  
 
               float s4s9_1 = e2x2IslandBCEB[i]/e3x3IslandBCEB[i];
               float s4s9_2 = e2x2IslandBCEB[j]/e3x3IslandBCEB[j];
 
               float s9s25_1 = e3x3IslandBCEB[i]/e5x5IslandBCEB[i];
               float s9s25_2 = e3x3IslandBCEB[j]/e5x5IslandBCEB[j];
 
               //float s9Esc_1 = e3x3IslandBCEB[i]/eIslandBCEB[i];
               //float s9Esc_2 = e3x3IslandBCEB[j]/eIslandBCEB[j];
 
               if (s4s9_1>selePi0S4S9GammaOneMin_ && s4s9_2>selePi0S4S9GammaTwoMin_ && s9s25_1>selePi0S9S25GammaOneMin_ && s9s25_2>selePi0S9S25GammaTwoMin_ && (et_belt/pt_pi0 < selePi0EtBeltIsoRatioMax_) ) 
             {
               //good pizero candidate
               if ( m_inv>(selePi0MinvMeanFixed_ - 2*selePi0MinvSigmaFixed_) && m_inv<(selePi0MinvMeanFixed_ + 2*selePi0MinvSigmaFixed_) ) 
                 {
                   //fill wxtals and mwxtals weights
                   cout<<" Pi0 Good candidate : minv = "<<m_inv<<endl;
                   for(int kk=0 ; kk<nIslandBCEBRecHits[i] ; kk++)
                 {
                   int ieta_xtal = ietaIslandBCEBRecHits[i][kk];
                   int iphi_xtal = iphiIslandBCEBRecHits[i][kk];
                   int sign = zsideIslandBCEBRecHits[i][kk]>0 ? 1 : 0;
                   wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] = wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] + eIslandBCEBRecHits[i][kk]/e3x3IslandBCEB[i];
                   mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] = mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] + (selePi0MinvMeanFixed_/m_inv)*(selePi0MinvMeanFixed_/m_inv)*(eIslandBCEBRecHits[i][kk]/e3x3IslandBCEB[i]);
                   cout<< "[Pi0FixedMassWindowCalibration] eta, phi, sign, e, e3x3, wxtals and mwxtals: "<<ieta_xtal<<" "<<iphi_xtal<<" "<<sign<<" "<<eIslandBCEBRecHits[i][kk]<<" "<<e3x3IslandBCEB[i]<<" "<<wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign]<<" "<<mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign]<<endl;
                 }
 
                   for(int kk=0 ; kk<nIslandBCEBRecHits[j] ; kk++)
                 {
                   int ieta_xtal = ietaIslandBCEBRecHits[j][kk];
                   int iphi_xtal = iphiIslandBCEBRecHits[j][kk];
                   int sign = zsideIslandBCEBRecHits[j][kk]>0 ? 1 : 0;
                   wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] = wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] + eIslandBCEBRecHits[j][kk]/e3x3IslandBCEB[j];
                   mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] = mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign] + (selePi0MinvMeanFixed_/m_inv)*(selePi0MinvMeanFixed_/m_inv)*(eIslandBCEBRecHits[j][kk]/e3x3IslandBCEB[j]);
                   cout<< "[Pi0FixedMassWindowCalibration] eta, phi, sign, e, e3x3, wxtals and mwxtals: "<<ieta_xtal<<" "<<iphi_xtal<<" "<<sign<<" "<<eIslandBCEBRecHits[j][kk]<<" "<<e3x3IslandBCEB[j]<<" "<<wxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign]<<" "<<mwxtals[abs(ieta_xtal)-1][iphi_xtal-1][sign]<<endl;
                 }
 
                 }
             }
             }
         }
                 
           
         } // End of the "j" loop over BCEB
     } // End of the "i" loop over BCEB
 
 
     } else {
       cout<< " Not enough ECAL Barrel Basic Clusters: "<<nIslandBCEB<<endl;
     }
 
 
 
 
 
 
   return kContinue;
 }

void Pi0FixedMassWindowCalibration::endOfJob ( )

virtual

Called at end of job.

Reimplemented from edm::EDLooperBase.

Definition at line 119 of file Pi0FixedMassWindowCalibration.cc.

References abs, barrelCells, gather_cfg::cout, EcalBarrel, EBDetId::ieta(), EBDetId::iphi(), newCalibs_barl, Pi0CalibXMLwriter::writeLine(), and EBDetId::zside().

 {
   
   std::cout << "[Pi0FixedMassWindowCalibration] endOfJob"<<endl;
 
   // Write new calibration constants
 
   Pi0CalibXMLwriter barrelWriter(EcalBarrel,99);
 
   std::vector<DetId>::const_iterator barrelIt=barrelCells.begin();
   for (; barrelIt!=barrelCells.end(); barrelIt++) {
     EBDetId eb(*barrelIt);
     int ieta = eb.ieta();
     int iphi = eb.iphi();
     int sign = eb.zside()>0 ? 1 : 0;
     barrelWriter.writeLine(eb,newCalibs_barl[abs(ieta)-1][iphi-1][sign]);
   }
 
 }

edm::EDLooper::Status Pi0FixedMassWindowCalibration::endOfLoop	(	const edm::EventSetup &	iSetup,
		unsigned int	iLoop
	)

virtual

Called at end of loop.

Implements edm::EDLooperBase.

Definition at line 161 of file Pi0FixedMassWindowCalibration.cc.

References abs, barrelCells, gather_cfg::cout, EcalBarrel, EBDetId::ieta(), EBDetId::iphi(), edm::EDLooperBase::kContinue, edm::EDLooperBase::kStop, mwxtals, newCalibs_barl, oldCalibs_barl, theMaxLoops, Pi0CalibXMLwriter::writeLine(), wxtals, and EBDetId::zside().

 {
   std::cout << "[Pi0FixedMassWindowCalibration] Ending loop " << iLoop<<std::endl;
 
   for (int sign=0; sign<2; sign++) {
     for (int ieta=0; ieta<85; ieta++) {
       for (int iphi=0; iphi<360; iphi++) {
 
         if (wxtals[ieta][iphi][sign] == 0 ) 
       {
         newCalibs_barl[ieta][iphi][sign]=oldCalibs_barl[ieta][iphi][sign];
       } else {
         newCalibs_barl[ieta][iphi][sign]=oldCalibs_barl[ieta][iphi][sign]*(mwxtals[ieta][iphi][sign]/wxtals[ieta][iphi][sign]);
       }
     cout<< " New calibration constant: ieta iphi sign - old,mwxtals,wxtals,new: "<<ieta<<" "<<iphi<<" "<<sign<<" - "<<oldCalibs_barl[ieta][iphi][sign]<<" "<<mwxtals[ieta][iphi][sign]<<" "<<wxtals[ieta][iphi][sign]<<" "<<newCalibs_barl[ieta][iphi][sign]<<endl;
 
       }
     }
   }
     
   Pi0CalibXMLwriter barrelWriter(EcalBarrel,iLoop+1);
 
   std::vector<DetId>::const_iterator barrelIt=barrelCells.begin();
   for (; barrelIt!=barrelCells.end(); barrelIt++) {
     EBDetId eb(*barrelIt);
     int ieta = eb.ieta();
     int iphi = eb.iphi();
     int sign = eb.zside()>0 ? 1 : 0;
     barrelWriter.writeLine(eb,newCalibs_barl[abs(ieta)-1][iphi-1][sign]);
     if (iphi==1) {
       std::cout << "Calib constant for barrel crystal "
         << " (" << ieta << "," << iphi << ") changed from "
         << oldCalibs_barl[abs(ieta)-1][iphi-1][sign] << " to "
         << newCalibs_barl[abs(ieta)-1][iphi-1][sign] << std::endl;
     }
   }
 
   // replace old calibration constants with new one
 
   for (int sign=0; sign<2; sign++) {
     for (int ieta=0; ieta<85; ieta++) {
       for (int iphi=0; iphi<360; iphi++) {
         oldCalibs_barl[ieta][iphi][sign]=newCalibs_barl[ieta][iphi][sign];
     newCalibs_barl[ieta][iphi][sign]=0;
       }
     }
   }
 
 
   if ( iLoop == theMaxLoops-1 || iLoop >= theMaxLoops ) return kStop;
   else return kContinue;
 }

virtual void Pi0FixedMassWindowCalibration::produce	(	edm::Event &	,
		const edm::EventSetup &
	)

inlinevirtual

Dummy implementation (job done in duringLoop)

Definition at line 70 of file Pi0FixedMassWindowCalibration.h.

Referenced by JSONExport.JsonExport::export(), HTMLExport.HTMLExport::export(), and HTMLExport.HTMLExportStatic::export().

70 {};

void Pi0FixedMassWindowCalibration::startingNewLoop ( unsigned int iLoop )

virtual

Called at beginning of loop.

Implements edm::EDLooperBase.

Definition at line 141 of file Pi0FixedMassWindowCalibration.cc.

References gather_cfg::cout, mwxtals, and wxtals.

 {
 
   for (int sign=0; sign<2; sign++) {
     for (int ieta=0; ieta<85; ieta++) {
       for (int iphi=0; iphi<360; iphi++) {
     wxtals[ieta][iphi][sign]=0.;
     mwxtals[ieta][iphi][sign]=0.;
       }
     }
   }
   std::cout << "[Pi0FixedMassWindowCalibration] Starting loop number " << iLoop<<std::endl;
 
 }