Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:14:25

0001 // -*- C++ -*-
0002 //
0003 // Package:    HTrackAssociator
0004 // Class:      HDetIdAssociator
0005 //
0006 /*
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Dmytro Kovalskyi
0015 // Modified for ECAL+HCAL by:  Michal Szleper
0016 //         Created:  Fri Apr 21 10:59:41 PDT 2006
0017 //
0018 //
0019 
0020 #include "Calibration/Tools/interface/DetIdAssociator.h"
0021 
0022 #include <memory>
0023 
0024 // surfaces is a vector of GlobalPoint representing outermost point on a cylinder
0025 std::vector<GlobalPoint> HDetIdAssociator::getTrajectory(const FreeTrajectoryState& ftsStart,
0026                                                          const std::vector<GlobalPoint>& surfaces) {
0027   check_setup();
0028   std::vector<GlobalPoint> trajectory;
0029   TrajectoryStateOnSurface tSOSDest;
0030   FreeTrajectoryState ftsCurrent = ftsStart;
0031 
0032   for (std::vector<GlobalPoint>::const_iterator surface_iter = surfaces.begin(); surface_iter != surfaces.end();
0033        surface_iter++) {
0034     // this stuff is some weird pointer, which destroy itself
0035     std::unique_ptr<Cylinder> cylinder =
0036         std::make_unique<Cylinder>(surface_iter->perp(), Surface::PositionType(0, 0, 0), Surface::RotationType());
0037     std::unique_ptr<Plane> forwardEndcap =
0038         std::make_unique<Plane>(Surface::PositionType(0, 0, surface_iter->z()), Surface::RotationType());
0039     std::unique_ptr<Plane> backwardEndcap =
0040         std::make_unique<Plane>(Surface::PositionType(0, 0, -surface_iter->z()), Surface::RotationType());
0041 
0042     LogTrace("StartingPoint") << "Propagate from "
0043                               << "\n"
0044                               << "\tx: " << ftsStart.position().x() << "\n"
0045                               << "\ty: " << ftsStart.position().y() << "\n"
0046                               << "\tz: " << ftsStart.position().z() << "\n"
0047                               << "\tmomentum eta: " << ftsStart.momentum().eta() << "\n"
0048                               << "\tmomentum phi: " << ftsStart.momentum().phi() << "\n"
0049                               << "\tmomentum: " << ftsStart.momentum().mag() << "\n";
0050 
0051     float tanTheta = ftsCurrent.momentum().perp() / ftsCurrent.momentum().z();
0052     float corner = surface_iter->perp() / surface_iter->z();
0053     /*
0054       std::cout<<"Propagate from "<< "\n"
0055         << "\tx: " << ftsCurrent.position().x()<< "\n"
0056         << "\ty: " << ftsCurrent.position().y()<< "\n"
0057         << "\tz: " << ftsCurrent.position().z()<< "\n"
0058         << "\tz: " << ftsCurrent.position().perp()<< "\n"
0059         << "\tz: " << tanTheta<<" "<< corner <<"\n"
0060         << "\tmomentum eta: " << ftsCurrent.momentum().eta()<< "\n"
0061         << "\tmomentum phi: " << ftsCurrent.momentum().phi()<< "\n"
0062         << "\tmomentum: " << ftsCurrent.momentum().mag()<<std::endl;
0063 */
0064     // First propage the track to the cylinder if |eta|<1, othewise to the encap
0065     // and correct depending on the result
0066     int ibar = 0;
0067     if (fabs(tanTheta) > corner) {
0068       tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
0069       //                  std::cout<<" Propagate to cylinder "<<std::endl;
0070     } else if (tanTheta > 0.) {
0071       tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap);
0072       ibar = 1;
0073     } else {
0074       tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap);
0075       ibar = -1;
0076     }
0077 
0078     //       std::cout<<" Trajectory valid? "<<tSOSDest.isValid()<<" First propagation in "<<ibar<<std::endl;
0079 
0080     if (!tSOSDest.isValid()) {
0081       // barrel
0082       if (ibar == 0) {
0083         if (tanTheta < 0)
0084           tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap);
0085         if (tanTheta >= 0)
0086           tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap);
0087       } else {
0088         tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
0089       }
0090     } else {
0091       // missed target
0092       if (abs(ibar) > 0) {
0093         if (tSOSDest.globalPosition().perp() > surface_iter->perp()) {
0094           tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
0095         }
0096       } else {
0097         if (tanTheta < 0)
0098           tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap);
0099         if (tanTheta >= 0)
0100           tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap);
0101       }
0102     }
0103 
0104     // If missed the target, propagate to again
0105     //      if ((!tSOSDest.isValid()) && point.perp() > surface_iter->perp())
0106     //  {tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);std::cout<<" Propagate again 1 "<<std::endl;}
0107     //         std::cout<<" Track is ok after repropagation to cylinder or not? "<<tSOSDest.isValid()<<std::endl;
0108     //      if ((!tSOSDest.isValid()) && ftsStart.momentum().eta()>0. && fabs(ftsStart.momentum().eta())>1.)
0109     //  {tSOSDest = ivProp_->propagate(ftsStart, *forwardEndcap);std::cout<<" Propagate again 2 "<<std::endl;}
0110     //       std::cout<<" Track is ok after repropagation forward or not? "<<tSOSDest.isValid()<<std::endl;
0111     //      if ((!tSOSDest.isValid()) && ftsStart.momentum().eta()<0.&&fabs(ftsStart.momentum().eta())>1.)
0112     //  {tSOSDest = ivProp_->propagate(ftsStart, *backwardEndcap);std::cout<<" Propagate again 3 "<<std::endl;}
0113     //       std::cout<<" Track is after repropagation backward ok or not? "<<tSOSDest.isValid()<<std::endl;
0114 
0115     if (!tSOSDest.isValid())
0116       return trajectory;
0117 
0118     //      std::cout<<" Propagate reach something"<<std::endl;
0119     LogTrace("SuccessfullPropagation") << "Great, I reached something."
0120                                        << "\n"
0121                                        << "\tx: " << tSOSDest.freeState()->position().x() << "\n"
0122                                        << "\ty: " << tSOSDest.freeState()->position().y() << "\n"
0123                                        << "\tz: " << tSOSDest.freeState()->position().z() << "\n"
0124                                        << "\teta: " << tSOSDest.freeState()->position().eta() << "\n"
0125                                        << "\tphi: " << tSOSDest.freeState()->position().phi() << "\n";
0126 
0127     //      std::cout<<" The position of trajectory "<<tSOSDest.freeState()->position().perp()<<" "<<tSOSDest.freeState()->position().z()<<std::endl;
0128 
0129     GlobalPoint point = tSOSDest.freeState()->position();
0130     point = tSOSDest.freeState()->position();
0131     ftsCurrent = *tSOSDest.freeState();
0132     trajectory.push_back(point);
0133   }
0134   return trajectory;
0135 }
0136 
0137 //------------------------------------------------------------------------------
0138 std::set<DetId> HDetIdAssociator::getDetIdsCloseToAPoint(const GlobalPoint& direction, const int idR) {
0139   std::set<DetId> set;
0140   check_setup();
0141   int nDets = 0;
0142   if (!theMap_)
0143     buildMap();
0144   LogTrace("MatchPoint") << "point (eta,phi): " << direction.eta() << "," << direction.phi() << "\n";
0145   int ieta = iEta(direction);
0146   int iphi = iPhi(direction);
0147 
0148   LogTrace("MatchPoint") << "(ieta,iphi): " << ieta << "," << iphi << "\n";
0149 
0150   if (ieta >= 0 && ieta < nEta_ && iphi >= 0 && iphi < nPhi_) {
0151     set = (*theMap_)[ieta][iphi];
0152     nDets++;
0153     if (idR > 0) {
0154       LogTrace("MatchPoint") << "Add neighbors (ieta,iphi): " << ieta << "," << iphi << "\n";
0155       //add neighbors
0156       int maxIEta = ieta + idR;
0157       int minIEta = ieta - idR;
0158       if (maxIEta >= nEta_)
0159         maxIEta = nEta_ - 1;
0160       if (minIEta < 0)
0161         minIEta = 0;
0162       int maxIPhi = iphi + idR;
0163       int minIPhi = iphi - idR;
0164       if (minIPhi < 0) {
0165         minIPhi += nPhi_;
0166         maxIPhi += nPhi_;
0167       }
0168       LogTrace("MatchPoint") << "\tieta (min,max): " << minIEta << "," << maxIEta << "\n";
0169       LogTrace("MatchPoint") << "\tiphi (min,max): " << minIPhi << "," << maxIPhi << "\n";
0170       for (int i = minIEta; i <= maxIEta; i++)
0171         for (int j = minIPhi; j <= maxIPhi; j++) {
0172           if (i == ieta && j == iphi)
0173             continue;  // already in the set
0174           set.insert((*theMap_)[i][j % nPhi_].begin(), (*theMap_)[i][j % nPhi_].end());
0175           nDets++;
0176         }
0177     }
0178   }
0179   //     if(set.size() > 0) {
0180   //   if (ieta+idR<55 && ieta-idR>14 && set.size() != (2*idR+1)*(2*idR+1)){
0181   //       std::cout<<" RRRA: "<<set.size()<<" DetIds in region "<<ieta<<" "<<iphi<<std::endl;
0182   //       for( std::set<DetId>::const_iterator itr=set.begin(); itr!=set.end(); itr++) {
0183   //         GlobalPoint point = getPosition(*itr);
0184   //         std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<" "<<iEta(point)<<" "<<iPhi(point)<<std::endl;
0185   //       }
0186   //     }
0187   //     else {
0188   //       std::cout <<" HDetIdAssociator::getDetIdsCloseToAPoint::There are strange days "<<std::endl;
0189   //     }
0190   return set;
0191 }
0192 
0193 //------------------------------------------------------------------------------
0194 int HDetIdAssociator::iEta(const GlobalPoint& point) {
0195   // unequal bin sizes for endcap, following HCAL geometry
0196   int iEta1 = int(point.eta() / etaBinSize_ + nEta_ / 2);
0197   if (point.eta() > 1.827 && point.eta() <= 1.830)
0198     return iEta1 - 1;
0199   else if (point.eta() > 1.914 && point.eta() <= 1.930)
0200     return iEta1 - 1;
0201   else if (point.eta() > 2.001 && point.eta() <= 2.043)
0202     return iEta1 - 1;
0203   else if (point.eta() > 2.088 && point.eta() <= 2.172)
0204     return iEta1 - 1;
0205   else if (point.eta() > 2.175 && point.eta() <= 2.262)
0206     return iEta1 - 1;
0207   else if (point.eta() > 2.262 && point.eta() <= 2.332)
0208     return iEta1 - 2;
0209   else if (point.eta() > 2.332 && point.eta() <= 2.349)
0210     return iEta1 - 1;
0211   else if (point.eta() > 2.349 && point.eta() <= 2.436)
0212     return iEta1 - 2;
0213   else if (point.eta() > 2.436 && point.eta() <= 2.500)
0214     return iEta1 - 3;
0215   else if (point.eta() > 2.500 && point.eta() <= 2.523)
0216     return iEta1 - 2;
0217   else if (point.eta() > 2.523 && point.eta() <= 2.610)
0218     return iEta1 - 3;
0219   else if (point.eta() > 2.610 && point.eta() <= 2.650)
0220     return iEta1 - 4;
0221   else if (point.eta() > 2.650 && point.eta() <= 2.697)
0222     return iEta1 - 3;
0223   else if (point.eta() > 2.697 && point.eta() <= 2.784)
0224     return iEta1 - 4;
0225   else if (point.eta() > 2.784 && point.eta() <= 2.868)
0226     return iEta1 - 5;
0227   else if (point.eta() > 2.868 && point.eta() <= 2.871)
0228     return iEta1 - 4;
0229   else if (point.eta() > 2.871 && point.eta() <= 2.958)
0230     return iEta1 - 5;
0231   else if (point.eta() > 2.958)
0232     return iEta1 - 6;
0233   else if (point.eta() < -1.827 && point.eta() >= -1.830)
0234     return iEta1 + 1;
0235   else if (point.eta() < -1.914 && point.eta() >= -1.930)
0236     return iEta1 + 1;
0237   else if (point.eta() < -2.001 && point.eta() >= -2.043)
0238     return iEta1 + 1;
0239   else if (point.eta() < -2.088 && point.eta() >= -2.172)
0240     return iEta1 + 1;
0241   else if (point.eta() < -2.175 && point.eta() >= -2.262)
0242     return iEta1 + 1;
0243   else if (point.eta() < -2.262 && point.eta() >= -2.332)
0244     return iEta1 + 2;
0245   else if (point.eta() < -2.332 && point.eta() >= -2.349)
0246     return iEta1 + 1;
0247   else if (point.eta() < -2.349 && point.eta() >= -2.436)
0248     return iEta1 + 2;
0249   else if (point.eta() < -2.436 && point.eta() >= -2.500)
0250     return iEta1 + 3;
0251   else if (point.eta() < -2.500 && point.eta() >= -2.523)
0252     return iEta1 + 2;
0253   else if (point.eta() < -2.523 && point.eta() >= -2.610)
0254     return iEta1 + 3;
0255   else if (point.eta() < -2.610 && point.eta() >= -2.650)
0256     return iEta1 + 4;
0257   else if (point.eta() < -2.650 && point.eta() >= -2.697)
0258     return iEta1 + 3;
0259   else if (point.eta() < -2.697 && point.eta() >= -2.784)
0260     return iEta1 + 4;
0261   else if (point.eta() < -2.784 && point.eta() >= -2.868)
0262     return iEta1 + 5;
0263   else if (point.eta() < -2.868 && point.eta() >= -2.871)
0264     return iEta1 + 4;
0265   else if (point.eta() < -2.871 && point.eta() >= -2.958)
0266     return iEta1 + 5;
0267   else if (point.eta() < -2.349)
0268     return iEta1 + 6;
0269   else
0270     return iEta1;
0271 }
0272 
0273 //------------------------------------------------------------------------------
0274 int HDetIdAssociator::iPhi(const GlobalPoint& point) {
0275   double pi = 4 * atan(1.);
0276   int iPhi1 = int((double(point.phi()) + pi) / (2 * pi) * nPhi_);
0277   return iPhi1;
0278 }
0279 
0280 //------------------------------------------------------------------------------
0281 void HDetIdAssociator::buildMap() {
0282   // modified version: take only detector central position
0283   check_setup();
0284   LogTrace("HDetIdAssociator") << "building map"
0285                                << "\n";
0286   if (theMap_)
0287     delete theMap_;
0288   theMap_ = new std::vector<std::vector<std::set<DetId> > >(nEta_, std::vector<std::set<DetId> >(nPhi_));
0289   int numberOfDetIdsOutsideEtaRange = 0;
0290   int numberOfDetIdsActive = 0;
0291   std::set<DetId> validIds = getASetOfValidDetIds();
0292   for (std::set<DetId>::const_iterator id_itr = validIds.begin(); id_itr != validIds.end(); id_itr++) {
0293     //      std::vector<GlobalPoint> points = getDetIdPoints(*id_itr);
0294     GlobalPoint point = getPosition(*id_itr);
0295     // reject fake DetIds (eta=0 - what are they anyway???)
0296     if (point.eta() == 0)
0297       continue;
0298 
0299     int ieta = iEta(point);
0300     int iphi = iPhi(point);
0301     int etaMax(-1);
0302     int etaMin(nEta_);
0303     int phiMax(-1);
0304     int phiMin(nPhi_);
0305     if (iphi >= nPhi_)
0306       iphi = iphi % nPhi_;
0307     assert(iphi >= 0);
0308     if (etaMin > ieta)
0309       etaMin = ieta;
0310     if (etaMax < ieta)
0311       etaMax = ieta;
0312     if (phiMin > iphi)
0313       phiMin = iphi;
0314     if (phiMax < iphi)
0315       phiMax = iphi;
0316     // for abs(eta)>1.8 one tower covers two phi segments
0317     if ((ieta > 54 || ieta < 15) && iphi % 2 == 0)
0318       phiMax++;
0319     if ((ieta > 54 || ieta < 15) && iphi % 2 == 1)
0320       phiMin--;
0321 
0322     if (etaMax < 0 || phiMax < 0 || etaMin >= nEta_ || phiMin >= nPhi_) {
0323       LogTrace("HDetIdAssociator") << "Out of range: DetId:" << id_itr->rawId() << "\n\teta (min,max): " << etaMin
0324                                    << "," << etaMax << "\n\tphi (min,max): " << phiMin << "," << phiMax
0325                                    << "\nTower id: " << id_itr->rawId() << "\n";
0326       numberOfDetIdsOutsideEtaRange++;
0327       continue;
0328     }
0329 
0330     if (phiMax - phiMin > phiMin + nPhi_ - phiMax) {
0331       phiMin += nPhi_;
0332       std::swap(phiMin, phiMax);
0333     }
0334     for (int ieta = etaMin; ieta <= etaMax; ieta++)
0335       for (int iphi = phiMin; iphi <= phiMax; iphi++)
0336         (*theMap_)[ieta][iphi % nPhi_].insert(*id_itr);
0337     numberOfDetIdsActive++;
0338   }
0339   LogTrace("HDetIdAssociator") << "Number of elements outside the allowed range ( |eta|>" << nEta_ / 2 * etaBinSize_
0340                                << "): " << numberOfDetIdsOutsideEtaRange << "\n";
0341   LogTrace("HDetIdAssociator") << "Number of active DetId's mapped: " << numberOfDetIdsActive << "\n";
0342 }
0343 
0344 //------------------------------------------------------------------------------
0345 std::set<DetId> HDetIdAssociator::getDetIdsInACone(const std::set<DetId>& inset,
0346                                                    const std::vector<GlobalPoint>& trajectory,
0347                                                    const double dR) {
0348   // modified version: if dR<0, returns 3x3 towers around the input one (Michal)
0349   check_setup();
0350   std::set<DetId> outset;
0351 
0352   if (dR >= 0) {
0353     for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
0354       for (std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end();
0355            point_iter++)
0356         if (nearElement(*point_iter, *id_iter, dR))
0357           outset.insert(*id_iter);
0358   } else {
0359     if (inset.size() != 1)
0360       return outset;
0361     std::set<DetId>::const_iterator id_inp = inset.begin();
0362     int ieta;
0363     int iphi;
0364     GlobalPoint point = getPosition(*id_inp);
0365     ieta = iEta(point);
0366     iphi = iPhi(point);
0367     for (int i = ieta - 1; i <= ieta + 1; i++) {
0368       for (int j = iphi - 1; j <= iphi + 1; j++) {
0369         //         if( i==ieta && j==iphi) continue;
0370         if (i < 0 || i >= nEta_)
0371           continue;
0372         int j2fill = j % nPhi_;
0373         if (j2fill < 0)
0374           j2fill += nPhi_;
0375         if ((*theMap_)[i][j2fill].empty())
0376           continue;
0377         outset.insert((*theMap_)[i][j2fill].begin(), (*theMap_)[i][j2fill].end());
0378       }
0379     }
0380   }
0381 
0382   //   if (outset.size() > 0) {
0383   //     std::cout<<" RRRA: DetIds in cone:"<<std::endl;
0384   //     for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
0385   //       GlobalPoint point = getPosition(*itr);
0386   //       std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
0387   //     }
0388   //   }
0389 
0390   return outset;
0391 }
0392 
0393 //------------------------------------------------------------------------------
0394 std::set<DetId> HDetIdAssociator::getCrossedDetIds(const std::set<DetId>& inset,
0395                                                    const std::vector<GlobalPoint>& trajectory) {
0396   check_setup();
0397   std::set<DetId> outset;
0398   for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
0399     for (std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end();
0400          point_iter++)
0401       if (insideElement(*point_iter, *id_iter))
0402         outset.insert(*id_iter);
0403   return outset;
0404 }
0405 
0406 //------------------------------------------------------------------------------
0407 std::set<DetId> HDetIdAssociator::getMaxEDetId(const std::set<DetId>& inset,
0408                                                edm::Handle<CaloTowerCollection> caloTowers) {
0409   // returns the most energetic tower in the NxN box (Michal)
0410   check_setup();
0411   std::set<DetId> outset;
0412   std::set<DetId>::const_iterator id_max = inset.begin();
0413   double Ehadmax = 0;
0414 
0415   for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
0416     DetId id(*id_iter);
0417     //     GlobalPoint point = getPosition(*id_iter);
0418     //     int ieta = iEta(point);
0419     //     int iphi = iPhi(point);
0420     CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
0421     if (tower != (*caloTowers).end() && tower->hadEnergy() > Ehadmax) {
0422       id_max = id_iter;
0423       Ehadmax = tower->hadEnergy();
0424     }
0425   }
0426 
0427   if (Ehadmax > 0)
0428     outset.insert(*id_max);
0429 
0430   //   if (outset.size() > 0) {
0431   //     std::cout<<" RRRA: Most energetic DetId:"<<std::endl;
0432   //     for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
0433   //       GlobalPoint point = getPosition(*itr);
0434   //       std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
0435   //     }
0436   //   }
0437 
0438   return outset;
0439 }
0440 
0441 //------------------------------------------------------------------------------
0442 std::set<DetId> HDetIdAssociator::getMaxEDetId(const std::set<DetId>& inset,
0443                                                edm::Handle<HBHERecHitCollection> recHits) {
0444   // returns the most energetic tower in the NxN box - from RecHits (Michal)
0445   check_setup();
0446   std::set<DetId> outset;
0447   std::set<DetId>::const_iterator id_max = inset.begin();
0448   double Ehadmax = 0;
0449 
0450   for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
0451     DetId id(*id_iter);
0452     //     GlobalPoint point = getPosition(*id_iter);
0453     //     int ieta = iEta(point);
0454     //     int iphi = iPhi(point);
0455     HBHERecHitCollection::const_iterator hit = (*recHits).find(id);
0456     if (hit != (*recHits).end() && hit->energy() > Ehadmax) {
0457       id_max = id_iter;
0458       Ehadmax = hit->energy();
0459     }
0460   }
0461 
0462   if (Ehadmax > 0)
0463     outset.insert(*id_max);
0464 
0465   //   if (outset.size() > 0) {
0466   //     std::cout<<" RRRA: Most energetic DetId:"<<std::endl;
0467   //     for( std::set<DetId>::const_iterator itr=outset.begin(); itr!=outset.end(); itr++) {
0468   //       GlobalPoint point = getPosition(*itr);
0469   //       std::cout << "DetId: " <<itr->rawId() <<" (eta,phi): " << point.eta() << "," << point.phi()<<std::endl;
0470   //     }
0471   //   }
0472 
0473   return outset;
0474 }