8 if(passedGeometry ==
nullptr)
9 throw(std::runtime_error(
"\n\n TBPositionCalc: wrong initialization.\n\n"));
10 theGeometry_ = passedGeometry;
13 param_LogWeighted_ = providedParameters.find(
"LogWeighted")->second;
14 param_X0_ = providedParameters.find(
"X0")->second;
15 param_T0_ = providedParameters.find(
"T0")->second;
16 param_W0_ = providedParameters.find(
"W0")->second;
23 if (theTestMap)
delete theTestMap;
28 std::vector<EBDetId> passedDetIds = upassedDetIds;
30 if(passedRecHitsMap ==
nullptr)
31 throw(std::runtime_error(
"\n\n TBPositionCalc::CalculateTBPos called uninitialized.\n\n"));
34 std::vector<EBDetId> validDetIds;
35 std::vector<EBDetId>::const_iterator iter;
36 for (iter = passedDetIds.begin(); iter != passedDetIds.end(); iter++) {
37 if (((*iter) !=
DetId(0))
38 && (passedRecHitsMap->
find(*iter) != passedRecHitsMap->
end()))
39 validDetIds.push_back(*iter);
42 passedDetIds = validDetIds;
45 CLHEP::Hep3Vector cmsPos = CalculateCMSPos(passedDetIds, myCrystal, passedRecHitsMap);
48 CLHEP::HepRotation *CMStoTB =
new CLHEP::HepRotation();
49 computeRotation(myCrystal, (*CMStoTB));
52 CLHEP::Hep3Vector tbPos = (*CMStoTB)*cmsPos;
65 std::vector<EBDetId>::const_iterator myIt;
66 for (myIt = passedDetIds.begin(); myIt != passedDetIds.end(); myIt++) {
69 thisEne = itt->energy();
76 edm::LogError(
"NegativeClusterEnergy") <<
"cluster with negative energy: " << eTot <<
", setting depth to 0.";
78 depth = param_X0_ * (param_T0_ +
log(eTot));
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++) {
98 double ener = itj->energy();
100 if (param_LogWeighted_) {
101 if(eTot<=0.) { weight = 0.; }
102 else { weight =
std::max(0., param_W0_ +
log( fabs(ener)/eTot) ); }
111 cluster_theta += weight*jth_pos.
theta();
112 cluster_phi += weight*jth_pos.
phi();
116 cluster_theta /= total_weight;
117 cluster_phi /= total_weight;
118 if (cluster_phi >
M_PI) { cluster_phi -= 2.*
M_PI; }
119 if (cluster_phi < -
M_PI){ cluster_phi += 2.*
M_PI; }
122 double radius =
sqrt(center_pos.
x()*center_pos.
x() + center_pos.
y()*center_pos.
y() + center_pos.
z()*center_pos.
z());
123 double xpos = radius*
cos(cluster_phi)*
sin(cluster_theta);
124 double ypos = radius*
sin(cluster_phi)*
sin(cluster_theta);
125 double zpos = radius*
cos(cluster_theta);
127 return CLHEP::Hep3Vector(xpos, ypos, zpos);
137 double myTheta = 999.;
138 theTestMap->findCrystalAngles(MyCrystal, myEta, myPhi);
139 myTheta = 2.0*atan(
exp(-myEta));
142 CLHEP::HepRotation * fromCMStoTB =
new CLHEP::HepRotation();
143 double angle1 = 90.*deg - myPhi;
144 CLHEP::HepRotationZ *
r1 =
new CLHEP::HepRotationZ(angle1);
145 double angle2 = myTheta;
146 CLHEP::HepRotationX *
r2 =
new CLHEP::HepRotationX(angle2);
147 double angle3 = 90.*deg;
148 CLHEP::HepRotationZ * r3 =
new CLHEP::HepRotationZ(angle3);
149 (*fromCMStoTB) *= (*r3);
150 (*fromCMStoTB) *= (*r2);
151 (*fromCMStoTB) *= (*r1);
153 CMStoTB = (*fromCMStoTB);
void computeRotation(int myCrystal, CLHEP::HepRotation &CMStoTB)
Sin< T >::type sin(const T &t)
CLHEP::Hep3Vector CalculateTBPos(const std::vector< EBDetId > &passedDetIds, int myCrystal, EcalRecHitCollection const *passedRecHitsMap)
Geom::Phi< T > phi() const
std::vector< EcalRecHit >::const_iterator const_iterator
CaloGeometry const * getGeometry()
Geom::Theta< T > theta() const
Cos< T >::type cos(const T &t)
const_iterator end() const
CLHEP::Hep3Vector CalculateCMSPos(const std::vector< EBDetId > &passedDetIds, int myCrystal, EcalRecHitCollection const *passedRecHitsMap)
A base class to handle the particular shape of Ecal Xtals. Taken from ORCA Calorimetry Code...
iterator find(key_type k)
static const int SMCRYSTALMODE