#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::one::EDProducerBase
typedef EDProducerBase	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=default

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

	EDProducer (const EDProducer &)=delete

SerialTaskQueue *	globalLuminosityBlocksQueue () final

SerialTaskQueue *	globalRunsQueue () final

bool	hasAbilityToProduceInBeginLumis () const final

bool	hasAbilityToProduceInBeginProcessBlocks () const final

bool	hasAbilityToProduceInBeginRuns () const final

bool	hasAbilityToProduceInEndLumis () const final

bool	hasAbilityToProduceInEndProcessBlocks () const final

bool	hasAbilityToProduceInEndRuns () const final

const EDProducer &	operator= (const EDProducer &)=delete

bool	wantsGlobalLuminosityBlocks () const final

bool	wantsGlobalRuns () const final

bool	wantsInputProcessBlocks () const final

bool	wantsProcessBlocks () const final

Public Member Functions inherited from edm::one::EDProducerBase
	EDProducerBase ()

ModuleDescription const &	moduleDescription () const

bool	wantsStreamLuminosityBlocks () const

bool	wantsStreamRuns () const

	~EDProducerBase () 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

std::vector< bool > const &	recordProvenanceList () 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)

TypeLabelList const &	typeLabelList () const
	used by the fwk to register the list of products of this module More...

	~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

std::vector< ESProxyIndex > const &	esGetTokenIndicesVector (edm::Transition iTrans) const

std::vector< ESRecordIndex > const &	esGetTokenRecordIndicesVector (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::array< std::vector< ModuleDescription const > , NumBranchTypes > &modulesAll, std::vector< ModuleProcessName > &modulesInPreviousProcesses, 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

void	selectInputProcessBlocks (ProductRegistry const &productRegistry, ProcessBlockHelperBase const &processBlockHelperBase)

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

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

void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Static Public Member Functions
static void	fillDescriptions (edm::ConfigurationDescriptions &descriptions)

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

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

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

Private Member Functions
void	beginJob () 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 CaloSubdetectorGeometry *, const MagneticField &)

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

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

const HcalSeverityLevelComputer *	theHcalSevLvlComputer

edm::ESGetToken< CaloGeometry, CaloGeometryRecord >	tok_geom_

edm::EDGetTokenT< HBHERecHitCollection >	tok_hbhe_

edm::ESGetToken< HcalChannelQuality, HcalChannelQualityRcd >	tok_hcalChStatus_

edm::ESGetToken< HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd >	tok_hcalSevLvlComputer_

edm::EDGetTokenT< HORecHitCollection >	tok_ho_

edm::EDGetTokenT< LumiScalersCollection >	tok_lumi_

edm::ESGetToken< MagneticField, IdealMagneticFieldRecord >	tok_magField_

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
Protected Member Functions inherited from edm::ProducerBase
template<Transition Tr = Transition::Event>
auto	produces (std::string instanceName) noexcept
	declare what type of product will make and with which optional label More...

template<Transition B>
BranchAliasSetter	produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)

template<BranchType B>
BranchAliasSetter	produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)

BranchAliasSetter	produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)

template<typename ProductType , Transition B>
BranchAliasSetterT< ProductType >	produces (std::string instanceName)

template<class ProductType >
BranchAliasSetterT< ProductType >	produces ()

template<typename ProductType , BranchType B>
BranchAliasSetterT< ProductType >	produces (std::string instanceName)

template<typename ProductType , BranchType B>
BranchAliasSetterT< ProductType >	produces ()

template<class ProductType >
BranchAliasSetterT< ProductType >	produces (std::string instanceName)

template<typename ProductType , Transition B>
BranchAliasSetterT< ProductType >	produces ()

template<Transition Tr = Transition::Event>
auto	produces () noexcept

ProducesCollector	producesCollector ()

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

template<BranchType B = InEvent>
EDConsumerBaseAdaptor< B >	consumes (edm::InputTag tag) noexcept

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<Transition Tr = Transition::Event>
constexpr auto	esConsumes ()

template<Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag tag)

template<Transition Tr = Transition::Event>
ESGetTokenGeneric	esConsumes (eventsetup::EventSetupRecordKey const &iRecord, eventsetup::DataKey const &iKey)
	Used with EventSetupRecord::doGet. More...

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)

void	resetItemsToGetFrom (BranchType iType)

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 143 of file AlCaHOCalibProducer.cc.

Member Typedef Documentation

◆ CovarianceMatrix

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

private

Definition at line 217 of file AlCaHOCalibProducer.cc.

◆ DirectionType

typedef Basic3DVector<float> AlCaHOCalibProducer::DirectionType

Definition at line 150 of file AlCaHOCalibProducer.cc.

◆ PositionType

typedef Basic3DVector<float> AlCaHOCalibProducer::PositionType

Definition at line 149 of file AlCaHOCalibProducer.cc.

◆ RotationType

typedef Basic3DVector<float> AlCaHOCalibProducer::RotationType

Definition at line 151 of file AlCaHOCalibProducer.cc.

Constructor & Destructor Documentation

◆ AlCaHOCalibProducer()

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

explicit

Definition at line 244 of file AlCaHOCalibProducer.cc.

References compareTotals::fs, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), ho_occupency, TFileService::kSharedResource, m_cosmic, m_hbinfo, m_occupancy, m_sigma, muonTags_, netamx, nphimx, runGCPTkAlMap::title, tok_geom_, tok_hbhe_, tok_hcalChStatus_, tok_hcalSevLvlComputer_, tok_ho_, tok_lumi_, tok_magField_, tok_muons_, tok_muonsCosmic_, tok_tower_, and tok_vertex_.

                                                                        {
   usesResource(TFileService::kSharedResource);
   //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"));
 
   tok_hcalChStatus_ = esConsumes<HcalChannelQuality, HcalChannelQualityRcd>(edm::ESInputTag("", "withTopo"));
   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
   tok_hcalSevLvlComputer_ = esConsumes<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd>();
   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
 
   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::~AlCaHOCalibProducer ( )

overridedefault

Member Function Documentation

◆ beginJob()

void AlCaHOCalibProducer::beginJob ( )

overrideprivatevirtual

Reimplemented from edm::one::EDProducerBase.

Definition at line 393 of file AlCaHOCalibProducer.cc.

References Nevents, Noccu, and nRuns.

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

◆ endJob()

void AlCaHOCalibProducer::endJob ( void )

overrideprivatevirtual

Reimplemented from edm::one::EDProducerBase.

Definition at line 400 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;
 }

◆ fillDescriptions()

void AlCaHOCalibProducer::fillDescriptions ( edm::ConfigurationDescriptions & descriptions )

static

Definition at line 288 of file AlCaHOCalibProducer.cc.

References edm::ConfigurationDescriptions::add(), submitPVResolutionJobs::desc, HLT_2022v12_cff::InputTag, and AlCaHLTBitMon_QueryRunRegistry::string.

                                                                                      {
   edm::ParameterSetDescription desc;
   desc.add<edm::InputTag>("hbheInput", edm::InputTag("hbhereco"));
   desc.addUntracked<bool>("hotime", false);
   desc.addUntracked<bool>("hbinfo", false);
   desc.addUntracked<double>("sigma", 1.0);
   desc.addUntracked<bool>("plotOccupancy", false);
   desc.addUntracked<bool>("CosmicData", false);
   desc.add<edm::InputTag>("hoInput", edm::InputTag("horeco"));
   desc.add<edm::InputTag>("towerInput", edm::InputTag("towerMaker"));
   desc.addUntracked<std::string>("RootFileName", "test.root");
   desc.addUntracked<double>("m_scale", 4.0);
   desc.addUntracked<bool>("debug", false);
   desc.addUntracked<edm::InputTag>("muons", edm::InputTag("muons"));
   desc.add<edm::InputTag>("vertexTags", edm::InputTag("offlinePrimaryVertices"));
   desc.add<edm::InputTag>("lumiTags", edm::InputTag("scalersRawToDigi"));
   descriptions.add("alcaHOCalibProducer", desc);
 }

◆ fillHOStore()

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 CaloSubdetectorGeometry *	gHO,
		const MagneticField &	magField
	)

private

Definition at line 409 of file AlCaHOCalibProducer.cc.

Referenced by produce().

                                                                      {
   // Get Hcal Severity Level Computer, so that the severity of each rechit flag/status may be determined
 
   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") {
     auto const& calotower = iEvent.getHandle(tok_tower_);
 
     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();
     }
 
     SteppingHelixPropagator myHelix(&magField, 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, &(magField), 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;
           }
 
           auto const& hbheht = iEvent.getHandle(tok_hbhe_);  // iEvent.getByType(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
 
         auto const& hoht = iEvent.getHandle(tok_ho_);
 
         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 = theHcalSevLvlComputer->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
 }

◆ findHOEtaPhi()

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

private

Definition at line 922 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;
   }
 }

◆ getFreeTrajectoryState()

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

private

Definition at line 1062 of file AlCaHOCalibProducer.cc.

References reco::TrackBase::charge(), DeadROC_duringRun::dir, submitPVResolutionJobs::err, 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);
   }
 }

◆ produce()

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

overrideprivatevirtual

Implements edm::one::EDProducerBase.

Definition at line 308 of file AlCaHOCalibProducer.cc.

References fillHOStore(), edm::EventSetup::getData(), CaloGeometry::getSubdetectorGeometry(), DetId::Hcal, HcalOuter, iEvent, HOCalibVariables::inslumi, edm::HandleBase::isValid(), BXlumiParameters_cfi::lumiScale, m_cosmic, eostools::move(), Nevents, Noccu, HOCalibVariables::nprim, HLT_2022v12_cff::primaryVertices, theHcalChStatus, theHcalSevLvlComputer, tok_geom_, tok_hcalChStatus_, tok_hcalSevLvlComputer_, tok_lumi_, tok_magField_, tok_muons_, tok_muonsCosmic_, and tok_vertex_.

                                                                                {
   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();
 
   theHcalChStatus = &iSetup.getData(tok_hcalChStatus_);
 
   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) {
     cosmicmuon = iEvent.getHandle(tok_muonsCosmic_);
     muonOK = (cosmicmuon.isValid() && !cosmicmuon->empty());
   } else {
     collisionmuon = iEvent.getHandle(tok_muons_);
     muonOK = (collisionmuon.isValid() && !collisionmuon->empty());
 
     if (iEvent.isRealData()) {
       auto const& primaryVertices = iEvent.getHandle(tok_vertex_);
       if (primaryVertices.isValid()) {
         tmpHOCalib.nprim = primaryVertices->size();
       }
 
       tmpHOCalib.inslumi = 0.;
 
       auto const& lumiScale = iEvent.getHandle(tok_lumi_);
 
       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;
 
     theHcalSevLvlComputer = &iSetup.getData(tok_hcalSevLvlComputer_);
 
     MagneticField const& magField = iSetup.getData(tok_magField_);
 
     const CaloGeometry& geo = iSetup.getData(tok_geom_);
     const CaloSubdetectorGeometry* gHO = geo.getSubdetectorGeometry(DetId::Hcal, HcalOuter);
 
     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, gHO, magField);
       }
     } 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, gHO, magField);
       }
     }
   }
 
   iEvent.put(std::move(hostore), "HOCalibVariableCollection");
 }