File indexing completed on 2023-03-17 10:42:49
0001 #ifndef Calibration_HcalCalibALgos_CommonUsefulStuff_h
0002 #define Calibration_HcalCalibALgos_CommonUsefulStuff_h
0003
0004 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0005 #include "DataFormats/DetId/interface/DetId.h"
0006 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0007 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0008 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0009
0010 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0011 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0012
0013 #include "FWCore/Utilities/interface/Exception.h"
0014
0015 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0016 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0017 #include "CommonTools/UtilAlgos/interface/DeltaR.h"
0018
0019
0020
0021
0022
0023 struct MaxHit_struct {
0024 int iphihitm;
0025 int ietahitm;
0026 int depthhit;
0027 float hitenergy;
0028 float dr;
0029 GlobalPoint posMax;
0030 MaxHit_struct() : iphihitm(0), ietahitm(0), depthhit(0), hitenergy(-100), dr(0) {}
0031 };
0032
0033 inline double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint) {
0034
0035
0036
0037 const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(), caloPoint.z());
0038
0039 const GlobalVector caloIntersectUnitVector = caloIntersectVector.unit();
0040
0041 const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
0042
0043 const GlobalVector rechitUnitVector = rechitVector.unit();
0044 double dotprod = caloIntersectUnitVector.dot(rechitUnitVector);
0045 double rechitdist = caloIntersectVector.mag() / dotprod;
0046
0047 const GlobalVector effectiveRechitVector = rechitdist * rechitUnitVector;
0048 const GlobalPoint effectiveRechitPoint(
0049 effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());
0050
0051 GlobalVector distance_vector = effectiveRechitPoint - caloPoint;
0052
0053 if (dotprod > 0.) {
0054 return distance_vector.mag();
0055 } else {
0056 return 999999.;
0057 }
0058 }
0059
0060 inline double getDistInPlaneTrackDir(const GlobalPoint caloPoint,
0061 const GlobalVector caloVector,
0062 const GlobalPoint rechitPoint) {
0063
0064
0065
0066 const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(),
0067 caloPoint.z());
0068
0069 const GlobalVector caloUnitVector = caloVector.unit();
0070 const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
0071 const GlobalVector rechitUnitVector = rechitVector.unit();
0072 double dotprod_denominator = caloUnitVector.dot(rechitUnitVector);
0073 double dotprod_numerator = caloUnitVector.dot(caloIntersectVector);
0074 double rechitdist = dotprod_numerator / dotprod_denominator;
0075 const GlobalVector effectiveRechitVector = rechitdist * rechitUnitVector;
0076 const GlobalPoint effectiveRechitPoint(
0077 effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());
0078 GlobalVector distance_vector = effectiveRechitPoint - caloPoint;
0079 if (dotprod_denominator > 0. && dotprod_numerator > 0.) {
0080 return distance_vector.mag();
0081 } else {
0082 return 999999.;
0083 }
0084 }
0085
0086 inline double getDistInPlane(const GlobalVector trackDirection,
0087 const GlobalPoint caloPoint,
0088 const GlobalPoint rechitPoint,
0089 double coneHeight) {
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110 const GlobalVector heightVector = trackDirection * coneHeight;
0111
0112
0113
0114 const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(), caloPoint.z());
0115
0116
0117 const GlobalVector coneBaseVector = heightVector + caloIntersectVector;
0118
0119
0120 const GlobalPoint coneBasePoint(coneBaseVector.x(), coneBaseVector.y(), coneBaseVector.z());
0121
0122
0123 const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
0124 const GlobalVector rechitDirection = rechitVector.unit();
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142 double rechitdist = trackDirection.dot(coneBaseVector) / trackDirection.dot(rechitDirection);
0143
0144
0145
0146
0147 const GlobalVector effectiveRechitVector = rechitdist * rechitDirection;
0148 const GlobalPoint effectiveRechitPoint(
0149 effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());
0150
0151 GlobalVector distance_vector = effectiveRechitPoint - coneBasePoint;
0152 return distance_vector.mag();
0153 }
0154
0155
0156 inline double ecalEnergyInCone(const GlobalPoint center,
0157 double radius,
0158 const EcalRecHitCollection ecalCol,
0159 const CaloGeometry* geo) {
0160 double eECALcone = 0;
0161 std::vector<int> usedHitsEcal;
0162 usedHitsEcal.clear();
0163
0164 for (std::vector<EcalRecHit>::const_iterator ehit = ecalCol.begin(); ehit != ecalCol.end(); ehit++) {
0165
0166 bool hitIsUsed = false;
0167 int hitHashedIndex = -10000;
0168 if (ehit->id().subdetId() == EcalBarrel) {
0169 EBDetId did(ehit->id());
0170 hitHashedIndex = did.hashedIndex();
0171 }
0172 if (ehit->id().subdetId() == EcalEndcap) {
0173 EEDetId did(ehit->id());
0174 hitHashedIndex = did.hashedIndex();
0175 }
0176 for (uint32_t i = 0; i < usedHitsEcal.size(); i++) {
0177 if (usedHitsEcal[i] == hitHashedIndex)
0178 hitIsUsed = true;
0179 }
0180 if (hitIsUsed)
0181 continue;
0182 usedHitsEcal.push_back(hitHashedIndex);
0183
0184
0185 const GlobalPoint& pos = geo->getPosition((*ehit).detid());
0186
0187 if (getDistInPlaneSimple(center, pos) < radius) {
0188 eECALcone += ehit->energy();
0189 }
0190 }
0191 return eECALcone;
0192 }
0193
0194
0195 inline double ecalEnergyInCone(const GlobalVector trackMom,
0196 const GlobalPoint center,
0197 double radius,
0198 const EcalRecHitCollection ecalCol,
0199 const CaloGeometry* geo) {
0200 double eECALcone = 0;
0201 std::vector<int> usedHitsEcal;
0202 usedHitsEcal.clear();
0203 for (std::vector<EcalRecHit>::const_iterator ehit = ecalCol.begin(); ehit != ecalCol.end(); ehit++) {
0204
0205 bool hitIsUsed = false;
0206 int hitHashedIndex = -10000;
0207 if (ehit->id().subdetId() == EcalBarrel) {
0208 EBDetId did(ehit->id());
0209 hitHashedIndex = did.hashedIndex();
0210 }
0211 if (ehit->id().subdetId() == EcalEndcap) {
0212 EEDetId did(ehit->id());
0213 hitHashedIndex = did.hashedIndex();
0214 }
0215 for (uint32_t i = 0; i < usedHitsEcal.size(); i++) {
0216 if (usedHitsEcal[i] == hitHashedIndex)
0217 hitIsUsed = true;
0218 }
0219 if (hitIsUsed)
0220 continue;
0221 usedHitsEcal.push_back(hitHashedIndex);
0222
0223
0224 const GlobalPoint& pos = geo->getPosition((*ehit).detid());
0225
0226 if (getDistInPlaneTrackDir(center, trackMom, pos) < radius) {
0227 eECALcone += ehit->energy();
0228 }
0229 }
0230 return eECALcone;
0231 }
0232
0233 #endif