Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:04 

0001 #include <iostream>
0002 #include <fstream>
0003 #include <algorithm>
0004 #include <numeric>
0005 
0006 #include "TString.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "Calibration/HcalCalibAlgos/interface/hcalCalibUtils.h"
0009 
0010 //#include "Calibration/HcalCalibAlgos/plugins/CommonUsefulStuff.h"
0011 #include "Calibration/HcalCalibAlgos/interface/CommonUsefulStuff.h"
0012 #include "DataFormats/DetId/interface/DetId.h"
0013 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0014 
0015 void sumDepths(std::vector<TCell>& selectCells) {
0016   // Assignes teh sum of the energy in cells with the same iEta, iPhi to the cell with depth=1.
0017   // All cells with depth>1 are removed form the container. If
0018   // the cell at depth=1 is not present: create it and follow the procedure.
0019 
0020   if (selectCells.empty())
0021     return;
0022 
0023   std::vector<TCell> selectCellsDepth1;
0024   std::vector<TCell> selectCellsHighDepth;
0025 
0026   //
0027   // NB: Here we add depth 3 for iEta==16 in *HE* to the value in the barrel
0028   // this approach is reflected in several of the following loops: make sure
0029   // to check when making changes.
0030   //
0031   // In some documents it is described as having depth 1, the mapping in CMSSW uses depth 3.
0032 
0033   for (std::vector<TCell>::iterator i_it = selectCells.begin(); i_it != selectCells.end(); ++i_it) {
0034     if (HcalDetId(i_it->id()).depth() == 1) {
0035       selectCellsDepth1.push_back(*i_it);
0036     } else {
0037       selectCellsHighDepth.push_back(*i_it);
0038     }
0039   }
0040 
0041   // case where depth 1 has zero energy, but higher depths with same (iEta, iPhi) have energy.
0042   // For iEta<15 there is one depth -> selectCellsHighDepth is empty and we do not get in the loop.
0043   for (std::vector<TCell>::iterator i_it2 = selectCellsHighDepth.begin(); i_it2 != selectCellsHighDepth.end();
0044        ++i_it2) {
0045     // protect against corrupt data
0046     if (HcalDetId(i_it2->id()).ietaAbs() < 15 && HcalDetId(i_it2->id()).depth() > 1) {
0047       edm::LogWarning("HcalCalib") << "ERROR!!! there are no HB cells with depth>1 for iEta<15!\n"
0048                                    << "Check the input data...\nHCalDetId: " << HcalDetId(i_it2->id());
0049       return;
0050     }
0051 
0052     bool foundDepthOne = false;
0053     for (std::vector<TCell>::iterator i_it = selectCellsDepth1.begin(); i_it != selectCellsDepth1.end(); ++i_it) {
0054       if (HcalDetId(i_it->id()).ieta() == HcalDetId(i_it2->id()).ieta() &&
0055           HcalDetId(i_it->id()).iphi() == HcalDetId(i_it2->id()).iphi())
0056         foundDepthOne = true;
0057       continue;
0058     }
0059     if (!foundDepthOne) {  // create entry for depth 1 with 0 energy
0060 
0061       UInt_t newId;
0062       if (abs(HcalDetId(i_it2->id()).ieta()) == 16)
0063         newId = HcalDetId(HcalBarrel, HcalDetId(i_it2->id()).ieta(), HcalDetId(i_it2->id()).iphi(), 1);
0064       else
0065         newId =
0066             HcalDetId(HcalDetId(i_it2->id()).subdet(), HcalDetId(i_it2->id()).ieta(), HcalDetId(i_it2->id()).iphi(), 1);
0067 
0068       selectCellsDepth1.push_back(TCell(newId, 0.0));
0069     }
0070   }
0071 
0072   for (std::vector<TCell>::iterator i_it = selectCellsDepth1.begin(); i_it != selectCellsDepth1.end(); ++i_it) {
0073     for (std::vector<TCell>::iterator i_it2 = selectCellsHighDepth.begin(); i_it2 != selectCellsHighDepth.end();
0074          ++i_it2) {
0075       if (HcalDetId(i_it->id()).ieta() == HcalDetId(i_it2->id()).ieta() &&
0076           HcalDetId(i_it->id()).iphi() == HcalDetId(i_it2->id()).iphi()) {
0077         i_it->SetE(i_it->e() + i_it2->e());
0078         i_it2->SetE(0.0);  // paranoid, aren't we...
0079       }
0080     }
0081   }
0082 
0083   // replace the original vectors with the new ones
0084   selectCells = selectCellsDepth1;
0085 
0086   return;
0087 }
0088 
0089 void combinePhi(std::vector<TCell>& selectCells) {
0090   // Map: NxN -> N cluster
0091   // Comine the targetE of cells with the same iEta
0092 
0093   if (selectCells.empty())
0094     return;
0095 
0096   // new container for the TCells
0097   // dummy cell id created with iEta; iPhi=1; depth
0098   // if combinePhi() is run after combining depths, depth=1
0099   std::vector<TCell> combinedCells;
0100 
0101   std::map<UInt_t, std::vector<Float_t> > etaSliceE;  // keyed by id of cell with iEta and **iPhi=1**
0102 
0103   // map the cells to the eta ring
0104   std::vector<TCell>::iterator i_it = selectCells.begin();
0105   for (; i_it != selectCells.end(); ++i_it) {
0106     DetId id = HcalDetId(i_it->id());
0107     UInt_t thisKey = HcalDetId(HcalDetId(id).subdet(), HcalDetId(id).ieta(), 1, HcalDetId(id).depth());
0108     etaSliceE[thisKey].push_back(i_it->e());
0109   }
0110 
0111   std::map<UInt_t, std::vector<Float_t> >::iterator m_it = etaSliceE.begin();
0112   for (; m_it != etaSliceE.end(); ++m_it) {
0113     combinedCells.push_back(TCell(m_it->first, accumulate(m_it->second.begin(), m_it->second.end(), 0.0)));
0114   }
0115 
0116   // replace the original TCell vector with the new one
0117   selectCells = combinedCells;
0118 }
0119 
0120 void combinePhi(std::vector<TCell>& selectCells, std::vector<TCell>& combinedCells) {
0121   // Map: NxN -> N cluster
0122   // Comine the targetE of cells with the same iEta
0123 
0124   if (selectCells.empty())
0125     return;
0126 
0127   std::map<UInt_t, std::vector<Float_t> > etaSliceE;  // keyed by id of cell with iEta and **iPhi=1**
0128 
0129   // map the cells to the eta ring
0130   std::vector<TCell>::iterator i_it = selectCells.begin();
0131   for (; i_it != selectCells.end(); ++i_it) {
0132     DetId id = HcalDetId(i_it->id());
0133     UInt_t thisKey = HcalDetId(HcalDetId(id).subdet(), HcalDetId(id).ieta(), 1, HcalDetId(id).depth());
0134     etaSliceE[thisKey].push_back(i_it->e());
0135   }
0136 
0137   std::map<UInt_t, std::vector<Float_t> >::iterator m_it = etaSliceE.begin();
0138   for (; m_it != etaSliceE.end(); ++m_it) {
0139     combinedCells.push_back(TCell(m_it->first, accumulate(m_it->second.begin(), m_it->second.end(), 0.0)));
0140   }
0141 }
0142 
0143 void getIEtaIPhiForHighestE(std::vector<TCell>& selectCells, Int_t& iEtaMostE, UInt_t& iPhiMostE) {
0144   std::vector<TCell> summedDepthsCells = selectCells;
0145 
0146   sumDepths(summedDepthsCells);
0147   std::vector<TCell>::iterator highCell = summedDepthsCells.begin();
0148 
0149   // sum depths locally to get highest energy tower
0150 
0151   Float_t highE = -999;
0152 
0153   for (std::vector<TCell>::iterator it = summedDepthsCells.begin(); it != summedDepthsCells.end(); ++it) {
0154     if (highE < it->e()) {
0155       highCell = it;
0156       highE = it->e();
0157     }
0158   }
0159 
0160   iEtaMostE = HcalDetId(highCell->id()).ieta();
0161   iPhiMostE = HcalDetId(highCell->id()).iphi();
0162 
0163   return;
0164 }
0165 
0166 //
0167 // Remove RecHits outside the 3x3 cluster and replace the  vector that will
0168 // be used in the minimization. Acts on "event" level.
0169 // This can not be done for iEta>20 due to segmentation => in principle the result should be restricted
0170 // to iEta<20. Attempted to minimize affect at the boundary without a sharp jump.
0171 
0172 void filterCells3x3(std::vector<TCell>& selectCells, Int_t iEtaMaxE, UInt_t iPhiMaxE) {
0173   std::vector<TCell> filteredCells;
0174 
0175   Int_t dEta, dPhi;
0176 
0177   for (std::vector<TCell>::iterator it = selectCells.begin(); it != selectCells.end(); ++it) {
0178     Bool_t passDEta = false;
0179     Bool_t passDPhi = false;
0180 
0181     dEta = HcalDetId(it->id()).ieta() - iEtaMaxE;
0182     dPhi = HcalDetId(it->id()).iphi() - iPhiMaxE;
0183 
0184     if (dPhi > 36)
0185       dPhi -= 72;
0186     if (dPhi < -36)
0187       dPhi += 72;
0188 
0189     if (abs(dEta) <= 1 || (iEtaMaxE * HcalDetId(it->id()).ieta() == -1))
0190       passDEta = true;
0191 
0192     if (abs(iEtaMaxE) <= 20) {
0193       if (abs(HcalDetId(it->id()).ieta()) <= 20) {
0194         if (abs(dPhi) <= 1)
0195           passDPhi = true;
0196       } else {
0197         // iPhi is labelled by odd numbers
0198         if (iPhiMaxE % 2 == 0) {
0199           if (abs(dPhi) <= 1)
0200             passDPhi = true;
0201         } else {
0202           if (dPhi == -2 || dPhi == 0)
0203             passDPhi = true;
0204         }
0205       }
0206 
0207     }  // if hottest cell with iEta<=20
0208 
0209     else {
0210       if (abs(HcalDetId(it->id()).ieta()) <= 20) {
0211         if (abs(dPhi) <= 1 || dPhi == 2)
0212           passDPhi = true;
0213       } else {
0214         if (abs(dPhi) <= 2)
0215           passDPhi = true;
0216       }
0217     }  // if hottest cell with iEta>20
0218 
0219     if (passDEta && passDPhi)
0220       filteredCells.push_back(*it);
0221   }
0222 
0223   selectCells = filteredCells;
0224 
0225   return;
0226 }
0227 
0228 //
0229 // Remove RecHits outside the 5x5 cluster and replace the  vector that will
0230 // be used in the minimization. Acts on "event" level
0231 // In principle the ntuple should be produced with 5x5 already precelected
0232 //
0233 // Size for iEta>20 is 3x3, but the segmentation changes by x2 in phi.
0234 // There is some bias in the selection of towers near the boundary
0235 
0236 void filterCells5x5(std::vector<TCell>& selectCells, Int_t iEtaMaxE, UInt_t iPhiMaxE) {
0237   std::vector<TCell> filteredCells;
0238 
0239   Int_t dEta, dPhi;
0240 
0241   for (std::vector<TCell>::iterator it = selectCells.begin(); it != selectCells.end(); ++it) {
0242     dEta = HcalDetId(it->id()).ieta() - iEtaMaxE;
0243     dPhi = HcalDetId(it->id()).iphi() - iPhiMaxE;
0244 
0245     if (dPhi > 36)
0246       dPhi -= 72;
0247     if (dPhi < -36)
0248       dPhi += 72;
0249 
0250     bool passDPhi = (abs(dPhi) < 3);
0251 
0252     bool passDEta = (abs(dEta) < 3 || (iEtaMaxE * HcalDetId(it->id()).ieta() == -2));
0253     // includes  +/- eta boundary
0254 
0255     if (passDPhi && passDEta)
0256       filteredCells.push_back(*it);
0257   }
0258 
0259   selectCells = filteredCells;
0260 
0261   return;
0262 }
0263 
0264 // this is for the problematic layer near the HB/HE boundary
0265 // sum depths 1,2 in towers 15,16
0266 
0267 void sumSmallDepths(std::vector<TCell>& selectCells) {
0268   if (selectCells.empty())
0269     return;
0270 
0271   std::vector<TCell> newCells;          // holds unaffected cells to which the modified ones are added
0272   std::vector<TCell> manipulatedCells;  // the ones that are combined
0273 
0274   for (std::vector<TCell>::iterator i_it = selectCells.begin(); i_it != selectCells.end(); ++i_it) {
0275     if ((HcalDetId(i_it->id()).ietaAbs() == 15 && HcalDetId(i_it->id()).depth() <= 2) ||
0276         (HcalDetId(i_it->id()).ietaAbs() == 16 && HcalDetId(i_it->id()).depth() <= 2)) {
0277       manipulatedCells.push_back(*i_it);
0278     } else {
0279       newCells.push_back(*i_it);
0280     }
0281   }
0282 
0283   // if the list is empty there is nothing to manipulate
0284   // leave the original vector unchanged
0285 
0286   if (manipulatedCells.empty()) {
0287     newCells.clear();
0288     return;
0289   }
0290 
0291   // See what cells are needed to hold the combined information:
0292   // Make holders for depth=1 for each (iEta,iPhi)
0293   // if a cell with these values is present in "manupulatedCells"
0294   std::vector<UInt_t> dummyIds;     // to keep track of kreated cells
0295   std::vector<TCell> createdCells;  // cells that need to be added or they exists;
0296 
0297   for (std::vector<TCell>::iterator i_it = manipulatedCells.begin(); i_it != manipulatedCells.end(); ++i_it) {
0298     UInt_t dummyId =
0299         HcalDetId(HcalDetId(i_it->id()).subdet(), HcalDetId(i_it->id()).ieta(), HcalDetId(i_it->id()).iphi(), 1);
0300     if (find(dummyIds.begin(), dummyIds.end(), dummyId) == dummyIds.end()) {
0301       dummyIds.push_back(dummyId);
0302       createdCells.push_back(TCell(dummyId, 0.0));
0303     }
0304   }
0305 
0306   for (std::vector<TCell>::iterator i_it = createdCells.begin(); i_it != createdCells.end(); ++i_it) {
0307     for (std::vector<TCell>::iterator i_it2 = manipulatedCells.begin(); i_it2 != manipulatedCells.end(); ++i_it2) {
0308       if (HcalDetId(i_it->id()).ieta() == HcalDetId(i_it2->id()).ieta() &&
0309           HcalDetId(i_it->id()).iphi() == HcalDetId(i_it2->id()).iphi() && HcalDetId(i_it2->id()).depth() <= 2) {
0310         i_it->SetE(i_it->e() + i_it2->e());
0311       }
0312     }
0313   }
0314 
0315   for (std::vector<TCell>::iterator i_it = createdCells.begin(); i_it != createdCells.end(); ++i_it) {
0316     newCells.push_back(*i_it);
0317   }
0318 
0319   // replace the original vectors with the new ones
0320   selectCells = newCells;
0321 
0322   return;
0323 }
0324 
0325 void filterCellsInCone(std::vector<TCell>& selectCells,
0326                        const GlobalPoint hitPositionHcal,
0327                        Float_t maxConeDist,
0328                        const CaloGeometry* theCaloGeometry) {
0329   std::vector<TCell> filteredCells;
0330 
0331   for (std::vector<TCell>::iterator it = selectCells.begin(); it != selectCells.end(); ++it) {
0332     GlobalPoint recHitPoint;
0333     DetId id = it->id();
0334     if (id.det() == DetId::Hcal) {
0335       recHitPoint = (static_cast<const HcalGeometry*>(theCaloGeometry->getSubdetectorGeometry(id)))->getPosition(id);
0336     } else {
0337       recHitPoint = GlobalPoint(theCaloGeometry->getPosition(id));
0338     }
0339 
0340     if (getDistInPlaneSimple(hitPositionHcal, recHitPoint) <= maxConeDist)
0341       filteredCells.push_back(*it);
0342   }
0343 
0344   selectCells = filteredCells;
0345 
0346   return;
0347 }
0348 
0349 // From Jim H. => keep till the code is included centrally
0350 /*
0351 double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint) {
0352   
0353   // Simplified version of getDistInPlane
0354   // Assume track direction is origin -> point of hcal intersection
0355   
0356   const GlobalVector caloIntersectVector(caloPoint.x(), 
0357                      caloPoint.y(), 
0358                      caloPoint.z());
0359 
0360   const GlobalVector caloIntersectUnitVector = caloIntersectVector.unit();
0361   
0362   const GlobalVector rechitVector(rechitPoint.x(),
0363                   rechitPoint.y(),
0364                   rechitPoint.z());
0365 
0366   const GlobalVector rechitUnitVector = rechitVector.unit();
0367 
0368   double dotprod = caloIntersectUnitVector.dot(rechitUnitVector);
0369   double rechitdist = caloIntersectVector.mag()/dotprod;
0370   
0371   
0372   const GlobalVector effectiveRechitVector = rechitdist*rechitUnitVector;
0373   const GlobalPoint effectiveRechitPoint(effectiveRechitVector.x(),
0374                      effectiveRechitVector.y(),
0375                      effectiveRechitVector.z());
0376   
0377   
0378   GlobalVector distance_vector = effectiveRechitPoint-caloPoint;
0379   
0380   if (dotprod > 0.)
0381     {
0382       return distance_vector.mag();
0383     }
0384   else
0385     {
0386       return 999999.;
0387     
0388     }
0389 
0390     
0391 }
0392 */
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!