CMS 3D CMS Logo

List of all members | Public Member Functions | Private Attributes
TBPositionCalc Class Reference

#include <TBPositionCalc.h>

Public Member Functions

CLHEP::Hep3Vector CalculateCMSPos (const std::vector< EBDetId > &passedDetIds, int myCrystal, EcalRecHitCollection const *passedRecHitsMap)
 
CLHEP::Hep3Vector CalculateTBPos (const std::vector< EBDetId > &passedDetIds, int myCrystal, EcalRecHitCollection const *passedRecHitsMap)
 
void computeRotation (int myCrystal, CLHEP::HepRotation &CMStoTB)
 
 TBPositionCalc (const std::map< std::string, double > &providedParameters, const std::string &mapFile, const CaloSubdetectorGeometry *passedGeometry)
 
 TBPositionCalc ()
 
 ~TBPositionCalc ()
 

Private Attributes

bool param_LogWeighted_
 
Double32_t param_T0_
 
Double32_t param_W0_
 
Double32_t param_X0_
 
const CaloSubdetectorGeometrytheGeometry_
 
EcalTBCrystalMaptheTestMap
 

Detailed Description

Definition at line 25 of file TBPositionCalc.h.

Constructor & Destructor Documentation

◆ TBPositionCalc() [1/2]

TBPositionCalc::TBPositionCalc ( const std::map< std::string, double > &  providedParameters,
const std::string &  mapFile,
const CaloSubdetectorGeometry passedGeometry 
)

Definition at line 5 of file TBPositionCalc.cc.

7  {
8  // barrel geometry initialization
9  if (passedGeometry == nullptr)
10  throw(std::runtime_error("\n\n TBPositionCalc: wrong initialization.\n\n"));
11  theGeometry_ = passedGeometry;
12 
13  // parameters initialization
14  param_LogWeighted_ = providedParameters.find("LogWeighted")->second;
15  param_X0_ = providedParameters.find("X0")->second;
16  param_T0_ = providedParameters.find("T0")->second;
17  param_W0_ = providedParameters.find("W0")->second;
18 
19  theTestMap = new EcalTBCrystalMap(fullMapName);
20 }
Double32_t param_X0_
Double32_t param_W0_
const CaloSubdetectorGeometry * theGeometry_
EcalTBCrystalMap * theTestMap
Double32_t param_T0_

◆ TBPositionCalc() [2/2]

TBPositionCalc::TBPositionCalc ( )
inline

Definition at line 31 of file TBPositionCalc.h.

31 {}

◆ ~TBPositionCalc()

TBPositionCalc::~TBPositionCalc ( )

Definition at line 22 of file TBPositionCalc.cc.

22  {
23  if (theTestMap)
24  delete theTestMap;
25 }
EcalTBCrystalMap * theTestMap

Member Function Documentation

◆ CalculateCMSPos()

CLHEP::Hep3Vector TBPositionCalc::CalculateCMSPos ( const std::vector< EBDetId > &  passedDetIds,
int  myCrystal,
EcalRecHitCollection const *  passedRecHitsMap 
)

Definition at line 59 of file TBPositionCalc.cc.

References funct::cos(), hcalRecHitTable_cff::depth, edm::SortedCollection< T, SORT >::find(), CrabHelper::log, M_PI, WZElectronSkims53X_cff::max, PV3DBase< T, PVType, FrameType >::phi(), CosmicsPD_Skims::radius, funct::sin(), EBDetId::SMCRYSTALMODE, mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::theta(), mps_merge::weight, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

61  {
62  // Calculate the total energy
63  double thisEne = 0;
64  double eTot = 0;
65  EBDetId myId;
66  std::vector<EBDetId>::const_iterator myIt;
67  for (myIt = passedDetIds.begin(); myIt != passedDetIds.end(); myIt++) {
68  myId = (*myIt);
69  EcalRecHitCollection::const_iterator itt = passedRecHitsMap->find(myId);
70  thisEne = itt->energy();
71  eTot += thisEne;
72  }
73 
74  // Calculate shower depth
75  float depth = 0.;
76  if (eTot <= 0.) {
77  edm::LogError("NegativeClusterEnergy") << "cluster with negative energy: " << eTot << ", setting depth to 0.";
78  } else {
79  depth = param_X0_ * (param_T0_ + log(eTot));
80  }
81 
82  // Get position of the central crystal from shower depth
83  EBDetId maxId_ = EBDetId(1, myCrystal, EBDetId::SMCRYSTALMODE);
84  auto center_cell = theGeometry_->getGeometry(maxId_);
85  GlobalPoint center_pos = center_cell->getPosition(depth);
86 
87  // Loop over the hits collection
88  double weight = 0;
89  double total_weight = 0;
90  double cluster_theta = 0;
91  double cluster_phi = 0;
92  std::vector<EBDetId>::const_iterator myIt2;
93  for (myIt2 = passedDetIds.begin(); myIt2 != passedDetIds.end(); myIt2++) {
94  // getting weights
95  myId = (*myIt2);
96  EcalRecHitCollection::const_iterator itj = passedRecHitsMap->find(myId);
97  double ener = itj->energy();
98 
99  if (param_LogWeighted_) {
100  if (eTot <= 0.) {
101  weight = 0.;
102  } else {
103  weight = std::max(0., param_W0_ + log(fabs(ener) / eTot));
104  }
105  } else {
106  weight = ener / eTot;
107  }
108  total_weight += weight;
109 
110  // weighted position of this detId
111  auto jth_cell = theGeometry_->getGeometry(myId);
112  GlobalPoint jth_pos = jth_cell->getPosition(depth);
113  cluster_theta += weight * jth_pos.theta();
114  cluster_phi += weight * jth_pos.phi();
115  }
116 
117  // normalizing
118  cluster_theta /= total_weight;
119  cluster_phi /= total_weight;
120  if (cluster_phi > M_PI) {
121  cluster_phi -= 2. * M_PI;
122  }
123  if (cluster_phi < -M_PI) {
124  cluster_phi += 2. * M_PI;
125  }
126 
127  // position in the cms frame
128  double radius =
129  sqrt(center_pos.x() * center_pos.x() + center_pos.y() * center_pos.y() + center_pos.z() * center_pos.z());
130  double xpos = radius * cos(cluster_phi) * sin(cluster_theta);
131  double ypos = radius * sin(cluster_phi) * sin(cluster_theta);
132  double zpos = radius * cos(cluster_theta);
133 
134  return CLHEP::Hep3Vector(xpos, ypos, zpos);
135 }
virtual CellMayOwnPtr getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
T z() const
Definition: PV3DBase.h:61
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
Double32_t param_X0_
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::vector< EcalRecHit >::const_iterator const_iterator
Double32_t param_W0_
Definition: weight.py:1
Log< level::Error, false > LogError
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
const CaloSubdetectorGeometry * theGeometry_
T sqrt(T t)
Definition: SSEVec.h:23
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
#define M_PI
Double32_t param_T0_
static const int SMCRYSTALMODE
Definition: EBDetId.h:159
Geom::Theta< T > theta() const
Definition: PV3DBase.h:72

◆ CalculateTBPos()

CLHEP::Hep3Vector TBPositionCalc::CalculateTBPos ( const std::vector< EBDetId > &  passedDetIds,
int  myCrystal,
EcalRecHitCollection const *  passedRecHitsMap 
)

Definition at line 27 of file TBPositionCalc.cc.

References edm::SortedCollection< T, SORT >::end(), and edm::SortedCollection< T, SORT >::find().

29  {
30  std::vector<EBDetId> passedDetIds = upassedDetIds;
31  // throw an error if the cluster was not initialized properly
32  if (passedRecHitsMap == nullptr)
33  throw(std::runtime_error("\n\n TBPositionCalc::CalculateTBPos called uninitialized.\n\n"));
34 
35  // check DetIds are nonzero
36  std::vector<EBDetId> validDetIds;
37  std::vector<EBDetId>::const_iterator iter;
38  for (iter = passedDetIds.begin(); iter != passedDetIds.end(); iter++) {
39  if (((*iter) != DetId(0)) && (passedRecHitsMap->find(*iter) != passedRecHitsMap->end()))
40  validDetIds.push_back(*iter);
41  }
42  passedDetIds.clear();
43  passedDetIds = validDetIds;
44 
45  // computing the position in the cms frame
46  CLHEP::Hep3Vector cmsPos = CalculateCMSPos(passedDetIds, myCrystal, passedRecHitsMap);
47 
48  // computing the rotation matrix (from CMS to TB)
49  CLHEP::HepRotation *CMStoTB = new CLHEP::HepRotation();
50  computeRotation(myCrystal, (*CMStoTB));
51 
52  // moving to testbeam frame
53  CLHEP::Hep3Vector tbPos = (*CMStoTB) * cmsPos;
54  delete CMStoTB;
55 
56  return tbPos;
57 }
void computeRotation(int myCrystal, CLHEP::HepRotation &CMStoTB)
CLHEP::Hep3Vector CalculateCMSPos(const std::vector< EBDetId > &passedDetIds, int myCrystal, EcalRecHitCollection const *passedRecHitsMap)
Definition: DetId.h:17

◆ computeRotation()

void TBPositionCalc::computeRotation ( int  myCrystal,
CLHEP::HepRotation &  CMStoTB 
)

Definition at line 138 of file TBPositionCalc.cc.

References JetChargeProducer_cfi::exp, and diffTwoXMLs::r2.

138  {
139  // taking eta/phi of the crystal
140 
141  double myEta = 999.;
142  double myPhi = 999.;
143  double myTheta = 999.;
144  theTestMap->findCrystalAngles(MyCrystal, myEta, myPhi);
145  myTheta = 2.0 * atan(exp(-myEta));
146 
147  // matrix
148  CLHEP::HepRotation *fromCMStoTB = new CLHEP::HepRotation();
149  double angle1 = 90. * CLHEP::deg - myPhi;
150  CLHEP::HepRotationZ *r1 = new CLHEP::HepRotationZ(angle1);
151  double angle2 = myTheta;
152  CLHEP::HepRotationX *r2 = new CLHEP::HepRotationX(angle2);
153  double angle3 = 90. * CLHEP::deg;
154  CLHEP::HepRotationZ *r3 = new CLHEP::HepRotationZ(angle3);
155  (*fromCMStoTB) *= (*r3);
156  (*fromCMStoTB) *= (*r2);
157  (*fromCMStoTB) *= (*r1);
158 
159  CMStoTB = (*fromCMStoTB);
160 
161  delete fromCMStoTB;
162  delete r1;
163  delete r2;
164  delete r3;
165 }
EcalTBCrystalMap * theTestMap
void findCrystalAngles(const int thisCrysIndex, double &thisEta, double &thisPhi)

Member Data Documentation

◆ param_LogWeighted_

bool TBPositionCalc::param_LogWeighted_
private

Definition at line 46 of file TBPositionCalc.h.

◆ param_T0_

Double32_t TBPositionCalc::param_T0_
private

Definition at line 48 of file TBPositionCalc.h.

◆ param_W0_

Double32_t TBPositionCalc::param_W0_
private

Definition at line 49 of file TBPositionCalc.h.

◆ param_X0_

Double32_t TBPositionCalc::param_X0_
private

Definition at line 47 of file TBPositionCalc.h.

◆ theGeometry_

const CaloSubdetectorGeometry* TBPositionCalc::theGeometry_
private

Definition at line 53 of file TBPositionCalc.h.

◆ theTestMap

EcalTBCrystalMap* TBPositionCalc::theTestMap
private

Definition at line 51 of file TBPositionCalc.h.