9 #include "Math/Interpolator.h"
13 : theEBthreshold(-1000.),
14 theEEthreshold(-1000.),
16 theUseEtEBTresholdFlag(
false),
17 theUseEtEETresholdFlag(
false),
18 theUseSymEBTresholdFlag(
false),
19 theUseSymEETresholdFlag(
false),
22 theHcalThreshold(-1000.),
23 theHBthreshold(-1000.),
24 theHESthreshold(-1000.),
25 theHEDthreshold(-1000.),
26 theHOthreshold0(-1000.),
27 theHOthresholdPlus1(-1000.),
28 theHOthresholdMinus1(-1000.),
29 theHOthresholdPlus2(-1000.),
30 theHOthresholdMinus2(-1000.),
31 theHF1threshold(-1000.),
32 theHF2threshold(-1000.),
33 theEBGrid(std::vector<double>(5,10.)),
34 theEBWeights(std::vector<double>(5,1.)),
35 theEEGrid(std::vector<double>(5,10.)),
36 theEEWeights(std::vector<double>(5,1.)),
37 theHBGrid(std::vector<double>(5,10.)),
38 theHBWeights(std::vector<double>(5,1.)),
39 theHESGrid(std::vector<double>(5,10.)),
40 theHESWeights(std::vector<double>(5,1.)),
41 theHEDGrid(std::vector<double>(5,10.)),
42 theHEDWeights(std::vector<double>(5,1.)),
43 theHOGrid(std::vector<double>(5,10.)),
44 theHOWeights(std::vector<double>(5,1.)),
45 theHF1Grid(std::vector<double>(5,10.)),
46 theHF1Weights(std::vector<double>(5,1.)),
47 theHF2Grid(std::vector<double>(5,10.)),
48 theHF2Weights(std::vector<double>(5,1.)),
58 theEBSumThreshold(-1000.),
59 theEESumThreshold(-1000.),
62 theTowerConstituentsMap(0),
65 theMomConstrMethod(0),
77 bool useSymEBTreshold,
78 bool useSymEETreshold,
81 double HBthreshold,
double HESthreshold,
double HEDthreshold,
82 double HOthreshold0,
double HOthresholdPlus1,
double HOthresholdMinus1,
83 double HOthresholdPlus2,
double HOthresholdMinus2,
84 double HF1threshold,
double HF2threshold,
85 double EBweight,
double EEweight,
86 double HBweight,
double HESweight,
double HEDweight,
87 double HOweight,
double HF1weight,
double HF2weight,
97 : theEBthreshold(EBthreshold),
98 theEEthreshold(EEthreshold),
100 theUseEtEBTresholdFlag(useEtEBTreshold),
101 theUseEtEETresholdFlag(useEtEETreshold),
102 theUseSymEBTresholdFlag(useSymEBTreshold),
103 theUseSymEETresholdFlag(useSymEETreshold),
105 theHcalThreshold(HcalThreshold),
106 theHBthreshold(HBthreshold),
107 theHESthreshold(HESthreshold),
108 theHEDthreshold(HEDthreshold),
109 theHOthreshold0(HOthreshold0),
110 theHOthresholdPlus1(HOthresholdPlus1),
111 theHOthresholdMinus1(HOthresholdMinus1),
112 theHOthresholdPlus2(HOthresholdPlus2),
113 theHOthresholdMinus2(HOthresholdMinus2),
114 theHF1threshold(HF1threshold),
115 theHF2threshold(HF2threshold),
116 theEBGrid(std::vector<double>(5,10.)),
117 theEBWeights(std::vector<double>(5,1.)),
118 theEEGrid(std::vector<double>(5,10.)),
119 theEEWeights(std::vector<double>(5,1.)),
120 theHBGrid(std::vector<double>(5,10.)),
121 theHBWeights(std::vector<double>(5,1.)),
122 theHESGrid(std::vector<double>(5,10.)),
123 theHESWeights(std::vector<double>(5,1.)),
124 theHEDGrid(std::vector<double>(5,10.)),
125 theHEDWeights(std::vector<double>(5,1.)),
126 theHOGrid(std::vector<double>(5,10.)),
127 theHOWeights(std::vector<double>(5,1.)),
128 theHF1Grid(std::vector<double>(5,10.)),
129 theHF1Weights(std::vector<double>(5,1.)),
130 theHF2Grid(std::vector<double>(5,10.)),
131 theHF2Weights(std::vector<double>(5,1.)),
132 theEBweight(EBweight),
133 theEEweight(EEweight),
134 theHBweight(HBweight),
135 theHESweight(HESweight),
136 theHEDweight(HEDweight),
137 theHOweight(HOweight),
138 theHF1weight(HF1weight),
139 theHF2weight(HF2weight),
140 theEcutTower(EcutTower),
141 theEBSumThreshold(EBSumThreshold),
142 theEESumThreshold(EESumThreshold),
145 theMomConstrMethod(momConstrMethod),
146 theMomHBDepth(momHBDepth),
147 theMomHEDepth(momHEDepth),
148 theMomEBDepth(momEBDepth),
149 theMomEEDepth(momEEDepth)
156 bool useEtEBTreshold,
157 bool useEtEETreshold,
158 bool useSymEBTreshold,
159 bool useSymEETreshold,
162 double HBthreshold,
double HESthreshold,
double HEDthreshold,
163 double HOthreshold0,
double HOthresholdPlus1,
double HOthresholdMinus1,
164 double HOthresholdPlus2,
double HOthresholdMinus2,
165 double HF1threshold,
double HF2threshold,
174 double EBweight,
double EEweight,
175 double HBweight,
double HESweight,
double HEDweight,
176 double HOweight,
double HF1weight,
double HF2weight,
186 : theEBthreshold(EBthreshold),
187 theEEthreshold(EEthreshold),
189 theUseEtEBTresholdFlag(useEtEBTreshold),
190 theUseEtEETresholdFlag(useEtEETreshold),
191 theUseSymEBTresholdFlag(useSymEBTreshold),
192 theUseSymEETresholdFlag(useSymEETreshold),
194 theHcalThreshold(HcalThreshold),
195 theHBthreshold(HBthreshold),
196 theHESthreshold(HESthreshold),
197 theHEDthreshold(HEDthreshold),
198 theHOthreshold0(HOthreshold0),
199 theHOthresholdPlus1(HOthresholdPlus1),
200 theHOthresholdMinus1(HOthresholdMinus1),
201 theHOthresholdPlus2(HOthresholdPlus2),
202 theHOthresholdMinus2(HOthresholdMinus2),
203 theHF1threshold(HF1threshold),
204 theHF2threshold(HF2threshold),
206 theEBWeights(EBWeights),
208 theEEWeights(EEWeights),
210 theHBWeights(HBWeights),
212 theHESWeights(HESWeights),
214 theHEDWeights(HEDWeights),
216 theHOWeights(HOWeights),
218 theHF1Weights(HF1Weights),
220 theHF2Weights(HF2Weights),
221 theEBweight(EBweight),
222 theEEweight(EEweight),
223 theHBweight(HBweight),
224 theHESweight(HESweight),
225 theHEDweight(HEDweight),
226 theHOweight(HOweight),
227 theHF1weight(HF1weight),
228 theHF2weight(HF2weight),
229 theEcutTower(EcutTower),
230 theEBSumThreshold(EBSumThreshold),
231 theEESumThreshold(EESumThreshold),
234 theMomConstrMethod(momConstrMethod),
235 theMomHBDepth(momHBDepth),
236 theMomHEDepth(momHEDepth),
237 theMomEBDepth(momEBDepth),
238 theMomEEDepth(momEEDepth)
258 hbheItr != hbhe.
end(); ++hbheItr)
264 hoItr != ho.
end(); ++hoItr)
270 hfItr != hf.
end(); ++hfItr)
276 ecItr != ec.
end(); ++ecItr)
286 ctcItr != ctc.
end(); ++ctcItr) {
295 for(MetaTowerMap::const_iterator mapItr =
theTowerMap.begin();
300 if ( (mapItr->second).metaConstituents.size()<1)
continue;
314 ctcItr != ctc.
end(); ++ctcItr) {
317 double newE_em = ctcItr->emEnergy();
318 double newE_had = ctcItr->hadEnergy();
319 double newE_outer = ctcItr->outerEnergy();
325 if (ctcItr->ietaAbs()>=30) {
326 double E_short = 0.5 * newE_had;
327 double E_long = newE_em + 0.5 * newE_had;
332 newE_em = E_long - E_short;
333 newE_had = 2.0 * E_short;
339 for (
unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
340 DetId constId = ctcItr->constituent(iConst);
347 for (
unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
348 DetId constId = ctcItr->constituent(iConst);
356 for (
unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
357 DetId constId = ctcItr->constituent(iConst);
362 if (ctcItr->ietaAbs()>16) newE_outer *= weight;
370 double newE_hadTot = (
theHOIsUsed && twrId.
ietaAbs()<16)? newE_had+newE_outer : newE_had;
375 double f_em = 1.0/cosh(emPoint.
eta());
376 double f_had = 1.0/cosh(hadPoint.
eta());
380 if (ctcItr->ietaAbs()<30) {
385 double newE_tot = newE_em + newE_had;
392 CaloTower rescaledTower(twrId, newE_em, newE_had, newE_outer, -1, -1, towerP4, emPoint, hadPoint);
394 rescaledTower.
setEcalTime(
int(ctcItr->ecalTime()*100.0 + 0.5) );
395 rescaledTower.
setHcalTime(
int(ctcItr->hcalTime()*100.0 + 0.5) );
398 for (
unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
399 contains.push_back(ctcItr->constituent(iConst));
439 if (towerDetId.
null())
return;
452 tower29.numBadHcalCells += 1;
455 else if (0.5*energy >= threshold) {
459 tower29.numRecHcalCells += 1;
463 tower29.numProbHcalCells += 1;
467 double e28 = 0.5 *
e;
468 double e29 = 0.5 *
e;
470 tower28.
E_had += e28;
472 std::pair<DetId,double> mc(detId,e28);
475 tower29.E_had += e29;
477 tower29.metaConstituents.push_back(mc);
485 tower29.hadSumTimeTimesE += ( e29 * recHit->
time() );
486 tower29.hadSumEForTime += e29;
491 tower29.E_outer += e29;
500 if (towerDetId.
null())
return;
514 bool passEmThreshold =
false;
519 else passEmThreshold = (energy >=
threshold);
525 else passEmThreshold = (energy >=
threshold);
550 std::pair<DetId,double> mc(detId,e);
568 else if (energy >= threshold) {
584 std::pair<DetId,double> mc(detId,e);
598 else if (energy >= threshold) {
599 if (hcalDetId.
depth() == 1) {
623 std::pair<DetId,double> mc(detId,e);
636 else if (energy >= threshold) {
662 std::pair<DetId,double> mc(detId,e);
716 std::pair<DetId, double> mc(detId, 0);
728 E(0),E_em(0),E_had(0),E_outer(0), emSumTimeTimesE(0), hadSumTimeTimesE(0), emSumEForTime(0), hadSumEForTime(0),
729 numBadEcalCells(0), numRecEcalCells(0), numProbEcalCells(0), numBadHcalCells(0), numRecHcalCells(0), numProbHcalCells(0) { }
733 MetaTowerMap::iterator itr =
theTowerMap.find(detId);
739 theTowerMap.insert(std::pair<CaloTowerDetId, CaloTowersCreationAlgo::MetaTower>(detId, t));
752 double E_had=mt.
E_had;
758 if (
id.ietaAbs()<=17) {
781 if (
id.ietaAbs()<=29 && E_em<ecalThres) {
784 std::vector<std::pair<DetId,double> > metaContains_noecal;
786 for (std::vector<std::pair<DetId,double> >::iterator
i=metaContains.begin();
i!=metaContains.end(); ++
i)
787 if (
i->first.det()!=
DetId::Ecal) metaContains_noecal.push_back(*
i);
788 metaContains.swap(metaContains_noecal);
798 std::vector<std::pair<DetId,double> > metaContains_nohcal;
800 for (std::vector<std::pair<DetId,double> >::iterator
i=metaContains.begin();
i!=metaContains.end(); ++
i)
801 if (
i->first.det()!=
DetId::Hcal) metaContains_nohcal.push_back(*
i);
802 metaContains.swap(metaContains_nohcal);
806 double E_had_tot = (
theHOIsUsed &&
id.ietaAbs()<16)? E_had+E_outer : E_had;
822 double pf=1.0/cosh(p.
eta());
832 if (
id.ietaAbs()<=29) {
835 double emPf = 1.0/cosh(emPoint.
eta());
838 if ( (E_had + E_outer) >0) {
840 double hadPf = 1.0/cosh(hadPoint.
eta());
849 double pf=1.0/cosh(p.
eta());
859 if (
id.ietaAbs()<=29) {
863 double sumPf = 1.0/cosh(emPoint.
eta());
870 double pf=1.0/cosh(p.
eta());
881 CaloTower retval(
id, E_em, E_had, E_outer, -1, -1, towerP4, emPoint, hadPoint);
894 unsigned int numBadEcalChan = 0;
910 for (std::vector<DetId>::iterator ac_it=allConstituents.begin();
911 ac_it!=allConstituents.end(); ++ac_it) {
915 int thisEcalSevLvl = -999;
936 numRecHcalChan, numRecEcalChan,
937 numProbHcalChan, numProbEcalChan);
939 double maxCellE = -999.0;
942 for (std::vector<std::pair<DetId,double> >::iterator
i=metaContains.begin();
i!=metaContains.end(); ++
i) {
944 contains.push_back(
i->first);
946 if (maxCellE < i->
second) {
954 maxCellE =
i->second;
958 maxCellE =
i->second;
1035 else if(hcalDetId.
ieta() < 0) {
1050 if(hcalDetId.
depth() == 1) {
1068 std::cout <<
"BAD CELL det " << det << std::endl;
1117 if (fracDepth<0) fracDepth=0;
1118 else if (fracDepth>1) fracDepth=1;
1120 if (fracDepth>0.0) {
1123 0.25*( cv[4].
y() + cv[5].
y() + cv[6].
y() + cv[7].
y() ),
1124 0.25*( cv[4].
z() + cv[5].
z() + cv[6].
z() + cv[7].
z() ) );
1125 point += fracDepth * (backPoint-
point);
1135 if (fracDepth<0) fracDepth=0;
1136 else if (fracDepth>1) fracDepth=1;
1138 if (fracDepth>0.0) {
1141 0.25*( cv[4].
y() + cv[5].
y() + cv[6].
y() + cv[7].
y() ),
1142 0.25*( cv[4].
z() + cv[5].
z() + cv[6].
z() + cv[7].
z() ) );
1143 point += fracDepth * (backPoint-
point);
1151 float fracDepth,
double hadE) {
1163 std::vector<std::pair<DetId,double> >::iterator mc_it = metaContains.begin();
1164 for (; mc_it!=metaContains.end(); ++mc_it) {
1179 return GlobalPoint(hadX/nConst, hadY/nConst, hadZ/nConst);
1190 if (fracDepth < 0) fracDepth = 0;
1191 else if (fracDepth > 1) fracDepth = 1;
1195 int iEta = towerId.
ieta();
1196 int iPhi = towerId.
iphi();
1200 if (towerId.
ietaAbs() <= 14) {
1205 else if (towerId.
ietaAbs() == 15) {
1210 else if (towerId.
ietaAbs() == 16) {
1215 else if (towerId.
ietaAbs() == 17) {
1225 else if (towerId.
ietaAbs() <= 29) {
1229 if (iEta == 29) iEta = 28;
1230 if (iEta == -29) iEta = -28;
1233 else if (towerId.
ietaAbs() >= 30) {
1261 0.25 * (cv[4].
y() + cv[5].
y() + cv[6].
y() + cv[7].
y()),
1262 0.25 * (cv[4].
z() + cv[5].
z() + cv[6].
z() + cv[7].
z()));
1264 point += fracDepth * (backPoint -
point);
1271 float fracDepth,
double emE) {
1281 std::vector<std::pair<DetId,double> >::iterator mc_it = metaContains.begin();
1282 for (; mc_it!=metaContains.end(); ++mc_it) {
1285 double e = mc_it->second;
1301 float fracDepth,
double emE) {
1308 double sumWeights = 0;
1310 double crystalThresh = 0.015 * emE;
1312 std::vector<std::pair<DetId,double> >::iterator mc_it = metaContains.begin();
1313 for (; mc_it!=metaContains.end(); ++mc_it) {
1314 if (mc_it->second < 0)
continue;
1315 if (mc_it->first.det() ==
DetId::Ecal && mc_it->second > crystalThresh) sumEmE += mc_it->second;
1318 for (mc_it = metaContains.begin(); mc_it!=metaContains.end(); ++mc_it) {
1320 if (mc_it->first.det() !=
DetId::Ecal || mc_it->second < crystalThresh)
continue;
1324 weight = 4.2 +
log(mc_it->second/sumEmE);
1332 return GlobalPoint(emX/sumWeights, emY/sumWeights, emZ/sumWeights);
1341 const float timeUnit = 0.01;
1343 if (time> 300.0)
return 30000;
1344 if (time< -300.0)
return -30000;
1346 return int(time/timeUnit + 0.5);
1367 for (std::vector<DetId>::iterator it = allChanInStatusCont.begin(); it!=allChanInStatusCont.end(); ++it) {
1401 const uint32_t recHitFlag = hit->
flags();
1469 int severityLevel = 999;
std::vector< double > theHBGrid
const EcalChannelStatus * theEcalChStatus
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
bool theUseEtEETresholdFlag
const HcalChannelQuality * theHcalChStatus
size_t constituentsSize() const
unsigned int theEcalAcceptSeverityLevelForRejectedHit
bool contains(EventRange const &lh, EventID const &rh)
HcalSubdetector subdet() const
get the subdetector
int ietaAbs() const
get the absolute value of the tower ieta
const DetId & detid() const
DetId constituent(size_t i) const
void setHF1EScale(double scale)
void setCaloTowerStatus(unsigned int numBadHcalChan, unsigned int numBadEcalChan, unsigned int numRecHcalChan, unsigned int numRecEcalChan, unsigned int numProbHcalChan, unsigned int numProbEcalChan)
std::vector< double > theHESGrid
std::vector< double > theHEDGrid
GlobalPoint emCrystalShwrPos(DetId detId, float fracDepth)
math::PtEtaPhiMLorentzVector PolarLorentzVector
Lorentz vector.
Geom::Phi< T > phi() const
Global3DPoint GlobalPoint
std::vector< T >::const_iterator const_iterator
void push_back(T const &t)
void setHF2EScale(double scale)
void finish(CaloTowerCollection &destCollection)
GlobalPoint emShwrLogWeightPos(std::vector< std::pair< DetId, double > > &metaContains, float fracDepth, double totEmE)
edm::Handle< EcalRecHitCollection > theEeHandle
std::vector< double > theHOWeights
MetaTower & find(const CaloTowerDetId &id)
looks for a given tower in the internal cache. If it can't find it, it makes it.
std::vector< double > theEEGrid
void rescale(const CaloTower *ct)
CaloTower convert(const CaloTowerDetId &id, const MetaTower &mt)
unsigned int ecalChanStatusForCaloTower(const CaloRecHit *hit)
void setHottestCellE(double e)
void setHBEScale(double scale)
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
std::vector< DetId > constituentsOf(const CaloTowerDetId &id) const
Get the constituent detids for this tower id ( not yet implemented )
std::vector< double > theEEWeights
GlobalPoint hadShwrPos(std::vector< std::pair< DetId, double > > &metaContains, float fracDepth, double hadE)
std::vector< double > theHESWeights
U second(std::pair< T, U > const &p)
const CaloSubdetectorGeometry * theTowerGeometry
virtual double energy() const
energy
int depth() const
get the tower depth
void addConstituents(const std::vector< DetId > &ids)
unsigned int hcalChanStatusForCaloTower(const CaloRecHit *hit)
CaloTowerDetId towerOf(const DetId &id) const
Get the tower id for this det id (or null if not known)
std::vector< double > theHF2Grid
std::vector< double > theHEDWeights
unsigned int useRejectedRecoveredEcalHits
void rescaleTowers(const CaloTowerCollection &ctInput, CaloTowerCollection &ctResult)
static const int SubdetId
std::vector< DetId > getAllChannels() const
bool recoveredRecHit(const DetId &myid, const uint32_t &myflag) const
const CaloGeometry * theGeometry
int ieta() const
get the cell ieta
int iphi() const
get the tower iphi
std::vector< double > theHOGrid
bool dropChannel(const uint32_t &mystatus) const
const CaloTowerConstituentsMap * theTowerConstituentsMap
double theHOthresholdPlus1
std::vector< double > theHF1Weights
int firstHEDoublePhiRing() const
GlobalPoint hadSegmentShwrPos(DetId detId, float fracDepth)
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
bool theRecoveredHcalHitsAreUsed
bool theHOIsUsed
only affects energy and ET calculation. HO is still recorded in the tower
unsigned int theEcalAcceptSeverityLevel
void process(const HBHERecHitCollection &hbhe)
int ietaAbs() const
get the absolute value of the cell ieta
unsigned int useRejectedRecoveredHcalHits
const_iterator end() const
double theHOthresholdMinus2
GlobalPoint emShwrPos(std::vector< std::pair< DetId, double > > &metaContains, float fracDepth, double totEmE)
std::vector< double > theHBWeights
const HcalSeverityLevelComputer * theHcalSevLvlComputer
void setHESEScale(double scale)
void setHEDEScale(double scale)
bool theUseSymEBTresholdFlag
void setEEEScale(double scale)
CaloTowerDetId id() const
Log< T >::type log(const T &t)
int compactTime(float time)
GlobalPoint hadShwPosFromCells(DetId frontCell, DetId backCell, float fracDepth)
bool null() const
is this a null id ?
int zside() const
get the z-side of the tower (1/-1)
void getThresholdAndWeight(const DetId &detId, double &threshold, double &weight) const
helper method to look up the appropriate threshold & weight
int getSeverityLevel(const DetId &myid, const uint32_t &myflag, const uint32_t &mystatus) const
const HcalTopology * theHcalTopology
void setHOEScale(double scale)
std::map< CaloTowerDetId, int > hcalDropChMap
double theHOthresholdMinus1
const EcalSeverityLevelAlgo * theEcalSevLvlAlgo
const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id.
std::vector< double > theHF1Grid
perl if(1 lt scalar(@::datatypes))
void assignHit(const CaloRecHit *recHit)
adds a single hit to the tower
std::vector< double > theEBGrid
static int severityLevel(const DetId, const EcalRecHitCollection &, const EcalChannelStatus &, float recHitEtThreshold=5., SpikeId spId=kSwissCross, float spIdThreshold=0.95, float recHitEnergyThresholdForTiming=2., float recHitEnergyThresholdForEE=15, float spIdThresholdIEta85=0.999)
std::vector< double > theHF2Weights
int ieta() const
get the tower ieta
bool theUseEtEBTresholdFlag
uint32_t getValue() const
unsigned int theHcalAcceptSeverityLevelForRejectedHit
const Item * getValues(DetId fId) const
Detector det() const
get the detector field from this detid
void setGeometry(const CaloTowerConstituentsMap *cttopo, const HcalTopology *htopo, const CaloGeometry *geo)
const GlobalPoint & getPosition() const
void setEBEScale(double scale)
std::vector< double > theEBWeights
double theHOthresholdPlus2
bool theRecoveredEcalHitsAreUsed
unsigned int theHcalAcceptSeverityLevel
bool theUseSymEETresholdFlag
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
double outerEnergy() const
const_iterator begin() const
virtual const CornersVec & getCorners() const =0
edm::Handle< EcalRecHitCollection > theEbHandle
math::PtEtaPhiMLorentzVector PolarLorentzVector
Lorentz vector.