#include <Calibration/AlCaHOCalibProducer/src/AlCaHOCalibProducer.cc>

Inheritance diagram for AlCaHOCalibProducer:

Public Types
typedef Basic3DVector< float >	DirectionType

typedef Basic3DVector< float >	PositionType

typedef Basic3DVector< float >	RotationType

Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

Public Types inherited from edm::ProducerBase
using	ModuleToResolverIndicies = std::unordered_multimap< std::string, std::tuple< edm::TypeID const , const char , edm::ProductResolverIndex >>

typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Public Member Functions
	AlCaHOCalibProducer (const edm::ParameterSet &)

	~AlCaHOCalibProducer () override

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

SerialTaskQueue *	globalLuminosityBlocksQueue ()

SerialTaskQueue *	globalRunsQueue ()

ModuleDescription const &	moduleDescription () const

	~EDProducer () override

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

std::vector< edm::ProductResolverIndex > const &	indiciesForPutProducts (BranchType iBranchType) const

	ProducerBase ()

std::vector< edm::ProductResolverIndex > const &	putTokenIndexToProductResolverIndex () const

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

std::function< void(BranchDescription const &)>	registrationCallback () const
	used by the fwk to register list of products More...

void	resolvePutIndicies (BranchType iBranchType, ModuleToResolverIndicies const &iIndicies, std::string const &moduleLabel)

	~ProducerBase () noexcept(false) 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

ESProxyIndex const *	esGetTokenIndices (edm::Transition iTrans) const

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)

void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Private Types
typedef math::Error< 5 >::type	CovarianceMatrix

Private Member Functions
void	beginJob () override

void	beginLuminosityBlock (const edm::LuminosityBlock &, const edm::EventSetup &) override

void	endJob () override

void	fillHOStore (const reco::TrackRef &ncosm, HOCalibVariables &tmpHOCalib, std::unique_ptr< HOCalibVariableCollection > &hostore, int Noccu_old, int indx, edm::Handle< reco::TrackCollection > cosmicmuon, edm::View< reco::Muon >::const_iterator muon1, const edm::Event &iEvent, const edm::EventSetup &iSetup)

void	findHOEtaPhi (int iphsect, int &ietaho, int &iphiho)

FreeTrajectoryState	getFreeTrajectoryState (const reco::Track &tk, const MagneticField *field, int itag, bool dir)

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

Private Attributes
std::map< std::string, bool >	fired

edm::ESHandle< HcalSeverityLevelComputer >	hcalSevLvlComputerHndl

TH2F *	ho_occupency [5]

int	iring

float	localxhor0

float	localxhor1

float	localyhor0

float	localyhor1

bool	m_cosmic

int	m_endTS

bool	m_hbinfo

bool	m_occupancy

double	m_sigma

int	m_startTS

edm::InputTag	muonTags_

const int	ncidmx = 5

const int	netabin = 16

const int	netamx = 32

int	Nevents

int	Noccu

const int	nphimx = 72

int	nRuns

unsigned int	Ntp

const double	rHOL0 = 382.0

const double	rHOL1 = 407.0

const HcalChannelQuality *	theHcalChStatus

edm::EDGetTokenT< HBHERecHitCollection >	tok_hbhe_

edm::EDGetTokenT< HORecHitCollection >	tok_ho_

edm::EDGetTokenT< LumiScalersCollection >	tok_lumi_

edm::EDGetTokenT< edm::View< reco::Muon > >	tok_muons_

edm::EDGetTokenT< reco::TrackCollection >	tok_muonsCosmic_

edm::EDGetTokenT< CaloTowerCollection >	tok_tower_

edm::EDGetTokenT< reco::VertexCollection >	tok_vertex_

float	xhor0

float	xhor1

float	yhor0

float	yhor1

Additional Inherited Members
Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

static bool	wantsGlobalLuminosityBlocks ()

static bool	wantsGlobalRuns ()

static bool	wantsStreamLuminosityBlocks ()

static bool	wantsStreamRuns ()

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 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

change magnetic field inside ../data/HOCosmicCalib_RecoLocalMuon.cff ../data/HOCosmicCalib_RecoLocalTracker.cff

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>=""> Missing towers : eta=5, phi=18-19 : eta = -5, phi =11-14

HO tile sizes Ring +-2 : width Tray 6:404.6, 5&4:347.6, 3:352.6, 2:364.6, 1:315.6 (phi ordering is opposite) lenght Tile 1:420.1, 2:545.1, 3:583.3, 4:626.0, 5:335.5

    (five tiles, 1 is close to Ring 1 and 5 is towardslc endcap)

Ring +-1 : width Tray 6:404.6, 5&4:347.6, 3:352.6, 2:364.6, 1:315.6 (same as Ring+-2) lenght Tile 1:391.5, 2:394.2, 3:411.0, 4:430.9, 5:454.0, 6:426.0 (1: near R0 and 6 near R2)

Ring 0 L1 : Width Tray (6:290.6, 5&4:345.6, 3:350.6, 2:362.6, 1:298.6 lenght 1:351.2, 2:353.8, 3:359.2, 4:189.1 (4 is towards Ring1)

Ring 0 L0 : Width Tray 6:266.6, 5&4:325.6, 3:330.6, 2:341.6, 1:272.6 length 1:331.5, 2:334.0, 3:339.0, 4:248.8 (4 is towards Ring1)

Definition at line 152 of file AlCaHOCalibProducer.cc.

Member Typedef Documentation

typedef math::Error<5>::type AlCaHOCalibProducer::CovarianceMatrix

private

Definition at line 220 of file AlCaHOCalibProducer.cc.

typedef Basic3DVector<float> AlCaHOCalibProducer::DirectionType

Definition at line 158 of file AlCaHOCalibProducer.cc.

typedef Basic3DVector<float> AlCaHOCalibProducer::PositionType

Definition at line 157 of file AlCaHOCalibProducer.cc.

typedef Basic3DVector<float> AlCaHOCalibProducer::RotationType

Definition at line 159 of file AlCaHOCalibProducer.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 248 of file AlCaHOCalibProducer.cc.

References edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), ho_occupency, m_cosmic, m_hbinfo, m_occupancy, m_sigma, TFileService::make(), muonTags_, netamx, nphimx, fftjetcommon_cfi::title, tok_hbhe_, tok_ho_, tok_lumi_, tok_muons_, tok_muonsCosmic_, tok_tower_, and tok_vertex_.

                                                                        {
   //register your products
 
   m_hbinfo = iConfig.getUntrackedParameter<bool>("hbinfo", false);
   m_sigma = iConfig.getUntrackedParameter<double>("sigma", 0.05);
   m_occupancy = iConfig.getUntrackedParameter<bool>("plotOccupancy", false);
   m_cosmic = iConfig.getUntrackedParameter<bool>("CosmicData", false);
 
   // keep InputTag muonTags_ since it is used below. - cowden
   muonTags_ = iConfig.getUntrackedParameter<edm::InputTag>("muons");
   tok_muonsCosmic_ = consumes<reco::TrackCollection>(muonTags_);
   tok_muons_ = consumes<edm::View<reco::Muon> >(muonTags_);
   tok_vertex_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexTags"));
   //  tok_lumi_ = consumes<LumiDetails ,edm::InLumi>(iConfig.getParameter<edm::InputTag>("lumiTags"));
   tok_lumi_ = consumes<LumiScalersCollection>(iConfig.getParameter<edm::InputTag>("lumiTags"));
   tok_ho_ = consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInput"));
   tok_hbhe_ = consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInput"));
   tok_tower_ = consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("towerInput"));
 
   produces<HOCalibVariableCollection>("HOCalibVariableCollection").setBranchAlias("HOCalibVariableCollection");
 
   if (m_occupancy) {
     edm::Service<TFileService> fs;
 
     char title[200];
 
     for (int ij = 0; ij < 5; ij++) {
       sprintf(title, "ho_occupency (>%i #sigma)", ij + 2);
       ho_occupency[ij] =
           fs->make<TH2F>(title, title, netamx + 1, -netamx - 0.5, netamx / 2 + 0.5, nphimx, 0.5, nphimx + 0.5);
     }
   }
 }

AlCaHOCalibProducer::~AlCaHOCalibProducer ( )

override

Definition at line 282 of file AlCaHOCalibProducer.cc.

                                           {
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)
 }

Member Function Documentation

void AlCaHOCalibProducer::beginJob ( void )

overrideprivatevirtual

Reimplemented from edm::EDProducer.

Definition at line 369 of file AlCaHOCalibProducer.cc.

References Nevents, Noccu, and nRuns.

                                    {
   Nevents = 0;
   nRuns = 0;
   Noccu = 0;
 }

void AlCaHOCalibProducer::beginLuminosityBlock	(	const edm::LuminosityBlock &	lBlock,
		const edm::EventSetup &	es
	)

overrideprivate

Definition at line 390 of file AlCaHOCalibProducer.cc.

References edm::EventSetup::get(), edm::ESHandle< T >::product(), and theHcalChStatus.

                                                                                                         {
   edm::ESHandle<HcalChannelQuality> hcalChStatus;
   es.get<HcalChannelQualityRcd>().get("withTopo", hcalChStatus);
   theHcalChStatus = hcalChStatus.product();
 
   return;
 }

void AlCaHOCalibProducer::endJob ( void )

overrideprivatevirtual

Reimplemented from edm::EDProducer.

Definition at line 376 of file AlCaHOCalibProducer.cc.

References ho_occupency, m_occupancy, SiStripPI::max, Nevents, and Noccu.

Referenced by o2olib.O2ORunMgr::executeJob().

                                  {
   if (m_occupancy) {
     for (int ij = 0; ij < 5; ij++) {
       ho_occupency[ij]->Scale(1. / std::max(1, Noccu));
     }
   }
   edm::LogInfo("HOCalib") << " AlCaHOCalibProducer processed event " << Nevents;
 }

void AlCaHOCalibProducer::fillHOStore	(	const reco::TrackRef &	ncosm,
		HOCalibVariables &	tmpHOCalib,
		std::unique_ptr< HOCalibVariableCollection > &	hostore,
		int	Noccu_old,
		int	indx,
		edm::Handle< reco::TrackCollection >	cosmicmuon,
		edm::View< reco::Muon >::const_iterator	muon1,
		const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

private

Definition at line 398 of file AlCaHOCalibProducer.cc.

Referenced by produce().

                                                                    {
   edm::ESHandle<CaloGeometry> pG;
   iSetup.get<CaloGeometryRecord>().get(pG);
   const CaloGeometry* geo = pG.product();
   const CaloSubdetectorGeometry* gHO = geo->getSubdetectorGeometry(DetId::Hcal, HcalOuter);
 
   // Get Hcal Severity Level Computer, so that the severity of each rechit flag/status may be determined
   //  edm::ESHandle<HcalSeverityLevelComputer> hcalSevLvlComputerHndl;
   iSetup.get<HcalSeverityLevelComputerRcd>().get(hcalSevLvlComputerHndl);
   const HcalSeverityLevelComputer* hcalSevLvlComputer = hcalSevLvlComputerHndl.product();
 
   int charge = ncosm->charge();
 
   double innerr = (*ncosm).innerPosition().Perp2();
   double outerr = (*ncosm).outerPosition().Perp2();
   int iiner = (innerr < outerr) ? 1 : 0;
 
   //---------------------------------------------------
   //             in_to_out  Dir         in_to_out  Dir
   //   StandAlone ^         ^     Cosmic    ^    |
   //              |         |               |    v
   //---------------------------------------------------Y=0
   //   StandAlone |         |     Cosmic    ^    |
   //              v         v               |    v
   //----------------------------------------------------
 
   double posx, posy, posz;
   double momx, momy, momz;
 
   if (iiner == 1) {
     posx = (*ncosm).innerPosition().X();
     posy = (*ncosm).innerPosition().Y();
     posz = (*ncosm).innerPosition().Z();
 
     momx = (*ncosm).innerMomentum().X();
     momy = (*ncosm).innerMomentum().Y();
     momz = (*ncosm).innerMomentum().Z();
 
   } else {
     posx = (*ncosm).outerPosition().X();
     posy = (*ncosm).outerPosition().Y();
     posz = (*ncosm).outerPosition().Z();
 
     momx = (*ncosm).outerMomentum().X();
     momy = (*ncosm).outerMomentum().Y();
     momz = (*ncosm).outerMomentum().Z();
   }
 
   PositionType trkpos(posx, posy, posz);
 
   CLHEP::Hep3Vector tmpmuon3v(posx, posy, posz);
   CLHEP::Hep3Vector tmpmuondir(momx, momy, momz);
 
   bool samedir = (tmpmuon3v.dot(tmpmuondir) > 0) ? true : false;
   for (int ij = 0; ij < 3; ij++) {
     tmpHOCalib.caloen[ij] = 0.0;
   }
   int inearbymuon = 0;
   localxhor0 = localyhor0 = 20000;  //GM for 22OCT07 data
 
   if (m_cosmic) {
     int ind(0);
     for (reco::TrackCollection::const_iterator ncosmcor = cosmicmuon->begin(); ncosmcor != cosmicmuon->end();
          ++ncosmcor, ++ind) {
       if (indx == ind)
         continue;
       CLHEP::Hep3Vector tmpmuon3vcor;
       CLHEP::Hep3Vector tmpmom3v;
       if (iiner == 1) {
         tmpmuon3vcor = CLHEP::Hep3Vector(
             (*ncosmcor).innerPosition().X(), (*ncosmcor).innerPosition().Y(), (*ncosmcor).innerPosition().Z());
         tmpmom3v = CLHEP::Hep3Vector(
             (*ncosmcor).innerMomentum().X(), (*ncosmcor).innerMomentum().Y(), (*ncosmcor).innerMomentum().Z());
       } else {
         tmpmuon3vcor = CLHEP::Hep3Vector(
             (*ncosmcor).outerPosition().X(), (*ncosmcor).outerPosition().Y(), (*ncosmcor).outerPosition().Z());
         tmpmom3v = CLHEP::Hep3Vector(
             (*ncosmcor).outerMomentum().X(), (*ncosmcor).outerMomentum().Y(), (*ncosmcor).outerMomentum().Z());
       }
 
       if (tmpmom3v.mag() < 0.2 || (*ncosmcor).ndof() < 5)
         continue;
 
       double angle = tmpmuon3v.angle(tmpmuon3vcor);
       if (angle < 7.5 * CLHEP::deg) {
         inearbymuon = 1;
       }  //  break;}
 
       //    if (muonTagsi_.label() =="cosmicMuons") {
       if (angle < 7.5 * CLHEP::deg) {
         tmpHOCalib.caloen[0] += 1.;
       }
       if (angle < 15.0 * CLHEP::deg) {
         tmpHOCalib.caloen[1] += 1.;
       }
       if (angle < 35.0 * CLHEP::deg) {
         tmpHOCalib.caloen[2] += 1.;
       }
     }
   } else {
     //            if (muonTags_.label() =="muons") {
     edm::Handle<CaloTowerCollection> calotower;
     iEvent.getByToken(tok_tower_, calotower);
 
     for (CaloTowerCollection::const_iterator calt = calotower->begin(); calt != calotower->end(); calt++) {
       //CMSSW_2_1_x const math::XYZVector towermom = (*calt).momentum();
       double ith = (*calt).momentum().theta();
       double iph = (*calt).momentum().phi();
 
       CLHEP::Hep3Vector calo3v(sin(ith) * cos(iph), sin(ith) * sin(iph), cos(ith));
 
       double angle = tmpmuon3v.angle(calo3v);
 
       if (angle < 7.5 * CLHEP::deg) {
         tmpHOCalib.caloen[0] += calt->emEnergy() + calt->hadEnergy();
       }
       if (angle < 15 * CLHEP::deg) {
         tmpHOCalib.caloen[1] += calt->emEnergy() + calt->hadEnergy();
       }
       if (angle < 35 * CLHEP::deg) {
         tmpHOCalib.caloen[2] += calt->emEnergy() + calt->hadEnergy();
       }
     }
   }
   if ((m_cosmic) || (tmpHOCalib.caloen[0] <= 10.0)) {
     GlobalPoint glbpt(posx, posy, posz);
 
     double mom = sqrt(momx * momx + momy * momy + momz * momz);
 
     momx /= mom;
     momy /= mom;
     momz /= mom;
 
     DirectionType trkdir(momx, momy, momz);
 
     tmpHOCalib.trkdr = (*ncosm).d0();
     tmpHOCalib.trkdz = (*ncosm).dz();
     tmpHOCalib.nmuon = (m_cosmic) ? cosmicmuon->size() : 1;
     tmpHOCalib.trkvx = glbpt.x();
     tmpHOCalib.trkvy = glbpt.y();
     tmpHOCalib.trkvz = glbpt.z();
     tmpHOCalib.trkmm = mom * charge;
     tmpHOCalib.trkth = trkdir.theta();
     tmpHOCalib.trkph = trkdir.phi();
     tmpHOCalib.isect2 = -2;
     tmpHOCalib.isect = -2;
     tmpHOCalib.hodx = -100;
     tmpHOCalib.hody = -100;
     tmpHOCalib.hoang = -2.0;
     tmpHOCalib.momatho = -2;
     tmpHOCalib.ndof = (inearbymuon == 0) ? (int)(*ncosm).ndof() : -(int)(*ncosm).ndof();
     tmpHOCalib.chisq = (*ncosm).normalizedChi2();  // max(1.,tmpHOCalib.ndof);
     if (!m_cosmic) {
       reco::MuonEnergy muonenr = muon1->calEnergy();
       reco::MuonIsolation iso03 = muon1->isolationR03();
       reco::MuonIsolation iso05 = muon1->isolationR05();
 
       tmpHOCalib.tkpt03 = iso03.sumPt;
       tmpHOCalib.ecal03 = iso05.sumPt;  // iso03.emEt+muonenr.em;
       tmpHOCalib.hcal03 = iso03.hadEt + muonenr.had;
     }
     tmpHOCalib.therr = 0.;
     tmpHOCalib.pherr = 0.;
     if (iiner == 1) {
       reco::TrackBase::CovarianceMatrix innercov = (*ncosm).innerStateCovariance();
       tmpHOCalib.therr = innercov(1, 1);  //thetaError();
       tmpHOCalib.pherr = innercov(2, 2);  //phi0Error();
     } else {
       reco::TrackBase::CovarianceMatrix outercov = (*ncosm).outerStateCovariance();
       tmpHOCalib.therr = outercov(1, 1);  //thetaError();
       tmpHOCalib.pherr = outercov(2, 2);  //phi0Error();
     }
     edm::ESHandle<MagneticField> theMagField;
     iSetup.get<IdealMagneticFieldRecord>().get(theMagField);
 
     SteppingHelixPropagator myHelix(&*theMagField, anyDirection);
     myHelix.setMaterialMode(false);
     myHelix.applyRadX0Correction(true);
     double phiho = trkpos.phi();
     if (phiho < 0)
       phiho += CLHEP::twopi;
 
     int iphisect_dt = int(6 * (phiho + 10.0 * CLHEP::deg) / CLHEP::pi);  //for u 18/12/06
     if (iphisect_dt >= 12)
       iphisect_dt = 0;
 
     int iphisect = -1;
     bool ipath = false;
     for (int kl = 0; kl <= 2; kl++) {
       int iphisecttmp = (kl < 2) ? iphisect_dt + kl : iphisect_dt - 1;
       if (iphisecttmp < 0)
         iphisecttmp = 11;
       if (iphisecttmp >= 12)
         iphisecttmp = 0;
 
       double phipos = iphisecttmp * CLHEP::pi / 6.;
       double phirot = phipos;
 
       GlobalVector xLocal(-sin(phirot), cos(phirot), 0.);
       GlobalVector yLocal(0., 0., 1.);
       GlobalVector zLocal = xLocal.cross(yLocal).unit();
       //    GlobalVector zLocal(cos(phirot), sin(phirot), 0.0);
 
       FreeTrajectoryState freetrajectorystate_ = getFreeTrajectoryState(*ncosm, &(*theMagField), iiner, samedir);
 
       Surface::RotationType rot(xLocal, yLocal, zLocal);
 
       for (int ik = 1; ik >= 0; ik--) {  //propagate track in two HO layers
 
         double radial = rHOL1;
         if (ik == 0)
           radial = rHOL0;
 
         Surface::PositionType pos(radial * cos(phipos), radial * sin(phipos), 0.);
         PlaneBuilder::ReturnType aPlane = PlaneBuilder().plane(pos, rot);
 
         auto aPlane2 = new Plane(pos, rot);
 
         SteppingHelixStateInfo steppingHelixstateinfo_;
         myHelix.propagate(SteppingHelixStateInfo(freetrajectorystate_), (*aPlane2), steppingHelixstateinfo_);
 
         if (steppingHelixstateinfo_.isValid()) {
           GlobalPoint hotrkpos2xx(steppingHelixstateinfo_.position().x(),
                                   steppingHelixstateinfo_.position().y(),
                                   steppingHelixstateinfo_.position().z());
 
           if (ik == 1) {
             HcalDetId ClosestCell = (HcalDetId)gHO->getClosestCell(hotrkpos2xx);
             int ixeta = ClosestCell.ieta();
             int ixphi = ClosestCell.iphi();
             tmpHOCalib.isect2 = 100 * std::abs(ixeta + 50) + std::abs(ixphi);
           }
 
           GlobalVector hotrkpos2(steppingHelixstateinfo_.position().x(),
                                  steppingHelixstateinfo_.position().y(),
                                  steppingHelixstateinfo_.position().z());
           CLHEP::Hep3Vector hotrkdir2(steppingHelixstateinfo_.momentum().x(),
                                       steppingHelixstateinfo_.momentum().y(),
                                       steppingHelixstateinfo_.momentum().z());
 
           LocalVector lclvt0 = (*aPlane).toLocal(hotrkpos2);
 
           double xx = lclvt0.x();
           double yy = lclvt0.y();
 
           if (ik == 1) {
             if ((std::abs(yy) < 130 && xx > -64.7 && xx < 138.2)                              //Ring-0
                 || (std::abs(yy) > 130 && std::abs(yy) < 700 && xx > -76.3 && xx < 140.5)) {  //Ring +-1,2
               ipath = true;  //Only look for tracks which as hits in layer 1
               iphisect = iphisecttmp;
             }
           }
 
           if (iphisect != iphisecttmp)
             continue;  //Look for ring-0 only when ring1 is accepted for that sector
 
           switch (ik) {
             case 0:
               xhor0 = xx;  //lclvt0.x();
               yhor0 = yy;  //lclvt0.y();
               break;
             case 1:
               xhor1 = xx;  //lclvt0.x();
               yhor1 = yy;  //lclvt0.y();
               tmpHOCalib.momatho = hotrkdir2.mag();
               tmpHOCalib.hoang = CLHEP::Hep3Vector(zLocal.x(), zLocal.y(), zLocal.z()).dot(hotrkdir2.unit());
               break;
             default:
               break;
           }
         } else {
           break;
         }
       }
       if (ipath)
         break;
     }
     if (ipath) {  //If muon crossed HO laeyrs
 
       int ietaho = 50;
       int iphiho = -1;
 
       for (int ij = 0; ij < 9; ij++) {
         tmpHOCalib.hosig[ij] = -100.0;
       }
       for (int ij = 0; ij < 18; ij++) {
         tmpHOCalib.hocorsig[ij] = -100.0;
       }
       for (int ij = 0; ij < 9; ij++) {
         tmpHOCalib.hbhesig[ij] = -100.0;
       }
       tmpHOCalib.hocro = -100;
       tmpHOCalib.htime = -1000;
 
       int isect = 0;
 
       findHOEtaPhi(iphisect, ietaho, iphiho);
 
       if (ietaho != 0 && iphiho != 0 && std::abs(iring) <= 2) {  //Muon passed through a tower
         isect = 100 * std::abs(ietaho + 50) + std::abs(iphiho);
         if (std::abs(ietaho) >= netabin || iphiho < 0)
           isect *= -1;  //Not extrapolated to any tower
         if (std::abs(ietaho) >= netabin)
           isect -= 1000000;  //not matched with eta
         if (iphiho < 0)
           isect -= 2000000;  //not matched with phi
         tmpHOCalib.isect = isect;
 
         tmpHOCalib.hodx = localxhor1;
         tmpHOCalib.hody = localyhor1;
 
         if (iring == 0) {
           tmpHOCalib.hocorsig[8] = localxhor0;
           tmpHOCalib.hocorsig[9] = localyhor0;
         }
 
         int etamn = -4;
         int etamx = 4;
         if (iring == 1) {
           etamn = 5;
           etamx = 10;
         }
         if (iring == 2) {
           etamn = 11;
           etamx = 16;
         }
         if (iring == -1) {
           etamn = -10;
           etamx = -5;
         }
         if (iring == -2) {
           etamn = -16;
           etamx = -11;
         }
 
         int phimn = 1;
         int phimx = 2;
         if (iring == 0) {
           phimx = 2 * int((iphiho + 1) / 2.);
           phimn = phimx - 1;
         } else {
           phimn = 3 * int((iphiho + 1) / 3.) - 1;
           phimx = phimn + 2;
         }
 
         if (phimn < 1)
           phimn += nphimx;
         if (phimx > 72)
           phimx -= nphimx;
 
         if (m_hbinfo) {
           for (int ij = 0; ij < 9; ij++) {
             tmpHOCalib.hbhesig[ij] = -100.0;
           }
 
           edm::Handle<HBHERecHitCollection> hbheht;  // iEvent.getByType(hbheht);
           iEvent.getByToken(tok_hbhe_, hbheht);
 
           if (!(*hbheht).empty()) {
             if ((*hbheht).empty())
               throw(int)(*hbheht).size();
 
             for (HBHERecHitCollection::const_iterator jk = (*hbheht).begin(); jk != (*hbheht).end(); jk++) {
               HcalDetId id = (*jk).id();
               int tmpeta = id.ieta();
               int tmpphi = id.iphi();
 
               int deta = tmpeta - ietaho;
               if (tmpeta < 0 && ietaho > 0)
                 deta += 1;
               if (tmpeta > 0 && ietaho < 0)
                 deta -= 1;
 
               //        if (tmpeta==-1 && ietaho== 1) deta = -1;
               //        if (tmpeta== 1 && ietaho==-1) deta =  1;
 
               int dphi = tmpphi - iphiho;
               if (dphi > nphimx / 2) {
                 dphi -= nphimx;
               }
               if (dphi < -nphimx / 2) {
                 dphi += nphimx;
               }
 
               //        if (phimn >phimx) {
               //          if (dphi==71) dphi=-1;
               //          if (dphi==-71) dphi=1;
               //        }
 
               if (m_occupancy) {
                 float signal = (*jk).energy();
                 //      int tmpeta1 = (tmpeta>0) ? tmpeta -1 : -tmpeta +14;
                 if (signal > -100 && Noccu == Noccu_old) {
                   for (int ij = 0; ij < 5; ij++) {
                     if (signal > (ij + 2) * m_sigma) {
                       ho_occupency[ij]->Fill(tmpeta, tmpphi);
                     }
                   }
                 }
               }
 
               int ipass2 = (std::abs(deta) <= 1 && std::abs(dphi) <= 1) ? 1 : 0;  //NEED correction in full CMS detector
               if (ipass2 == 0)
                 continue;
 
               float signal = (*jk).energy();
 
               if (3 * (deta + 1) + dphi + 1 < 9)
                 tmpHOCalib.hbhesig[3 * (deta + 1) + dphi + 1] = signal;
             }
           }
         }  //m_hbinfo #endif
 
         edm::Handle<HORecHitCollection> hoht;
         iEvent.getByToken(tok_ho_, hoht);
 
         if (!(*hoht).empty()) {
           for (HORecHitCollection::const_iterator jk = (*hoht).begin(); jk != (*hoht).end(); jk++) {
             HcalDetId id = (*jk).id();
             int tmpeta = id.ieta();
             int tmpphi = id.iphi();
 
             int ipass1 = 0;
             if (tmpeta >= etamn && tmpeta <= etamx) {
               if (phimn < phimx) {
                 ipass1 = (tmpphi >= phimn && tmpphi <= phimx) ? 1 : 0;
               } else {
                 ipass1 = (tmpphi == 71 || tmpphi == 72 || tmpphi == 1) ? 1 : 0;
               }
             }
 
             int deta = tmpeta - ietaho;
             int dphi = tmpphi - iphiho;
 
             if (tmpeta < 0 && ietaho > 0)
               deta += 1;
             if (tmpeta > 0 && ietaho < 0)
               deta -= 1;
             //        if (tmpeta==-1 && ietaho== 1) deta = -1;
             //        if (tmpeta== 1 && ietaho==-1) deta =  1;
 
             if (dphi > nphimx / 2) {
               dphi -= nphimx;
             }
             if (dphi < -nphimx / 2) {
               dphi += nphimx;
             }
             //        if (phimn>phimx) {
             //      if (dphi==71) dphi=-1;
             //      if (dphi==-71) dphi=1;
             //        }
 
             float signal = (*jk).energy();
 
             int ipass2 = (std::abs(deta) <= 1 && std::abs(dphi) <= 1) ? 1 : 0;
 
             if (ipass1 == 0 && ipass2 == 0)
               continue;
 
             if (ipass1 == 1) {
               int tmpdph = tmpphi - phimn;
               if (tmpdph < 0)
                 tmpdph = 2;  //only case of iphi==1, where phimn=71
 
               int ilog = 2 * (tmpeta - etamn) + tmpdph;
               if (iring != 0) {
                 if (iring > 0) {
                   ilog = 3 * (tmpeta - etamn) + tmpdph;  //Again CMS correction
                 } else {
                   ilog = 3 * (etamx - tmpeta) + tmpdph;  //Again CMS correction
                 }
               }
               if (ilog > -1 && ilog < 18) {
                 tmpHOCalib.hocorsig[ilog] = signal;
               }
             }
 
             if (ipass2 == 1) {
               if (3 * (deta + 1) + dphi + 1 < 9) {
                 tmpHOCalib.hosig[3 * (deta + 1) + dphi + 1] = signal;  //Again CMS azimuthal near phi 1&72
               }
             }
 
             if (deta == 0 && dphi == 0) {
               tmpHOCalib.htime = (*jk).time();
               tmpHOCalib.hoflag = (*jk).flags();
 
               // Get Channel Quality information for the given detID
               unsigned theStatusValue = theHcalChStatus->getValues(id)->getValue();
               // Now get severity of problems for the given detID, based on the rechit flag word and the channel quality status value
               int hitSeverity = hcalSevLvlComputer->getSeverityLevel(id, (*jk).flags(), theStatusValue);
               tmpHOCalib.hoflag = hitSeverity;
               int crphi = tmpphi + 6;
               if (crphi > 72)
                 crphi -= 72;
 
               for (HORecHitCollection::const_iterator jcr = (*hoht).begin(); jcr != (*hoht).end(); jcr++) {
                 const HORecHit reccr = (const HORecHit)(*jcr);
                 HcalDetId idcr = reccr.id();
                 int etacr = idcr.ieta();
                 int phicr = idcr.iphi();
                 if (tmpeta == etacr && crphi == phicr) {
                   tmpHOCalib.hocro = reccr.energy();
                 }
               }
             }
           }
         }
       }
 
       //GMA   Npass++;
       if (Noccu == Noccu_old)
         Noccu++;
       hostore->push_back(tmpHOCalib);
     }  // if (ipath)
   }    // Cut on calo energy
 }

void AlCaHOCalibProducer::findHOEtaPhi	(	int	iphsect,
		int &	ietaho,
		int &	iphiho
	)

private

Definition at line 924 of file AlCaHOCalibProducer.cc.

References funct::abs(), iring, localxhor0, localxhor1, localyhor0, localyhor1, netabin, nphimx, xhor0, xhor1, yhor0, and yhor1.

Referenced by fillHOStore().

                                                                              {
   //18/12/06 : use only position, not angle phi
 
   const double etalow[16] = {0.025,
                              35.195,
                              70.625,
                              106.595,
                              141.565,
                              180.765,
                              220.235,
                              261.385,
                              304.525,
                              349.975,
                              410.025,
                              452.085,
                              506.645,
                              565.025,
                              627.725,
                              660.25};
   const double etahgh[16] = {35.145,
                              70.575,
                              106.545,
                              125.505,
                              180.715,
                              220.185,
                              261.335,
                              304.475,
                              349.925,
                              392.575,
                              452.035,
                              506.595,
                              564.975,
                              627.675,
                              661.075,
                              700.25};
 
   const double philow[6] = {-76.27, -35.11, 0.35, 35.81, 71.77, 108.93};  //Ring+/-1 & 2
   const double phihgh[6] = {-35.81, -0.35, 35.11, 71.07, 108.23, 140.49};
 
   const double philow00[6] = {-60.27, -32.91, 0.35, 33.61, 67.37, 102.23};  //Ring0 L0
   const double phihgh00[6] = {-33.61, -0.35, 32.91, 66.67, 101.53, 129.49};
 
   const double philow01[6] = {-64.67, -34.91, 0.35, 35.61, 71.37, 108.33};  //Ring0 L1
   const double phihgh01[6] = {-35.61, -0.35, 34.91, 70.67, 107.63, 138.19};
 
   iring = -10;
 
   double tmpdy = std::abs(yhor1);
   for (int ij = 0; ij < netabin; ij++) {
     if (tmpdy > etalow[ij] && tmpdy < etahgh[ij]) {
       ietaho = ij + 1;
       float tmp1 = fabs(tmpdy - etalow[ij]);
       float tmp2 = fabs(tmpdy - etahgh[ij]);
 
       localyhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
       if (yhor1 < 0)
         localyhor1 *= -1.;
 
       if (ij < 4)
         iring = 0;
       if (ij >= 4 && ij < 10)
         iring = 1;
       if (ij >= 10 && ij < netabin)
         iring = 2;
       break;
     }
   }
 
   int tmpphi = 0;
   int tmpphi0 = 0;
 
   if (ietaho > 4) {  //Ring 1 and 2
     for (int ij = 0; ij < 6; ij++) {
       if (xhor1 > philow[ij] && xhor1 < phihgh[ij]) {
         tmpphi = ij + 1;
         float tmp1 = fabs(xhor1 - philow[ij]);
         float tmp2 = fabs(xhor1 - phihgh[ij]);
         localxhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
         break;
       }
     }
   } else {  //Ring 0
     for (int ij = 0; ij < 6; ij++) {
       if (xhor1 > philow01[ij] && xhor1 < phihgh01[ij]) {
         tmpphi = ij + 1;
         float tmp1 = fabs(xhor1 - philow01[ij]);
         float tmp2 = fabs(xhor1 - phihgh01[ij]);
         localxhor1 = (tmp1 < tmp2) ? -tmp1 : tmp2;
         break;
       }
     }
 
     for (int ij = 0; ij < 6; ij++) {
       if (xhor0 > philow00[ij] && xhor0 < phihgh00[ij]) {
         tmpphi0 = ij + 1;
         float tmp1 = fabs(xhor0 - philow00[ij]);
         float tmp2 = fabs(xhor0 - phihgh00[ij]);
         localxhor0 = (tmp1 < tmp2) ? -tmp1 : tmp2;
         if (tmpphi != tmpphi0)
           localxhor0 += 10000.;
         break;
       }
     }
 
     double tmpdy = std::abs(yhor0);
     for (int ij = 0; ij < 4; ij++) {
       if (tmpdy > etalow[ij] && tmpdy < etahgh[ij]) {
         float tmp1 = fabs(tmpdy - etalow[ij]);
         float tmp2 = fabs(tmpdy - etahgh[ij]);
         localyhor0 = (tmp1 < tmp2) ? -tmp1 : tmp2;
         if (yhor0 < 0)
           localyhor0 *= -1.;
         if (ij + 1 != ietaho)
           localyhor0 += 10000.;
         break;
       }
     }
   }
 
   if (tmpphi != 0) {
     iphiho = 6 * iphisect - 2 + tmpphi;
     if (iphiho <= 0)
       iphiho += nphimx;
     if (iphiho > nphimx)
       iphiho -= nphimx;
   }
 
   //  isect2 = 15*iring+iphisect+1;
 
   if (yhor1 < 0) {
     if (std::abs(ietaho) > netabin) {  //Initialised with 50
       ietaho += 1;
     } else {
       ietaho *= -1;
     }
     //    isect2 *=-1;
     iring *= -1;
   }
 }

FreeTrajectoryState AlCaHOCalibProducer::getFreeTrajectoryState	(	const reco::Track &	tk,
		const MagneticField *	field,
		int	itag,
		bool	dir
	)

private

Definition at line 1064 of file AlCaHOCalibProducer.cc.

References reco::TrackBase::charge(), DEFINE_FWK_MODULE, reco::Track::extra(), reco::Track::innerMomentum(), reco::Track::innerPosition(), reco::Track::outerPx(), reco::Track::outerPy(), reco::Track::outerPz(), reco::Track::outerX(), reco::Track::outerY(), and reco::Track::outerZ().

Referenced by fillHOStore().

                                                                           {
   if (iiner == 0) {
     GlobalPoint gpos(tk.outerX(), tk.outerY(), tk.outerZ());
     GlobalVector gmom(tk.outerPx(), tk.outerPy(), tk.outerPz());
     if (dir)
       gmom *= -1.;
     GlobalTrajectoryParameters par(gpos, gmom, tk.charge(), field);
     CurvilinearTrajectoryError err(tk.extra()->outerStateCovariance());
     return FreeTrajectoryState(par, err);
   } else {
     GlobalPoint gpos(tk.innerPosition().X(), tk.innerPosition().Y(), tk.innerPosition().Z());
     GlobalVector gmom(tk.innerMomentum().X(), tk.innerMomentum().Y(), tk.innerMomentum().Z());
     if (dir)
       gmom *= -1.;
     GlobalTrajectoryParameters par(gpos, -gmom, tk.charge(), field);
     CurvilinearTrajectoryError err(tk.extra()->innerStateCovariance());
     return FreeTrajectoryState(par, err);
   }
 }

void AlCaHOCalibProducer::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

overrideprivate

Definition at line 292 of file AlCaHOCalibProducer.cc.

References edm::EventID::event(), fillHOStore(), edm::Event::getByToken(), edm::EventBase::id(), HOCalibVariables::inslumi, edm::EventBase::isRealData(), edm::HandleBase::isValid(), m_cosmic, eostools::move(), Nevents, Noccu, HOCalibVariables::nprim, jets_cff::primaryVertices, edm::Event::put(), edm::EventID::run(), tok_lumi_, tok_muons_, tok_muonsCosmic_, and tok_vertex_.

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

                                                                                {
   int irun = iEvent.id().run();
   //  int ilumi = iEvent.luminosityBlock();
 
   Nevents++;
 
   if (Nevents % 5000 == 1)
     edm::LogInfo("HOCalib") << "AlCaHOCalibProducer Processing event # " << Nevents << " " << Noccu << " " << irun
                             << " " << iEvent.id().event();
 
   auto hostore = std::make_unique<HOCalibVariableCollection>();
 
   edm::Handle<reco::TrackCollection> cosmicmuon;
   edm::Handle<edm::View<reco::Muon> > collisionmuon;
 
   bool muonOK(true);
   HOCalibVariables tmpHOCalib;
   tmpHOCalib.nprim = -1;
   tmpHOCalib.inslumi = -1.;
 
   if (m_cosmic) {
     iEvent.getByToken(tok_muonsCosmic_, cosmicmuon);
     muonOK = (cosmicmuon.isValid() && !cosmicmuon->empty());
   } else {
     iEvent.getByToken(tok_muons_, collisionmuon);
     muonOK = (collisionmuon.isValid() && !collisionmuon->empty());
 
     if (iEvent.isRealData()) {
       edm::Handle<reco::VertexCollection> primaryVertices;
       iEvent.getByToken(tok_vertex_, primaryVertices);
       if (primaryVertices.isValid()) {
         tmpHOCalib.nprim = primaryVertices->size();
       }
 
       tmpHOCalib.inslumi = 0.;
 
       edm::Handle<LumiScalersCollection> lumiScale;
       iEvent.getByToken(tok_lumi_, lumiScale);
 
       if (lumiScale.isValid()) {
         if (lumiScale->empty()) {
           edm::LogError("HOCalib") << "lumiScale collection is empty";
         } else {
           tmpHOCalib.inslumi = lumiScale->begin()->pileup();
         }
       }
     }
   }
 
   if (muonOK) {
     int Noccu_old = Noccu;
     edm::View<reco::Muon>::const_iterator muon1;
     if (m_cosmic) {
       int indx(0);
       for (reco::TrackCollection::const_iterator ncosm = cosmicmuon->begin(); ncosm != cosmicmuon->end();
            ++ncosm, ++indx) {
         if ((*ncosm).ndof() < 15)
           continue;
         if ((*ncosm).normalizedChi2() > 30.0)
           continue;
         reco::TrackRef tRef = reco::TrackRef(cosmicmuon, indx);
         fillHOStore(tRef, tmpHOCalib, hostore, Noccu_old, indx, cosmicmuon, muon1, iEvent, iSetup);
       }
     } else {
       for (muon1 = collisionmuon->begin(); muon1 < collisionmuon->end(); muon1++) {
         if ((!muon1->isGlobalMuon()) || (!muon1->isTrackerMuon()))
           continue;
         reco::TrackRef ncosm = muon1->innerTrack();
         fillHOStore(ncosm, tmpHOCalib, hostore, Noccu_old, 0, cosmicmuon, muon1, iEvent, iSetup);
       }
     }
   }
 
   iEvent.put(std::move(hostore), "HOCalibVariableCollection");
 }