EcalTBAnalysisCoreTools/src/TBPositionCalc.cc

0001 #include "RecoTBCalo/EcalTBAnalysisCoreTools/interface/TBPositionCalc.h"
0002
0003 using namespace std;
0004
0005 TBPositionCalc::TBPositionCalc(const std::map<std::string, double> &providedParameters,
0006                                std::string const &fullMapName,
0007                                const CaloSubdetectorGeometry *passedGeometry) {
0008   // barrel geometry initialization

0009   if (passedGeometry == nullptr)
0010     throw(std::runtime_error("\n\n TBPositionCalc: wrong initialization.\n\n"));
0011   theGeometry_ = passedGeometry;
0012
0013   // parameters initialization

0014   param_LogWeighted_ = providedParameters.find("LogWeighted")->second;
0015   param_X0_ = providedParameters.find("X0")->second;
0016   param_T0_ = providedParameters.find("T0")->second;
0017   param_W0_ = providedParameters.find("W0")->second;
0018
0019   theTestMap = new EcalTBCrystalMap(fullMapName);
0020 }
0021
0022 TBPositionCalc::~TBPositionCalc() {
0023   if (theTestMap)
0024     delete theTestMap;
0025 }
0026
0027 CLHEP::Hep3Vector TBPositionCalc::CalculateTBPos(const std::vector<EBDetId> &upassedDetIds,
0028                                                  int myCrystal,
0029                                                  EcalRecHitCollection const *passedRecHitsMap) {
0030   std::vector<EBDetId> passedDetIds = upassedDetIds;
0031   // throw an error if the cluster was not initialized properly

0032   if (passedRecHitsMap == nullptr)
0033     throw(std::runtime_error("\n\n TBPositionCalc::CalculateTBPos called uninitialized.\n\n"));
0034
0035   // check DetIds are nonzero

0036   std::vector<EBDetId> validDetIds;
0037   std::vector<EBDetId>::const_iterator iter;
0038   for (iter = passedDetIds.begin(); iter != passedDetIds.end(); iter++) {
0039     if (((*iter) != DetId(0)) && (passedRecHitsMap->find(*iter) != passedRecHitsMap->end()))
0040       validDetIds.push_back(*iter);
0041   }
0042   passedDetIds.clear();
0043   passedDetIds = validDetIds;
0044
0045   // computing the position in the cms frame

0046   CLHEP::Hep3Vector cmsPos = CalculateCMSPos(passedDetIds, myCrystal, passedRecHitsMap);
0047
0048   // computing the rotation matrix (from CMS to TB)

0049   CLHEP::HepRotation *CMStoTB = new CLHEP::HepRotation();
0050   computeRotation(myCrystal, (*CMStoTB));
0051
0052   // moving to testbeam frame

0053   CLHEP::Hep3Vector tbPos = (*CMStoTB) * cmsPos;
0054   delete CMStoTB;
0055
0056   return tbPos;
0057 }
0058
0059 CLHEP::Hep3Vector TBPositionCalc::CalculateCMSPos(const std::vector<EBDetId> &passedDetIds,
0060                                                   int myCrystal,
0061                                                   EcalRecHitCollection const *passedRecHitsMap) {
0062   // Calculate the total energy

0063   double thisEne = 0;
0064   double eTot = 0;
0065   EBDetId myId;
0066   std::vector<EBDetId>::const_iterator myIt;
0067   for (myIt = passedDetIds.begin(); myIt != passedDetIds.end(); myIt++) {
0068     myId = (*myIt);
0069     EcalRecHitCollection::const_iterator itt = passedRecHitsMap->find(myId);
0070     thisEne = itt->energy();
0071     eTot += thisEne;
0072   }
0073
0074   // Calculate shower depth

0075   float depth = 0.;
0076   if (eTot <= 0.) {
0077     edm::LogError("NegativeClusterEnergy") << "cluster with negative energy: " << eTot << ", setting depth to 0.";
0078   } else {
0079     depth = param_X0_ * (param_T0_ + log(eTot));
0080   }
0081
0082   // Get position of the central crystal from shower depth

0083   EBDetId maxId_ = EBDetId(1, myCrystal, EBDetId::SMCRYSTALMODE);
0084   auto center_cell = theGeometry_->getGeometry(maxId_);
0085   GlobalPoint center_pos = center_cell->getPosition(depth);
0086
0087   // Loop over the hits collection

0088   double weight = 0;
0089   double total_weight = 0;
0090   double cluster_theta = 0;
0091   double cluster_phi = 0;
0092   std::vector<EBDetId>::const_iterator myIt2;
0093   for (myIt2 = passedDetIds.begin(); myIt2 != passedDetIds.end(); myIt2++) {
0094     // getting weights

0095     myId = (*myIt2);
0096     EcalRecHitCollection::const_iterator itj = passedRecHitsMap->find(myId);
0097     double ener = itj->energy();
0098
0099     if (param_LogWeighted_) {
0100       if (eTot <= 0.) {
0101         weight = 0.;
0102       } else {
0103         weight = std::max(0., param_W0_ + log(fabs(ener) / eTot));
0104       }
0105     } else {
0106       weight = ener / eTot;
0107     }
0108     total_weight += weight;
0109
0110     // weighted position of this detId

0111     auto jth_cell = theGeometry_->getGeometry(myId);
0112     GlobalPoint jth_pos = jth_cell->getPosition(depth);
0113     cluster_theta += weight * jth_pos.theta();
0114     cluster_phi += weight * jth_pos.phi();
0115   }
0116
0117   // normalizing

0118   cluster_theta /= total_weight;
0119   cluster_phi /= total_weight;
0120   if (cluster_phi > M_PI) {
0121     cluster_phi -= 2. * M_PI;
0122   }
0123   if (cluster_phi < -M_PI) {
0124     cluster_phi += 2. * M_PI;
0125   }
0126
0127   // position in the cms frame

0128   double radius =
0129       sqrt(center_pos.x() * center_pos.x() + center_pos.y() * center_pos.y() + center_pos.z() * center_pos.z());
0130   double xpos = radius * cos(cluster_phi) * sin(cluster_theta);
0131   double ypos = radius * sin(cluster_phi) * sin(cluster_theta);
0132   double zpos = radius * cos(cluster_theta);
0133
0134   return CLHEP::Hep3Vector(xpos, ypos, zpos);
0135 }
0136
0137 // rotation matrix to move from the CMS reference frame to the test beam one

0138 void TBPositionCalc::computeRotation(int MyCrystal, CLHEP::HepRotation &CMStoTB) {
0139   // taking eta/phi of the crystal

0140
0141   double myEta = 999.;
0142   double myPhi = 999.;
0143   double myTheta = 999.;
0144   theTestMap->findCrystalAngles(MyCrystal, myEta, myPhi);
0145   myTheta = 2.0 * atan(exp(-myEta));
0146
0147   // matrix

0148   CLHEP::HepRotation *fromCMStoTB = new CLHEP::HepRotation();
0149   double angle1 = 90. * CLHEP::deg - myPhi;
0150   CLHEP::HepRotationZ *r1 = new CLHEP::HepRotationZ(angle1);
0151   double angle2 = myTheta;
0152   CLHEP::HepRotationX *r2 = new CLHEP::HepRotationX(angle2);
0153   double angle3 = 90. * CLHEP::deg;
0154   CLHEP::HepRotationZ *r3 = new CLHEP::HepRotationZ(angle3);
0155   (*fromCMStoTB) *= (*r3);
0156   (*fromCMStoTB) *= (*r2);
0157   (*fromCMStoTB) *= (*r1);
0158
0159   CMStoTB = (*fromCMStoTB);
0160
0161   delete fromCMStoTB;
0162   delete r1;
0163   delete r2;
0164   delete r3;
0165 }