CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/Geometry/CaloTopology/src/HcalTopology.cc

Go to the documentation of this file.
00001 #include "Geometry/CaloTopology/interface/HcalTopology.h"
00002 #include <cmath>
00003 #include <iostream>
00004 #include <cassert>
00005 #include <algorithm>
00006 #include "FWCore/Utilities/interface/Exception.h"
00007 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00008 #include "DataFormats/HcalDetId/interface/HcalTrigTowerDetId.h"
00009 #include "DataFormats/HcalDetId/interface/HcalCalibDetId.h"
00010 
00011 static const int IPHI_MAX=72;
00012 
00013 HcalTopology::HcalTopology(HcalTopologyMode::Mode mode, int maxDepthHB, int maxDepthHE, HcalTopologyMode::TriggerMode tmode) :
00014   excludeHB_(false),
00015   excludeHE_(false),
00016   excludeHO_(false),
00017   excludeHF_(false),
00018   mode_(mode),
00019   triggerMode_(tmode),
00020   firstHBRing_(1),
00021   lastHBRing_(16),
00022   firstHERing_(16),
00023   lastHERing_(29),
00024   firstHFRing_(29),
00025   lastHFRing_(41),
00026   firstHORing_(1),
00027   lastHORing_(15),
00028   firstHEDoublePhiRing_((mode==HcalTopologyMode::H2 || mode==HcalTopologyMode::H2HE)?(22):(21)),
00029   firstHFQuadPhiRing_(40),
00030   firstHETripleDepthRing_((mode==HcalTopologyMode::H2 || mode==HcalTopologyMode::H2HE)?(24):(27)),
00031   singlePhiBins_(72),
00032   doublePhiBins_(36),
00033   maxDepthHB_(maxDepthHB),
00034   maxDepthHE_(maxDepthHE),
00035   HBSize_(kHBSizePreLS1),
00036   HESize_(kHESizePreLS1),
00037   HOSize_(kHOSizePreLS1),
00038   HFSize_(kHFSizePreLS1),
00039   HTSize_(kHTSizePreLS1),
00040   numberOfShapes_(( mode==HcalTopologyMode::SLHC ) ? 500 : 87 ) {
00041 
00042   if (mode_==HcalTopologyMode::LHC) {
00043     topoVersion_=0; //DL
00044     HBSize_= kHBSizePreLS1; // qie-per-fiber * fiber/rm * rm/rbx * rbx/barrel * barrel/hcal
00045     HESize_= kHESizePreLS1; // qie-per-fiber * fiber/rm * rm/rbx * rbx/endcap * endcap/hcal
00046     HOSize_= kHOSizePreLS1; // ieta * iphi * 2
00047     HFSize_= kHFSizePreLS1; // phi * eta * depth * pm 
00048     HTSize_= kHTSizePreLS1; 
00049   } else if (mode_==HcalTopologyMode::SLHC) { // need to know more eventually
00050     HBSize_= maxDepthHB*16*72*2;
00051     HESize_= maxDepthHE*(29-16+1)*72*2;
00052     HOSize_= 15*72*2; // ieta * iphi * 2
00053     HFSize_= 72*13*2*2; // phi * eta * depth * pm 
00054     HTSize_= kHTSizePreLS1;  //no clue!
00055 
00056     topoVersion_=10;
00057   }
00058     
00059 }
00060 
00061 bool HcalTopology::valid(const DetId& id) const {
00062   assert(id.det()==DetId::Hcal);
00063   return validHcal(id);
00064 }
00065 
00066 bool HcalTopology::validHcal(const HcalDetId& id) const {
00067   // check the raw rules
00068   bool ok=validRaw(id);
00069 
00070   ok=ok && !isExcluded(id);
00071 
00072   return ok;
00073 }
00074 
00075 bool HcalTopology::isExcluded(const HcalDetId& id) const {
00076   bool exed=false;
00077   // first, check the full detector exclusions...  (fast)
00078   switch (id.subdet()) {
00079   case(HcalBarrel):  exed=excludeHB_; break;
00080   case(HcalEndcap):  exed=excludeHE_; break;
00081   case(HcalOuter):   exed=excludeHO_; break;
00082   case(HcalForward): exed=excludeHF_; break;
00083   default: exed=false;
00084   }
00085   // next, check the list (slower)
00086   if (!exed && !exclusionList_.empty()) {
00087     std::vector<HcalDetId>::const_iterator i=std::lower_bound(exclusionList_.begin(),exclusionList_.end(),id);
00088     if (i!=exclusionList_.end() && *i==id) exed=true;
00089   }
00090   return exed;
00091 }
00092 
00093 void HcalTopology::exclude(const HcalDetId& id) {
00094   std::vector<HcalDetId>::iterator i=std::lower_bound(exclusionList_.begin(),exclusionList_.end(),id);
00095   if (i==exclusionList_.end() || *i!=id) {
00096     exclusionList_.insert(i,id);
00097   }
00098 }
00099 
00100 void HcalTopology::excludeSubdetector(HcalSubdetector subdet) {
00101   switch (subdet) {
00102   case(HcalBarrel):  excludeHB_=true; break;
00103   case(HcalEndcap):  excludeHE_=true; break;
00104   case(HcalOuter):   excludeHO_=true; break;
00105   case(HcalForward): excludeHF_=true; break;
00106   default: break;
00107   }
00108 }
00109 
00110 std::vector<DetId> HcalTopology::east(const DetId& id) const {
00111   std::vector<DetId> vNeighborsDetId;
00112   HcalDetId neighbors[2];
00113   for (int i=0;i<decIEta(HcalDetId(id),neighbors);i++)
00114     vNeighborsDetId.push_back(DetId(neighbors[i].rawId()));
00115   return vNeighborsDetId;
00116 }
00117 
00118 std::vector<DetId> HcalTopology::west(const DetId& id) const {
00119   std::vector<DetId> vNeighborsDetId;
00120   HcalDetId neighbors[2];
00121   for (int i=0;i<incIEta(HcalDetId(id),neighbors);i++)
00122     vNeighborsDetId.push_back(DetId(neighbors[i].rawId()));
00123   return  vNeighborsDetId;
00124 }
00125 
00126 std::vector<DetId> HcalTopology::north(const DetId& id) const {
00127   std::vector<DetId> vNeighborsDetId;
00128   HcalDetId neighbor;
00129   if (incIPhi(HcalDetId(id),neighbor))
00130     vNeighborsDetId.push_back(DetId(neighbor.rawId()));
00131   return  vNeighborsDetId;
00132 }
00133 
00134 std::vector<DetId> HcalTopology::south(const DetId& id) const {
00135   std::vector<DetId> vNeighborsDetId;
00136   HcalDetId neighbor;
00137   if (decIPhi(HcalDetId(id),neighbor))
00138     vNeighborsDetId.push_back(DetId(neighbor.rawId()));
00139   return  vNeighborsDetId;
00140 }
00141 
00142 std::vector<DetId> HcalTopology::up(const DetId& id) const {
00143   HcalDetId neighbor = id;
00144   //A.N.
00145   //  incrementDepth(neighbor);
00146   std::vector<DetId> vNeighborsDetId;
00147   if(incrementDepth(neighbor)) 
00148   {
00149     vNeighborsDetId.push_back(neighbor);
00150   }
00151   return  vNeighborsDetId;
00152 }
00153 
00154 std::vector<DetId> HcalTopology::down(const DetId& id) const {
00155   std::cout << "HcalTopology::down() not yet implemented" << std::endl; 
00156   std::vector<DetId> vNeighborsDetId;
00157   return  vNeighborsDetId;
00158 }
00159 
00160 int HcalTopology::exclude(HcalSubdetector subdet, int ieta1, int ieta2, int iphi1, int iphi2, int depth1, int depth2) {
00161 
00162   bool exed=false;
00163   // first, check the full detector exclusions...  (fast)
00164   switch (subdet) {
00165   case(HcalBarrel): exed=excludeHB_; break;
00166   case(HcalEndcap): exed=excludeHE_; break;
00167   case(HcalOuter): exed=excludeHO_; break;
00168   case(HcalForward): exed=excludeHF_; break;
00169   default: exed=false;
00170   }
00171   if (exed) return 0; // if the whole detector is excluded...
00172 
00173   int ieta_l=std::min(ieta1,ieta2);
00174   int ieta_h=std::max(ieta1,ieta2);
00175   int iphi_l=std::min(iphi1,iphi2);
00176   int iphi_h=std::max(iphi1,iphi2);
00177   int depth_l=std::min(depth1,depth2);
00178   int depth_h=std::max(depth1,depth2);
00179 
00180   int n=0;
00181   for (int ieta=ieta_l; ieta<=ieta_h; ieta++) 
00182     for (int iphi=iphi_l; iphi<=iphi_h; iphi++) 
00183       for (int depth=depth_l; depth<=depth_h; depth++) {
00184         HcalDetId id(subdet,ieta,iphi,depth);
00185         if (validRaw(id)) { // use 'validRaw' to include check validity in "uncut" detector
00186           exclude(id);  
00187           n++;
00188         }
00189       }
00190   return n;
00191 }
00192 
00219 bool HcalTopology::validDetIdPreLS1(const HcalDetId& id) const {
00220   const HcalSubdetector sd (id.subdet());
00221   const int             ie (id.ietaAbs());
00222   const int             ip (id.iphi());
00223   const int             dp (id.depth());
00224 
00225   return ( ( ip >=  1         ) &&
00226            ( ip <= IPHI_MAX   ) &&
00227            ( dp >=  1         ) &&
00228            ( ie >=  1         ) &&
00229            ( ( ( sd == HcalBarrel ) &&
00230                ( ( ( ie <= 14         ) &&
00231                    ( dp ==  1         )    ) ||
00232                  ( ( ( ie == 15 ) || ( ie == 16 ) ) && 
00233                    ( dp <= 2          )                ) ) ) ||
00234              (  ( sd == HcalEndcap ) &&
00235                 ( ( ( ie == firstHERing() ) &&
00236                     ( dp ==  3 )          ) ||
00237                   ( ( ie == 17 ) &&
00238                     ( dp ==  1 )          ) ||
00239                   ( ( ie >= 18 ) &&
00240                     ( ie <= 20 ) &&
00241                     ( dp <=  2 )          ) ||
00242                   ( ( ie >= 21 ) &&
00243                     ( ie <= 26 ) &&
00244                     ( dp <=  2 ) &&
00245                     ( ip%2 == 1 )         ) ||
00246                   ( ( ie >= 27 ) &&
00247                     ( ie <= 28 ) &&
00248                     ( dp <=  3 ) &&
00249                     ( ip%2 == 1 )         ) ||
00250                   ( ( ie == 29 ) &&
00251                     ( dp <=  2 ) &&
00252                     ( ip%2 == 1 )         )          )      ) ||
00253              (  ( sd == HcalOuter ) &&
00254                 ( ie <= 15 ) &&
00255                 ( dp ==  4 )           ) ||
00256              (  ( sd == HcalForward ) &&
00257                 ( dp <=  2 )          &&
00258                 ( ( ( ie >= firstHFRing() ) &&
00259                     ( ie <  firstHFQuadPhiRing() ) &&
00260                     ( ip%2 == 1 )    ) ||
00261                   ( ( ie >= firstHFQuadPhiRing() ) &&
00262                     ( ie <= lastHFRing() ) &&
00263                     ( ip%4 == 3 )         )  ) ) ) ) ;
00264 }
00265 
00267 bool HcalTopology::validRaw(const HcalDetId& id) const {
00268   bool ok=true;
00269   int ieta=id.ieta();
00270   int aieta=id.ietaAbs();
00271   int depth=id.depth();
00272   int iphi=id.iphi();
00273   if ((ieta==0 || iphi<=0 || iphi>IPHI_MAX) || aieta>lastHFRing()) return false; // outer limits
00274  
00275   if (ok) {
00276     HcalSubdetector subdet=id.subdet();
00277     if (subdet==HcalBarrel) {
00278       if (mode_==HcalTopologyMode::SLHC || mode_==HcalTopologyMode::H2HE) {
00279         if ((aieta>lastHBRing() || depth>maxDepthHB_ || (aieta==lastHBRing() && depth > 2))) ok=false;
00280       } else {
00281         if (aieta>lastHBRing() || depth>2 || (aieta<=14 && depth>1)) ok=false;
00282       }
00283     } else if (subdet==HcalEndcap) {
00284       if (mode_==HcalTopologyMode::SLHC || mode_==HcalTopologyMode::H2HE) {
00285         if (depth>maxDepthHE_ || aieta<firstHERing() || aieta>lastHERing() || (aieta==firstHERing() && depth<3) || (aieta>=firstHEDoublePhiRing() && (iphi%2)==0)) ok=false;
00286       } else {
00287         if (depth>3 || aieta<firstHERing() || aieta>lastHERing() || (aieta==firstHERing() && depth!=3) || (aieta==17 && depth!=1 && mode_!=HcalTopologyMode::H2) || // special case at H2
00288             (((aieta>=17 && aieta<firstHETripleDepthRing()) || aieta==lastHERing()) && depth>2) ||
00289             (aieta>=firstHEDoublePhiRing() && (iphi%2)==0)) ok=false;
00290       }
00291     } else if (subdet==HcalOuter) {
00292       if (aieta>lastHORing() || iphi>IPHI_MAX || depth!=4) ok=false;
00293     } else if (subdet==HcalForward) {
00294       if (aieta<firstHFRing() || aieta>lastHFRing() || ((iphi%2)==0) || (depth>2) ||  (aieta>=firstHFQuadPhiRing() && ((iphi+1)%4)!=0)) ok=false;
00295     } else {
00296       ok=false;
00297     }
00298   }
00299   return ok;
00300 }
00301 
00302   
00303 bool HcalTopology::incIPhi(const HcalDetId& id, HcalDetId &neighbor) const {
00304   bool ok=valid(id);
00305   if (ok) {
00306     switch (id.subdet()) {
00307     case (HcalBarrel):
00308     case (HcalOuter):
00309       if (id.iphi()==IPHI_MAX) neighbor=HcalDetId(id.subdet(),id.ieta(),1,id.depth()); 
00310       else neighbor=HcalDetId(id.subdet(),id.ieta(),id.iphi()+1,id.depth()); 
00311       break;
00312     case (HcalEndcap):
00313       if (id.ietaAbs()>=firstHEDoublePhiRing()) {
00314         if (id.iphi()==IPHI_MAX-1) neighbor=HcalDetId(id.subdet(),id.ieta(),1,id.depth()); 
00315         else neighbor=HcalDetId(id.subdet(),id.ieta(),id.iphi()+2,id.depth()); 
00316       } else {
00317         if (id.iphi()==IPHI_MAX) neighbor=HcalDetId(id.subdet(),id.ieta(),1,id.depth()); 
00318         else neighbor=HcalDetId(id.subdet(),id.ieta(),id.iphi()+1,id.depth()); 
00319       } 
00320       break;
00321     case (HcalForward):
00322       if (id.ietaAbs()>=firstHFQuadPhiRing()) {
00323         if (id.iphi()==IPHI_MAX-1) neighbor=HcalDetId(id.subdet(),id.ieta(),3,id.depth()); 
00324         else neighbor=HcalDetId(id.subdet(),id.ieta(),id.iphi()+4,id.depth()); 
00325       } else {
00326         if (id.iphi()==IPHI_MAX-1) neighbor=HcalDetId(id.subdet(),id.ieta(),1,id.depth()); 
00327         else neighbor=HcalDetId(id.subdet(),id.ieta(),id.iphi()+2,id.depth()); 
00328       }
00329       break;
00330     default: ok=false;
00331     }
00332   } 
00333   return ok;
00334 }
00335 
00337 bool HcalTopology::decIPhi(const HcalDetId& id, HcalDetId &neighbor) const {
00338   bool ok=valid(id);
00339   if (ok) {
00340     switch (id.subdet()) {
00341     case (HcalBarrel):
00342     case (HcalOuter):
00343       if (id.iphi()==1) neighbor=HcalDetId(id.subdet(),id.ieta(),IPHI_MAX,id.depth()); 
00344       else neighbor=HcalDetId(id.subdet(),id.ieta(),id.iphi()-1,id.depth()); 
00345       break;
00346     case (HcalEndcap):
00347       if (id.ietaAbs()>=firstHEDoublePhiRing()) {
00348         if (id.iphi()==1) neighbor=HcalDetId(id.subdet(),id.ieta(),IPHI_MAX-1,id.depth()); 
00349         else neighbor=HcalDetId(id.subdet(),id.ieta(),id.iphi()-2,id.depth()); 
00350       } else {
00351         if (id.iphi()==1) neighbor=HcalDetId(id.subdet(),id.ieta(),IPHI_MAX,id.depth()); 
00352         else neighbor=HcalDetId(id.subdet(),id.ieta(),id.iphi()-1,id.depth()); 
00353       }
00354       break;
00355     case (HcalForward):
00356       if (id.ietaAbs()>=firstHFQuadPhiRing()) {
00357         if (id.iphi()==3) neighbor=HcalDetId(id.subdet(),id.ieta(),IPHI_MAX-1,id.depth()); 
00358         else neighbor=HcalDetId(id.subdet(),id.ieta(),id.iphi()-4,id.depth()); 
00359       } else {
00360         if (id.iphi()==1) neighbor=HcalDetId(id.subdet(),id.ieta(),IPHI_MAX-1,id.depth()); 
00361         else neighbor=HcalDetId(id.subdet(),id.ieta(),id.iphi()-2,id.depth()); 
00362       }
00363       break;
00364     default: ok=false;
00365     }
00366   } 
00367   return ok;
00368 }
00369 
00370 int HcalTopology::incIEta(const HcalDetId& id, HcalDetId neighbors[2]) const {
00371   if (id.zside()==1) return incAIEta(id,neighbors);
00372   else return decAIEta(id,neighbors);
00373 }
00374 
00375 int HcalTopology::decIEta(const HcalDetId& id, HcalDetId neighbors[2]) const {
00376   if (id.zside()==1) return decAIEta(id,neighbors);
00377   else return incAIEta(id,neighbors);
00378 }
00379 
00381 int HcalTopology::incAIEta(const HcalDetId& id, HcalDetId neighbors[2]) const {
00382   int n=1;
00383   int aieta=id.ietaAbs();
00384 
00385   if (aieta==firstHEDoublePhiRing()-1 && (id.iphi()%2)==0) 
00386     neighbors[0]=HcalDetId(id.subdet(),(aieta+1)*id.zside(),id.iphi()-1,id.depth());
00387   else if (aieta==firstHFQuadPhiRing()-1 && ((id.iphi()+1)%4)!=0) 
00388     neighbors[0]=HcalDetId(id.subdet(),(aieta+1)*id.zside(),((id.iphi()==1)?(71):(id.iphi()-2)),id.depth());
00389   else if (aieta==lastHBRing()) 
00390     neighbors[0]=HcalDetId(HcalEndcap,(aieta+1)*id.zside(),id.iphi(),1);
00391   else if (aieta==lastHERing()) 
00392     neighbors[0]=HcalDetId(HcalForward,(aieta+1)*id.zside(),id.iphi(),1);
00393   else
00394     neighbors[0]=HcalDetId(id.subdet(),(aieta+1)*id.zside(),id.iphi(),id.depth());
00395     
00396   if (!valid(neighbors[0])) n=0;
00397   return n;
00398 }
00399 
00401 int HcalTopology::decAIEta(const HcalDetId& id, HcalDetId neighbors[2]) const {
00402   int n=1;
00403   int aieta=id.ietaAbs();
00404 
00405   if (aieta==firstHEDoublePhiRing()) { 
00406     n=2;
00407     neighbors[0]=HcalDetId(id.subdet(),(aieta-1)*id.zside(),id.iphi(),id.depth());
00408     neighbors[1]=HcalDetId(id.subdet(),(aieta-1)*id.zside(),id.iphi()+1,id.depth());
00409   } else if (aieta==firstHFQuadPhiRing()) {
00410     n=2;
00411     neighbors[0]=HcalDetId(id.subdet(),(aieta-1)*id.zside(),id.iphi(),id.depth());
00412     if (id.iphi()==IPHI_MAX-1) neighbors[1]=HcalDetId(id.subdet(),(aieta-1)*id.zside(),1,id.depth());
00413     else neighbors[1]=HcalDetId(id.subdet(),(aieta-1)*id.zside(),id.iphi()+2,id.depth());
00414   } else if (aieta==1) {
00415     neighbors[0]=HcalDetId(id.subdet(),-aieta*id.zside(),id.iphi(),id.depth());
00416   } else if (aieta==lastHBRing()+1) {
00417     neighbors[0]=HcalDetId(HcalBarrel,(aieta-1)*id.zside(),id.iphi(),id.depth());
00418   } else if (aieta==lastHERing()+1) {
00419     neighbors[0]=HcalDetId(HcalEndcap,(aieta-1)*id.zside(),id.iphi(),id.depth());
00420   } else
00421     neighbors[0]=HcalDetId(id.subdet(),(aieta-1)*id.zside(),id.iphi(),id.depth());
00422   
00423   if (!valid(neighbors[0]) && n==2) {
00424     if (!valid(neighbors[1])) n=0;
00425     else {
00426       n=1;
00427       neighbors[0]=neighbors[1];
00428     }
00429   }
00430   if (n==2 && !valid(neighbors[1])) n=1;
00431   if (n==1 && !valid(neighbors[0])) n=0;
00432 
00433   return n;
00434 }
00435 
00436 
00437 void HcalTopology::depthBinInformation(HcalSubdetector subdet, int etaRing,
00438                                        int & nDepthBins, int & startingBin) const {
00439 
00440   if(subdet == HcalBarrel) {
00441     if (mode_==HcalTopologyMode::SLHC || mode_==HcalTopologyMode::H2HE) {
00442       startingBin = 1;
00443       if (etaRing==lastHBRing()) {
00444         nDepthBins = 2;
00445       } else {
00446         nDepthBins = maxDepthHB_;
00447       }
00448     } else {
00449       if (etaRing<=14) {
00450         nDepthBins = 1;
00451         startingBin = 1;
00452       } else {
00453         nDepthBins = 2;
00454         startingBin = 1;
00455       }
00456     }
00457   } else if(subdet == HcalEndcap) {
00458     if (mode_==HcalTopologyMode::SLHC || mode_==HcalTopologyMode::H2HE) {
00459       if (etaRing==firstHERing()) {
00460         nDepthBins  = maxDepthHE_ - 2;
00461         startingBin = 3;
00462       } else {
00463         nDepthBins  = maxDepthHE_;
00464         startingBin = 1;
00465       }
00466     } else {
00467       if (etaRing==firstHERing()) {
00468         nDepthBins = 1;
00469         startingBin = 3;
00470       } else if (etaRing==17) {
00471         nDepthBins = 1;
00472         startingBin = 1;
00473       } else if (etaRing==lastHERing()) {
00474         nDepthBins = 2;
00475         startingBin = 1;
00476       } else {
00477         nDepthBins = (etaRing >= firstHETripleDepthRing()) ? 3 : 2;
00478         startingBin = 1;
00479       }
00480     }
00481   } else if(subdet == HcalForward) {
00482     nDepthBins  = 2;
00483     startingBin = 1;
00484   } else if(subdet == HcalOuter) {
00485     nDepthBins = 1;
00486     startingBin = 4;
00487   } else {
00488     std::cerr << "Bad HCAL subdetector " << subdet << std::endl;
00489   }
00490 }
00491 
00492 
00493 bool HcalTopology::incrementDepth(HcalDetId & detId) const
00494 {
00495   HcalSubdetector subdet = detId.subdet();
00496   int ieta = detId.ieta();
00497   int etaRing = detId.ietaAbs();
00498   int depth = detId.depth();
00499   int nDepthBins, startingBin;
00500   depthBinInformation(subdet, etaRing, nDepthBins, startingBin);
00501 
00502   // see if the new depth bin exists
00503   ++depth;
00504   if(depth > nDepthBins) {
00505     // handle on a case-by-case basis
00506     if(subdet == HcalBarrel && etaRing < lastHORing())  {
00507       // HO
00508       subdet = HcalOuter;
00509       depth = 4;
00510     } else if(subdet == HcalBarrel && etaRing == lastHBRing()) {
00511       // overlap
00512       subdet = HcalEndcap;
00513     } else if(subdet == HcalEndcap && etaRing ==  lastHERing()-1) {
00514       // guard ring HF29 is behind HE 28
00515       subdet = HcalForward;
00516       (ieta > 0) ? ++ieta : --ieta;
00517       depth = 1;
00518     } else if(subdet == HcalEndcap && etaRing ==  lastHERing()) {
00519       // split cells go to bigger granularity.  Ring 29 -> 28
00520       (ieta > 0) ? --ieta : ++ieta;
00521     } else {
00522       // no more chances
00523       detId = HcalDetId();
00524       return false;
00525     }
00526   }
00527   detId = HcalDetId(subdet, ieta, detId.iphi(), depth);
00528   //A.N.  
00529   // assert(validRaw(detId));
00530   return validRaw(detId);
00531   //A.N.  return true;
00532 }
00533 
00534 
00535 int HcalTopology::nPhiBins(int etaRing) const {
00536   int lastPhiBin=singlePhiBins_;
00537   if (etaRing>= firstHFQuadPhiRing()) lastPhiBin=doublePhiBins_/2;
00538   else if (etaRing>= firstHEDoublePhiRing()) lastPhiBin=doublePhiBins_;
00539   return lastPhiBin;
00540 }
00541 
00542 void HcalTopology::getDepthSegmentation(unsigned ring, std::vector<int> & readoutDepths) const
00543 {
00544   // if it doesn't exist, return the first entry with a lower index.  So if we only
00545   // have entries for 1 and 17, any input from 1-16 should return the entry for ring 1
00546   SegmentationMap::const_iterator pos = depthSegmentation_.upper_bound(ring);
00547   if(pos == depthSegmentation_.begin()) {
00548     throw cms::Exception("HcalTopology") << "No depth segmentation found for ring" << ring;
00549   }
00550   --pos;
00551     // pos now refers to the last element with key <= ring.
00552   readoutDepths = pos->second;
00553 }
00554 
00555 void HcalTopology::setDepthSegmentation(unsigned ring, const std::vector<int> & readoutDepths)
00556 {
00557   depthSegmentation_[ring] = readoutDepths;
00558 }
00559 
00560 std::pair<int, int> HcalTopology::segmentBoundaries(unsigned ring, unsigned depth) const {
00561   std::vector<int> readoutDepths;
00562   getDepthSegmentation(ring, readoutDepths);
00563   int d1 = std::lower_bound(readoutDepths.begin(), readoutDepths.end(), depth) - readoutDepths.begin();
00564   int d2 = std::upper_bound(readoutDepths.begin(), readoutDepths.end(), depth) - readoutDepths.begin();
00565   return std::pair<int, int>(d1, d2);
00566 }
00567 
00568 unsigned int HcalTopology::detId2denseIdPreLS1 (const DetId& id) const {
00569 
00570   HcalDetId hid(id);
00571   const HcalSubdetector sd (hid.subdet()  ) ;
00572   const int             ip (hid.iphi()    ) ;
00573   const int             ie (hid.ietaAbs() ) ;
00574   const int             dp (hid.depth()   ) ;
00575   const int             zn (hid.zside() < 0 ? 1 : 0 ) ;
00576   unsigned int  retval = ( ( sd == HcalBarrel ) ?
00577                            ( ip - 1 )*18 + dp - 1 + ie - ( ie<16 ? 1 : 0 ) + zn*kHBhalf :
00578                            ( ( sd == HcalEndcap ) ?
00579                              2*kHBhalf + ( ip - 1 )*8 + ( ip/2 )*20 +
00580                              ( ( ie==16 || ie==17 ) ? ie - 16 :
00581                                ( ( ie>=18 && ie<=20 ) ? 2 + 2*( ie - 18 ) + dp - 1 :
00582                                  ( ( ie>=21 && ie<=26 ) ? 8 + 2*( ie - 21 ) + dp - 1 :
00583                                    ( ( ie>=27 && ie<=28 ) ? 20 + 3*( ie - 27 ) + dp - 1 :
00584                                      26 + 2*( ie - 29 ) + dp - 1 ) ) ) ) + zn*kHEhalf :
00585                              ( ( sd == HcalOuter ) ?
00586                                2*kHBhalf + 2*kHEhalf + ( ip - 1 )*15 + ( ie - 1 ) + zn*kHOhalf :
00587                                ( ( sd == HcalForward ) ?
00588                                  2*kHBhalf + 2*kHEhalf + 2*kHOhalf + 
00589                                  ( ( ip - 1 )/4 )*4 + ( ( ip - 1 )/2 )*22 + 
00590                                  2*( ie - 29 ) + ( dp - 1 ) + zn*kHFhalf : 0xFFFFFFFFu ) ) ) ) ; 
00591   return retval;
00592 }
00593 
00594 
00595 unsigned int HcalTopology::detId2denseIdHB(const DetId& id) const {
00596   HcalDetId hid(id);
00597   const int             ip (hid.iphi()    ) ;
00598   const int             ie (hid.ietaAbs() ) ;
00599   const int             dp (hid.depth()   ) ;
00600   const int             zn (hid.zside() < 0 ? 1 : 0 ) ;
00601   unsigned int  retval = 0xFFFFFFFFu;
00602   if (topoVersion_==0) {
00603     retval=( ip - 1 )*18 + dp - 1 + ie - ( ie<16 ? 1 : 0 ) + zn*kHBhalf;
00604   } else if (topoVersion_==10) {
00605     retval=dp-1 + maxDepthHB_*(ip-1)+maxDepthHB_*72*(hid.ieta()-1+33*zn);
00606   }
00607   return retval;
00608 }
00609 
00610 unsigned int HcalTopology::detId2denseIdHE(const DetId& id) const {
00611   HcalDetId hid(id);
00612   const int             ip (hid.iphi()    ) ;
00613   const int             ie (hid.ietaAbs() ) ;
00614   const int             dp (hid.depth()   ) ;
00615   const int             zn (hid.zside() < 0 ? 1 : 0 ) ;
00616   unsigned int  retval =  0xFFFFFFFFu;
00617   if (topoVersion_==0) {
00618     retval=( ip - 1 )*8 + ( ip/2 )*20 +
00619       ( ( ie==16 || ie==17 ) ? ie - 16 :
00620         ( ( ie>=18 && ie<=20 ) ? 2 + 2*( ie - 18 ) + dp - 1 :
00621           ( ( ie>=21 && ie<=26 ) ? 8 + 2*( ie - 21 ) + dp - 1 :
00622             ( ( ie>=27 && ie<=28 ) ? 20 + 3*( ie - 27 ) + dp - 1 :
00623               26 + 2*( ie - 29 ) + dp - 1 ) ) ) ) + zn*kHEhalf;
00624   } else if (topoVersion_==10) {
00625     retval=(dp-1)+maxDepthHE_*(ip-1)+maxDepthHE_*72*(hid.ieta()-16+zn*(14+29+16));
00626   }
00627   return retval;
00628 }
00629 
00630 unsigned int HcalTopology::detId2denseIdHO(const DetId& id) const {
00631   HcalDetId hid(id);
00632   const int             ip (hid.iphi()    ) ;
00633   const int             ie (hid.ietaAbs() ) ;
00634   const int             zn (hid.zside() < 0 ? 1 : 0 ) ;
00635 
00636   unsigned int  retval = 0xFFFFFFFFu;
00637   if (topoVersion_==0) {
00638     retval=( ip - 1 )*15 + ( ie - 1 ) + zn*kHOhalf;
00639   } else if (topoVersion_==10) {
00640     retval=( ip - 1 )*15 + ( ie - 1 ) + zn*kHOhalf;
00641   }
00642   return retval;
00643 }
00644 
00645 unsigned int HcalTopology::detId2denseIdHF(const DetId& id) const {
00646   HcalDetId hid(id);
00647   const int             ip (hid.iphi()    ) ;
00648   const int             ie (hid.ietaAbs() ) ;
00649   const int             dp (hid.depth()   ) ;
00650   const int             zn (hid.zside() < 0 ? 1 : 0 ) ;
00651 
00652   unsigned int  retval = 0xFFFFFFFFu;
00653   if (topoVersion_==0) {
00654     retval = ( ( ip - 1 )/4 )*4 + ( ( ip - 1 )/2 )*22 + 
00655       2*( ie - 29 ) + ( dp - 1 ) + zn*kHFhalf;
00656   } else if (topoVersion_==10) {
00657     retval = ( ( ip - 1 )/4 )*4 + ( ( ip - 1 )/2 )*22 + 
00658       2*( ie - 29 ) + ( dp - 1 ) + zn*kHFhalf;
00659   }
00660   return retval;
00661 }
00662 
00663 unsigned int HcalTopology::detId2denseIdHT(const DetId& id) const {
00664   HcalTrigTowerDetId tid(id); 
00665   int zside = tid.zside();
00666   unsigned int ietaAbs = tid.ietaAbs();
00667   unsigned int iphi = tid.iphi();
00668 
00669   unsigned int index;
00670   if ((iphi-1)%4==0) index = (iphi-1)*32 + (ietaAbs-1) - (12*((iphi-1)/4));
00671   else               index = (iphi-1)*28 + (ietaAbs-1) + (4*(((iphi-1)/4)+1));
00672   
00673   if (zside == -1) index += kHThalf;
00674 
00675   return index;
00676 }
00677 
00678 unsigned int HcalTopology::detId2denseIdCALIB(const DetId& id) const {
00679   HcalCalibDetId tid(id);
00680   int    channel = tid.cboxChannel();
00681   int ieta = tid.ieta();
00682   int iphi = tid.iphi();
00683   int zside = tid.zside();
00684   unsigned int index=0xFFFFFFFFu;
00685       
00686   if (tid.calibFlavor()==HcalCalibDetId::CalibrationBox) {
00687         
00688     HcalSubdetector subDet = tid.hcalSubdet();
00689         
00690     if (subDet==HcalBarrel) {
00691       //std::cout<<"CALIB_HB:  ";
00692       //dphi = 4 (18 phi values), 3 channel types (0,1,2), eta = -1 or 1
00693       //total of 18*3*2=108 channels
00694       index = ((iphi+1)/4-1) + 18*channel + 27*(ieta+1);
00695     } else if (subDet==HcalEndcap) {
00696       //std::cout<<"CALIB_HE:  ";
00697       //dphi = 4 (18 phi values), 6 channel types (0,1,3,4,5,6), eta = -1 or 1
00698       //total of 18*6*2=216 channels
00699       if (channel>2) channel-=1;
00700       index = ((iphi+1)/4-1) + 18*channel + 54*(ieta+1) + 108;
00701     } else if (subDet==HcalForward) {
00702       //std::cout<<"CALIB_HF:  ";
00703       //dphi = 18 (4 phi values), 3 channel types (0,1,8), eta = -1 or 1
00704       if (channel==8) channel = 2;
00705       //total channels 4*3*2=24
00706       index = (iphi-1)/18 + 4*channel + 6*(ieta+1) + 324;
00707     } else if (subDet==HcalOuter) {
00708       //std::cout<<"CALIB_HO:  ";
00709       //there are 5 special calib crosstalk channels, one in each ring
00710       if (channel==7) {
00711         channel = 2;
00712         index = (ieta+2) + 420;
00713       //for HOM/HOP dphi = 6 (12 phi values),  2 channel types (0,1), eta = -2,-1 or 1,2
00714       //for HO0/YB0 dphi = 12 (6 phi values),  2 channel types (0,1), eta = 0
00715       } else{
00716         if (ieta<0) index      = ((iphi+1)/12-1) + 36*channel + 6*(ieta+2) + 348;
00717         else if (ieta>0) index = ((iphi+1)/12-1) + 36*channel + 6*(ieta+2) + 6 + 348;
00718         else index             = ((iphi+1)/6-1)  + 36*channel + 6*(ieta+2) + 348;
00719       }
00720     } else {
00721       std::cout << "HCAL Det Id not valid!" << std::endl;
00722       index = 0;
00723     }
00724         
00725   } else if (tid.calibFlavor()==HcalCalibDetId::HOCrosstalk) {
00726     //std::cout<<"HX:  ";
00727     //for YB0/HO0 phi is grouped in 6 groups of 6 with dphi=2 but the transitions are 1 or 3
00728     // in such a way that the %36 operation yeilds unique values for every iphi
00729     if (abs(ieta)==4)  index = ((iphi-1)%36) + (((zside+1)*36)/2) + 72 + 425;   //ieta = 1 YB0/HO0;
00730     else               index = (iphi-1) + (36*(zside+1)*2) + 425;  //ieta = 0 for HO2M/HO1M ieta=2 for HO1P/HO2P;
00731   }
00732   //std::cout << "  " << ieta << "  " << zside << "  " << iphi << "  " << depth << "  " << index << std::endl;
00733   return index;
00734 }
00735 
00736 
00737 unsigned int HcalTopology::detId2denseId(const DetId& id) const {
00738   unsigned int retval(0);
00739   if (topoVersion_==0) { // pre-LS1
00740     retval = detId2denseIdPreLS1(id);
00741   } else if (topoVersion_==10) {
00742     HcalDetId hid(id);
00743     if (hid.subdet()==HcalBarrel) {
00744       retval=(hid.depth()-1)+maxDepthHB_*(hid.iphi()-1);
00745       if (hid.ieta()>0) {
00746         retval+=maxDepthHB_*72*(hid.ieta()-1);
00747       } else {
00748         retval+=maxDepthHB_*72*(32+hid.ieta());
00749       }
00750     } else if (hid.subdet()==HcalEndcap) {
00751       retval=HBSize_;
00752       retval+=(hid.depth()-1)+maxDepthHE_*(hid.iphi()-1);
00753       if (hid.ieta()>0) {
00754         retval+=maxDepthHE_*72*(hid.ieta()-16);
00755       } else {
00756         retval+=maxDepthHE_*72*((14+29)+hid.ieta());
00757       }
00758     } else if (hid.subdet()==HcalOuter) {
00759       retval=HBSize_+HESize_;
00760       if   (hid.ieta()>0) retval+=(hid.iphi()-1)+72*(hid.ieta()-1);
00761       else retval+=(hid.iphi()-1)+72*(30+hid.ieta());
00762     } else if (hid.subdet()==HcalForward) { 
00763       retval=HBSize_+HESize_+HOSize_;
00764       retval+=hid.depth()-1+2*(hid.iphi()-1);
00765       if (hid.ieta()>0) retval+=2*72*(hid.ieta()-29);
00766       else retval+=2*72*((41+13)+hid.ieta());
00767     } else {
00768       return 0xFFFFFFFu;
00769     }
00770   }
00771   return retval;
00772 }
00773 
00774 DetId HcalTopology::denseId2detId(unsigned int denseid) const {
00775 
00776   HcalSubdetector sd ( HcalBarrel ) ;
00777   int ie ( 0 ) ;
00778   int ip ( 0 ) ;
00779   int dp ( 0 ) ;
00780   int in ( denseid ) ;
00781   int iz ( 1 ) ;
00782   if (topoVersion_==0) { //DL// pre-LS1
00783     if (denseid < kSizeForDenseIndexingPreLS1) {
00784       if ( in > 2*( kHBhalf + kHEhalf + kHOhalf ) - 1 ) { // HF
00785         sd  = HcalForward ;
00786         in -= 2*( kHBhalf + kHEhalf + kHOhalf ) ; 
00787         iz  = ( in<kHFhalf ? 1 : -1 ) ;
00788         in %= kHFhalf ; 
00789         ip  = 4*( in/48 ) ;
00790         in %= 48 ;
00791         ip += 1 + ( in>21 ? 2 : 0 ) ;
00792         if( 3 == ip%4 ) in -= 22 ;
00793         ie  = 29 + in/2 ;
00794         dp  = 1 + in%2 ;
00795       } else if ( in > 2*( kHBhalf + kHEhalf ) - 1 ) { // HO
00796         sd  = HcalOuter ;
00797         in -= 2*( kHBhalf + kHEhalf ) ; 
00798         iz  = ( in<kHOhalf ? 1 : -1 ) ;
00799         in %= kHOhalf ; 
00800         dp  = 4 ;
00801         ip  = 1 + in/15 ;
00802         ie  = 1 + ( in - 15*( ip - 1 ) ) ;
00803       } else if ( in > 2*kHBhalf - 1 ) { // Endcap
00804         sd  = HcalEndcap ;
00805         in -= 2*kHBhalf ;
00806         iz  = ( in<kHEhalf ? 1 : -1 ) ;
00807         in %= kHEhalf ; 
00808         ip  = 2*( in/36 ) ;
00809         in %= 36 ;
00810         ip += 1 + in/28 ;
00811         if( 0 == ip%2 ) in %= 28 ;
00812         ie  = 15 + ( in<2 ? 1 + in : 2 + 
00813                      ( in<20 ? 1 + ( in - 2 )/2 : 9 +
00814                        ( in<26 ? 1 + ( in - 20 )/3 : 3 ) ) ) ;
00815         dp  = ( in<1 ? 3 :
00816                 ( in<2 ? 1 : 
00817                   ( in<20 ? 1 + ( in - 2 )%2 : 
00818                     ( in<26 ? 1 + ( in - 20 )%3 : 
00819                       ( 1 + ( in - 26 )%2 ) ) ) ) ) ;
00820       } else { // barrel
00821         iz  = ( in<kHBhalf ? 1 : -1 ) ;
00822         in %= kHBhalf ; 
00823         ip = in/18 + 1 ;
00824         in %= 18 ;
00825         if ( in < 14 ) {
00826           dp = 1 ;
00827           ie = in + 1 ;
00828         } else {
00829           in %= 14 ;
00830           dp =  1 + in%2 ;
00831           ie = 15 + in/2 ;
00832         }
00833       }
00834     }
00835   } else if (topoVersion_==10) {
00836     if (denseid < ncells()) {
00837       if (denseid >= (HBSize_+HESize_+HOSize_)) {
00838         sd  = HcalForward ;
00839         in -= (HBSize_+HESize_+HOSize_);
00840         dp  = (in%2) + 1;
00841         ip  = (in - dp + 1)%144;
00842         ip  = (ip/2) + 1;
00843         ie  = (in - dp + 1 - 2*(ip -1))/144;
00844         if (ie > 12) {ie  = 54 -ie; iz = -1;}
00845         else         {ie += 29;     iz =  1;}
00846       } else if (denseid >= (HBSize_+HESize_)) {
00847         sd  = HcalOuter ;
00848         in -= (HBSize_+HESize_);
00849         dp  = 4;
00850         ip  = (in%72) + 1;
00851         ie  = (in - ip + 1)/72;
00852         if (ie > 14) {ie  = 30 -ie; iz = -1;}
00853         else         {ie += 1;      iz =  1;}
00854       } else if (denseid >= (HBSize_)) {
00855         sd  = HcalEndcap ;
00856         in -= (HBSize_);
00857         dp  = (in%maxDepthHE_)+1;
00858         ip  = (in - dp + 1)%(maxDepthHE_*72);
00859         ip  = (ip/maxDepthHE_) + 1;
00860         ie  = (in - dp + 1 - maxDepthHE_*(ip-1))/(72*maxDepthHE_);
00861         if (ie > 13) {ie  = 43 - ie; iz = -1;}
00862         else         {ie += 16;      iz =  1;}
00863       } else {
00864         sd  = HcalBarrel ;
00865         dp  = (in%maxDepthHB_)+1;
00866         ip  = (in - dp + 1)%(maxDepthHB_*72);
00867         ip  = (ip/maxDepthHB_) + 1;
00868         ie  = (in - dp + 1 - maxDepthHB_*(ip-1))/(72*maxDepthHB_);
00869         if (ie > 15) {ie  = 32 - ie; iz = -1;}
00870         else         {ie += 1;       iz =  1;}
00871       } 
00872     }
00873   }
00874   return HcalDetId( sd, iz*int(ie), ip, dp );
00875 }
00876 
00877 unsigned int HcalTopology::ncells() const {
00878   return HBSize_+HESize_+HOSize_+HFSize_;
00879 }
00880 
00881 int HcalTopology::topoVersion() const {
00882   return topoVersion_;
00883 }