Reference Manual

Go to the documentation of this file.
00001 
00002 // File: HcalNumberingFromDDD.cc
00003 // Description: Usage of DDD to get to numbering scheme for hadron calorimeter
00005 
00006 #include "Geometry/HcalCommonData/interface/HcalNumberingFromDDD.h"
00007 
00008 #include "DetectorDescription/Base/interface/DDException.h"
00009 #include "DetectorDescription/Core/interface/DDFilter.h"
00010 #include "DetectorDescription/Core/interface/DDFilteredView.h"
00011 #include "DetectorDescription/Core/interface/DDSolid.h"
00012 #include "DetectorDescription/Core/interface/DDValue.h"
00013 
00014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00015 
00016 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00017 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00018 #include <iostream>
00019 
00020 //#define DebugLog
00021 
00022 HcalNumberingFromDDD::HcalNumberingFromDDD(std::string & name,
00023                                            const DDCompactView & cpv) {
00024   edm::LogInfo("HCalGeom") << "Creating HcalNumberingFromDDD";
00025   initialize(name, cpv);
00026 }
00027 
00028 HcalNumberingFromDDD::~HcalNumberingFromDDD() {
00029   edm::LogInfo("HCalGeom") << "Deleting HcalNumberingFromDDD";
00030 }
00031 
00032 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det,
00033                                                           CLHEP::Hep3Vector point,
00034                                                           int depth,
00035                                                           int lay) const {
00036 
00037 
00038   double hx  = point.x();
00039   double hy  = point.y();
00040   double hz  = point.z();
00041   double hR  = sqrt(hx*hx+hy*hy+hz*hz);
00042   double htheta = (hR == 0. ? 0. : acos(std::max(std::min(hz/hR,1.0),-1.0)));
00043   double hsintheta = sin(htheta);
00044   double hphi = (hR*hsintheta == 0. ? 0. :atan2(hy,hx));
00045   double heta = (fabs(hsintheta) == 1.? 0. : -log(fabs(tan(htheta/2.))) );
00046 
00047   int    hsubdet=0;
00048   double etaR;
00049 
00050   //First eta index
00051   if (det == 5) { // Forward HCal
00052     hsubdet = static_cast<int>(HcalForward);
00053     hR      = sqrt(hx*hx+hy*hy);
00054     etaR    = (heta >= 0. ? hR : -hR);
00055   } else { // Barrel or Endcap
00056     etaR    = heta;
00057     if (det == 3) {
00058       hsubdet = static_cast<int>(HcalBarrel);
00059     } else {
00060       hsubdet = static_cast<int>(HcalEndcap);
00061     }
00062   }
00063 
00064 #ifdef DebugLog
00065   LogDebug("HCalGeom") << "HcalNumberingFromDDD: point = " << point << " det "
00066                        << hsubdet << " eta/R " << etaR << " phi " << hphi;
00067 #endif
00068   HcalNumberingFromDDD::HcalID tmp = unitID(hsubdet,etaR,hphi,depth,lay);
00069   return tmp;
00070 
00071 }
00072 
00073 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(double eta,double fi,
00074                                                           int depth, 
00075                                                           int lay) const {
00076 
00077   int    ieta = 0;
00078   double heta = fabs(eta);
00079   for (int i = 0; i < nEta; i++)
00080     if (heta > etaTable[i]) ieta = i + 1;
00081   int    hsubdet=0;
00082   double etaR;
00083   if (ieta <= etaMin[1]) {
00084     if ((ieta <= etaMin[1] && depth==3) || ieta > etaMax[0]) {
00085       hsubdet = static_cast<int>(HcalEndcap);
00086     } else {
00087       hsubdet = static_cast<int>(HcalBarrel);
00088     }
00089     etaR    = eta;
00090   } else {
00091     hsubdet = static_cast<int>(HcalForward);
00092     double theta = 2.*atan(exp(-heta));
00093     double hR    = zVcal*tan(theta);
00094     etaR    = (eta >= 0. ? hR : -hR);
00095   }
00096 
00097   HcalNumberingFromDDD::HcalID tmp = unitID(hsubdet,etaR,fi,depth,lay);
00098   return tmp;
00099 }
00100 
00101 
00102 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det,
00103                                                           double etaR,
00104                                                           double phi,
00105                                                           int depth,
00106                                                           int lay) const {
00107 
00108   int ieta=0;
00109   double fioff, fibin;
00110   double hetaR = fabs(etaR);
00111 
00112   //First eta index
00113   if (det == static_cast<int>(HcalForward)) { // Forward HCal
00114     fioff   = phioff[2];
00115     ieta    = etaMax[2];
00116     for (int i = nR-1; i > 0; i--)
00117       if (hetaR < rTable[i]) ieta = etaMin[2] + nR - i - 1;
00118     fibin   = phibin[nEta+ieta-etaMin[2]-1];
00119     if  (ieta > etaMax[2]-2 ) {   // HF double-phi  
00120       fioff += 0.5*fibin;
00121     }
00122   } else { // Barrel or Endcap
00123     ieta  = 1;
00124     for (int i = 0; i < nEta-1; i++)
00125       if (hetaR > etaTable[i]) ieta = i + 1;
00126     if (det == static_cast<int>(HcalBarrel)) {
00127       fioff   = phioff[0];
00128       if (ieta > etaMax[0])  ieta = etaMax[0];
00129       if (lay == 18 && nOff.size() > 13) {
00130         if (hetaR > etaHO[1] && ieta == nOff[13]) ieta++;
00131       }
00132     } else {
00133       fioff   = phioff[1];
00134       if (ieta <= etaMin[1]) ieta = etaMin[1];
00135     }
00136     fibin = phibin[ieta-1];
00137   }
00138 
00139   int    nphi  = int((CLHEP::twopi+0.1*fibin)/fibin);
00140   int    zside = etaR>0 ? 1: 0;
00141   double hphi  = phi+fioff;
00142   if (hphi < 0)    hphi += CLHEP::twopi;
00143   int    iphi  = int(hphi/fibin) + 1;
00144   if (iphi > nphi) iphi = 1;
00145 
00146 #ifdef DebugLog
00147   LogDebug("HCalGeom") << "HcalNumberingFromDDD: etaR = " << etaR << " : "
00148                        << zside << "/" << ieta << " phi " << hphi << " : "
00149                        << iphi;
00150 #endif
00151   HcalNumberingFromDDD::HcalID tmp = unitID(det,zside,depth,ieta,iphi,lay);
00152   return tmp;
00153 
00154 }
00155 
00156 HcalNumberingFromDDD::HcalID HcalNumberingFromDDD::unitID(int det, int zside,
00157                                                           int depth, int etaR,
00158                                                           int phi, 
00159                                                           int lay) const {
00160 
00161   //Modify the depth index
00162   if (det == static_cast<int>(HcalForward)) { // Forward HCal
00163   } else {
00164     if (lay >= 0) {
00165       double fibin = phibin[etaR-1];
00166       int   depth0 = depth1[etaR-1];
00167       int   kphi   = phi + int((phioff[3]+0.1)/fibin);
00168       kphi         = (kphi-1)%4 + 1;
00169       if (etaR == nOff[0] && (kphi == 2 || kphi == 3)) depth0--;
00170       if (lay <= depth2[etaR-1]) {
00171         if (lay <= depth0) depth = 1;
00172         else               depth = 2;
00173       } else if (lay <= depth3[etaR-1]) {
00174         depth = 3;
00175       } else               depth = 4;
00176     } else if (det == static_cast<int>(HcalBarrel)) {
00177       if (depth==3) depth = 2;
00178     }
00179     if (det != static_cast<int>(HcalBarrel)) {
00180       if (etaR <= etaMin[1]) depth = 3;
00181     }
00182   }
00183   if (etaR == nOff[1] && depth > 2 && det == static_cast<int>(HcalEndcap))
00184     etaR = nOff[1]-1;
00185   if (det == static_cast<int>(HcalBarrel) && depth == 4) {
00186     det = static_cast<int>(HcalOuter);
00187   }
00188 
00189   int units     = unitPhi(det, etaR);
00190   int iphi_skip = phi;
00191   if      (units==2) iphi_skip  = (phi-1)*2+1;
00192   else if (units==4) iphi_skip  = (phi-1)*4-1;
00193   if (iphi_skip < 0) iphi_skip += 72;
00194 
00195 #ifdef DebugLog
00196   LogDebug("HCalGeom") << "HcalNumberingFromDDD: phi units=" <<  units  
00197                        <<  "  iphi_skip=" << iphi_skip; 
00198 #endif
00199   HcalNumberingFromDDD::HcalID tmp(det,zside,depth,etaR,phi,iphi_skip,lay);
00200 
00201 #ifdef DebugLog
00202   LogDebug("HCalGeom") << "HcalNumberingFromDDD: det = " << det << " " 
00203                        << tmp.subdet << " zside = " << tmp.zside << " depth = "
00204                        << tmp.depth << " eta/R = " << tmp.etaR << " phi = " 
00205                        << tmp.phi << " layer = " << tmp.lay;
00206 #endif
00207   return tmp;
00208 }
00209 
00210 HcalCellType::HcalCell HcalNumberingFromDDD::cell(int det, int zside, 
00211                                                   int depth, int etaR,
00212                                                   int iphi, bool corr) const {
00213 
00214   int idet = det;
00215   double etaMn = etaMin[0];
00216   double etaMx = etaMax[0];
00217   if (idet==static_cast<int>(HcalEndcap)) {
00218     etaMn = etaMin[1]; etaMx = etaMax[1];
00219   } else if (idet==static_cast<int>(HcalForward)) {
00220     etaMn = etaMin[2]; etaMx = etaMax[2];
00221   }
00222   if (corr) {
00223     if (etaR >= nOff[2] && depth == 3 && idet == static_cast<int>(HcalBarrel))
00224       idet = static_cast<int>(HcalEndcap);
00225   }
00226   double eta = 0, deta = 0, phi = 0, dphi = 0, rz = 0, drz = 0;
00227   bool   ok = false, flagrz = true;
00228   if ((idet==static_cast<int>(HcalBarrel)||idet==static_cast<int>(HcalEndcap)||
00229        idet==static_cast<int>(HcalOuter)||idet==static_cast<int>(HcalForward))
00230       && etaR >=etaMn && etaR <= etaMx)
00231     ok = true;
00232   if (idet == static_cast<int>(HcalEndcap)) {
00233     if      (depth < 3 && etaR <= etaMin[1]) ok = false;
00234     else if (depth > 2 && etaR == nOff[1])   ok = false;
00235   }
00236   if (ok) {
00237     int maxlay = (int)(rHB.size());
00238     if (idet == static_cast<int>(HcalEndcap)) maxlay = (int)(zHE.size());
00239     eta  = getEta(idet, etaR, zside, depth);
00240     deta = deltaEta(idet, etaR, depth);
00241     double fibin, fioff;
00242     if      (idet == static_cast<int>(HcalBarrel)||
00243              idet == static_cast<int>(HcalOuter)) {
00244       fioff = phioff[0];
00245       fibin = phibin[etaR-1];
00246     } else if (idet == static_cast<int>(HcalEndcap)) {
00247       fioff = phioff[1];
00248       fibin = phibin[etaR-1];
00249     } else {
00250       fioff = phioff[2];
00251       fibin = phibin[nEta+etaR-etaMin[2]-1];
00252       if  (etaR > etaMax[2]-2 ) fioff += 0.5*fibin; 
00253     }
00254     phi  = fioff + (iphi - 0.5)*fibin;
00255     dphi = 0.5*fibin;
00256     if (idet == static_cast<int>(HcalForward)) {
00257       int ir = nR + etaMin[2] - etaR - 1;
00258       if (ir > 0 && ir < nR) {
00259         rz     = 0.5*(rTable[ir]+rTable[ir-1]);
00260         drz    = 0.5*(rTable[ir]-rTable[ir-1]);
00261       } else {
00262         ok     = false;
00263 #ifdef DebugLog
00264         LogDebug("HCalGeom") << "HcalNumberingFromDDD: wrong eta " << etaR 
00265                              << " ("  << ir << "/" << nR << ") Detector "
00266                              << idet;
00267 #endif
00268       }
00269       if (depth != 1 && depth != 2) {
00270         ok     = false;
00271 #ifdef DebugLog
00272         LogDebug("HCalGeom") << "HcalNumberingFromDDD: wrong depth " << depth
00273                              << " in Detector " << idet;
00274 #endif
00275       }
00276     } else if (etaR <= nEta) {
00277       int depth0 = depth1[etaR-1];
00278       int kphi   = iphi + int((phioff[3]+0.1)/fibin);
00279       kphi       = (kphi-1)%4 + 1;
00280       if (etaR == nOff[0] && (kphi == 2 || kphi == 3)) depth0--;
00281       int laymin, laymax;
00282       if (depth == 1) {
00283         laymin = 1;
00284         if (idet==static_cast<int>(HcalEndcap)) laymin = 2;
00285         laymax = depth0;
00286         if (nOff.size() > 12) {
00287           if (etaR == nOff[6]) {
00288             laymin = nOff[7];
00289             laymax = nOff[8];
00290           } else if (etaR == nOff[9]) {
00291             laymin = nOff[10];
00292           }
00293         }
00294       } else if (depth == 2) {
00295         laymin = depth0+1;
00296         laymax = depth2[etaR-1];
00297         if (etaR==etaMax[0] && idet==static_cast<int>(HcalBarrel) &&
00298             nOff.size()>3) laymax = nOff[3];
00299         if (nOff.size() > 12) {
00300           if (etaR == nOff[9]) laymax = nOff[11];
00301           if (etaR == nOff[6]) laymax = nOff[12];
00302         }
00303       } else  if (depth == 3) {
00304         laymin = depth2[etaR-1]+1;
00305         laymax = depth3[etaR-1];
00306         if (etaR<=etaMin[1] && idet==static_cast<int>(HcalEndcap)) {
00307           if (nOff.size() > 4) laymin = nOff[4];
00308           if (nOff.size() > 5) laymax = nOff[5];
00309         }
00310       } else {
00311         laymin = depth3[etaR-1]+1;
00312         laymax = maxlay;
00313       }
00314       if (idet == static_cast<int>(HcalOuter) && nOff.size() > 13) {
00315         if (etaR > nOff[13] && laymin <= laymax) laymin = laymax;
00316       }
00317       double d1=0, d2=0;
00318       if (laymin <= maxlay && laymax <= maxlay && laymin <= laymax) {
00319         if (idet == static_cast<int>(HcalEndcap)) {
00320           flagrz = false;
00321           if (depth == 1 || laymin <= 1) d1 = zHE[laymin-1] - dzHE[laymin-1];
00322           else                           d1 = zHE[laymin-2] + dzHE[laymin-2];
00323           d2     = zHE[laymax-1] + dzHE[laymax-1];
00324         } else {
00325           if (idet == static_cast<int>(HcalOuter) ||
00326               depth == 1 || laymin <=1) d1 = rHB[laymin-1] - drHB[laymin-1];
00327           else                          d1 = rHB[laymin-2] + drHB[laymin-1];
00328           d2     = rHB[laymax-1] + drHB[laymax-1];
00329         }
00330         rz     = 0.5*(d2+d1);
00331         drz    = 0.5*(d2-d1);
00332       } else {
00333         ok = false;
00334 #ifdef DebugLog
00335         LogDebug("HCalGeom") << "HcalNumberingFromDDD: wrong depth " << depth
00336                              << " (Layer minimum " << laymin << " maximum " 
00337                              << laymax << " maxLay " << maxlay << ")";
00338 #endif
00339       }
00340     } else {
00341       ok = false;
00342 #ifdef DebugLog
00343       LogDebug("HCalGeom") << "HcalNumberingFromDDD: wrong eta " << etaR
00344                            << "/" << nEta << " Detector " << idet;
00345 #endif
00346     }
00347   } else {
00348     ok = false;
00349 #ifdef DebugLog
00350     LogDebug("HCalGeom") << "HcalNumberingFromDDD: wrong eta " << etaR 
00351                          << " det " << idet;
00352 #endif
00353   }
00354   HcalCellType::HcalCell tmp(ok,eta,deta,phi,dphi,rz,drz,flagrz);
00355 
00356 #ifdef DebugLog
00357   LogDebug("HCalGeom") << "HcalNumberingFromDDD: det/side/depth/etaR/phi "
00358                        << det  << "/" << zside << "/" << depth << "/" << etaR
00359                        << "/" << iphi << " Cell Flag " << tmp.ok << " " 
00360                        << tmp.eta << " " << tmp.deta << " phi " << tmp.phi 
00361                        << " " << tmp.dphi << " r(z) " << tmp.rz  << " " 
00362                        << tmp.drz << " " << tmp.flagrz;
00363 #endif
00364   return tmp;
00365 }
00366 
00367 std::vector<double> HcalNumberingFromDDD::getEtaTable() const {
00368 
00369   std::vector<double> tmp = etaTable;
00370   return tmp;
00371 }
00372 
00373 unsigned int HcalNumberingFromDDD::numberOfCells(HcalSubdetector subdet) const{
00374 
00375   unsigned int num = 0;
00376   std::vector<HcalCellType> cellTypes = HcalCellTypes(subdet);
00377   for (unsigned int i=0; i<cellTypes.size(); i++) {
00378     num += (unsigned int)(cellTypes[i].nPhiBins());
00379     if (cellTypes[i].nHalves() > 1) 
00380       num += (unsigned int)(cellTypes[i].nPhiBins());
00381     num -= (unsigned int)(cellTypes[i].nPhiMissingBins());
00382   }
00383 #ifdef DebugLog
00384   LogDebug ("HCalGeom") << "HcalNumberingFromDDD:numberOfCells " 
00385                         << cellTypes.size()  << " " << num 
00386                         << " for subdetector " << subdet;
00387 #endif
00388   return num;
00389 }
00390 
00391 std::vector<HcalCellType> HcalNumberingFromDDD::HcalCellTypes() const{
00392 
00393   std::vector<HcalCellType> cellTypes =HcalCellTypes(HcalBarrel);
00394 #ifdef DebugLog
00395   LogDebug ("HCalGeom") << "HcalNumberingFromDDD: " << cellTypes.size()
00396                         << " cells of type HCal Barrel";
00397   for (unsigned int i=0; i<cellTypes.size(); i++)
00398     LogDebug ("HCalGeom") << "Cell " << i << " " << cellTypes[i];
00399 #endif
00400 
00401   std::vector<HcalCellType> hoCells   =HcalCellTypes(HcalOuter);
00402 #ifdef DebugLog
00403   LogDebug ("HCalGeom") << "HcalNumberingFromDDD: " << hoCells.size()
00404                         << " cells of type HCal Outer";
00405   for (unsigned int i=0; i<hoCells.size(); i++)
00406     LogDebug ("HCalGeom") << "Cell " << i << " " << hoCells[i];
00407 #endif
00408   cellTypes.insert(cellTypes.end(), hoCells.begin(), hoCells.end());
00409 
00410   std::vector<HcalCellType> heCells   =HcalCellTypes(HcalEndcap);
00411 #ifdef DebugLog
00412   LogDebug ("HCalGeom") << "HcalNumberingFromDDD: " << heCells.size()
00413                         << " cells of type HCal Endcap";
00414   for (unsigned int i=0; i<heCells.size(); i++)
00415     LogDebug ("HCalGeom") << "Cell " << i << " " << heCells[i];
00416 #endif
00417   cellTypes.insert(cellTypes.end(), heCells.begin(), heCells.end());
00418 
00419   std::vector<HcalCellType> hfCells   =HcalCellTypes(HcalForward);
00420 #ifdef DebugLog
00421   LogDebug ("HCalGeom") << "HcalNumberingFromDDD: " << hfCells.size()
00422                         << " cells of type HCal Forward";
00423   for (unsigned int i=0; i<hfCells.size(); i++)
00424     LogDebug ("HCalGeom") << "Cell " << i << " " << hfCells[i];
00425 #endif
00426   cellTypes.insert(cellTypes.end(), hfCells.begin(), hfCells.end());
00427 
00428   return cellTypes;
00429 }
00430 
00431 std::vector<HcalCellType> HcalNumberingFromDDD::HcalCellTypes(HcalSubdetector subdet) const {
00432 
00433   std::vector<HcalCellType> cellTypes;
00434   if (subdet == HcalForward) {
00435     if (dzVcal < 0) return cellTypes;
00436   }
00437 
00438   int    dmin, dmax, indx, nz, nmod;
00439   double hsize = 0;
00440   switch(subdet) {
00441   case HcalEndcap:
00442     dmin = 1; dmax = 3; indx = 1; nz = nzHE; nmod = nmodHE;
00443     break;
00444   case HcalForward:
00445     dmin = 1; dmax = 2; indx = 2; nz = 2; nmod = 18; 
00446     break;
00447   case HcalOuter:
00448     dmin = 4; dmax = 4; indx = 0; nz = nzHB; nmod = nmodHB;
00449     break;
00450   default:
00451     dmin = 1; dmax = 3; indx = 0; nz = nzHB; nmod = nmodHB;
00452     break;
00453   }
00454 
00455   int phi = 1, zside  = 1;
00456   bool cor = false;
00457 
00458   // Get the Cells 
00459   int subdet0 = static_cast<int>(subdet);
00460   for (int depth=dmin; depth<=dmax; depth++) {
00461     int    shift = getShift(subdet, depth);
00462     double gain  = getGain (subdet, depth);
00463     if (subdet == HcalForward) {
00464       if (depth == 1) hsize = dzVcal;
00465       else            hsize = dzVcal-0.5*dlShort;
00466     }
00467     for (int eta=etaMin[indx]; eta<= etaMax[indx]; eta++) {
00468       HcalCellType::HcalCell temp1 = cell(subdet0,zside,depth,eta,phi,cor);
00469       if (temp1.ok) {
00470         int units = unitPhi (subdet0, eta);
00471         HcalCellType temp2(subdet, eta, phi, depth, temp1,
00472                                          shift, gain, nz, nmod, hsize, units);
00473         if (subdet == HcalOuter && nOff.size() > 17) {
00474           if (eta == nOff[15]) {
00475             std::vector<int> missPlus, missMinus;
00476             int kk = 18;
00477             for (int miss=0; miss<nOff[16]; miss++) {
00478               missPlus.push_back(nOff[kk]);
00479               kk++;
00480             }
00481             for (int miss=0; miss<nOff[17]; miss++) {
00482               missMinus.push_back(nOff[kk]);
00483               kk++;
00484             }
00485             temp2.setMissingPhi(missPlus, missMinus);
00486           }
00487         }
00488         cellTypes.push_back(temp2);
00489       }
00490     }
00491   }
00492   return cellTypes;
00493 }
00494 
00495 void HcalNumberingFromDDD::printTile() {
00496  
00497   std::cout << "Tile Information for HB:\n" << "========================\n\n";
00498   for (int eta=etaMin[0]; eta<= etaMax[0]; eta++) {
00499     int dmax = 1;
00500     if (depth1[eta-1] < 17) dmax = 2;
00501     for (int depth=1; depth<=dmax; depth++) 
00502       tileHB(eta, depth);
00503   }
00504 
00505   std::cout << "\nTile Information for HE:\n" <<"========================\n\n";
00506   for (int eta=etaMin[1]; eta<= etaMax[1]; eta++) {
00507     int dmin=1, dmax=3;
00508     if (eta == etaMin[1]) {
00509       dmin = 3;
00510     } else if (depth1[eta-1] > 18) {
00511       dmax = 1;
00512     } else if (depth2[eta-1] > 18) {
00513       dmax = 2;
00514     }
00515     for (int depth=dmin; depth<=dmax; depth++)
00516       tileHE(eta, depth);
00517   }
00518 }
00519 
00520 double HcalNumberingFromDDD::getEta(int det, int etaR, int zside,
00521                                     int depth) const {
00522 
00523   double tmp = 0;
00524   if (det == static_cast<int>(HcalForward)) {
00525     int ir = nR + etaMin[2] - etaR - 1;
00526     if (ir > 0 && ir < nR) {
00527       double z = zVcal;
00528       if (depth != 1) z += dlShort;
00529       tmp = 0.5*(getEta(rTable[ir-1],z)+getEta(rTable[ir],z));
00530     }
00531   } else {
00532     if (etaR > 0 && etaR < nEta) {
00533       if (etaR == nOff[1]-1 && depth > 2) {
00534         tmp = 0.5*(etaTable[etaR+1]+etaTable[etaR-1]);
00535       } else if (det == static_cast<int>(HcalOuter) && nOff.size() > 13) {
00536         if (etaR == nOff[13]) {
00537           tmp = 0.5*(etaHO[0]+etaTable[etaR-1]);
00538         } else if (etaR == nOff[13]+1) {
00539           tmp = 0.5*(etaTable[etaR]+etaHO[1]);
00540         } else if (etaR == nOff[14]) {
00541           tmp = 0.5*(etaHO[2]+etaTable[etaR-1]);
00542         } else if (etaR == nOff[14]+1) {
00543           tmp = 0.5*(etaTable[etaR]+etaHO[3]);
00544         } else {
00545           tmp = 0.5*(etaTable[etaR]+etaTable[etaR-1]);
00546         }
00547       } else {
00548         tmp = 0.5*(etaTable[etaR]+etaTable[etaR-1]);
00549       }
00550     }
00551   } 
00552   if (zside == 0) tmp = -tmp;
00553 #ifdef DebugLog
00554   LogDebug("HCalGeom") << "HcalNumberingFromDDD::getEta " << etaR << " " 
00555                        << zside << " " << depth << " ==> " << tmp;
00556 #endif
00557   return tmp;
00558 }
00559  
00560 double HcalNumberingFromDDD::getEta(double r, double z) const {
00561 
00562   double tmp = 0;
00563   if (z != 0) tmp = -log(tan(0.5*atan(r/z)));
00564 #ifdef DebugLog
00565   LogDebug("HCalGeom") << "HcalNumberingFromDDD::getEta " << r << " " << z 
00566                        << " ==> " << tmp;
00567 #endif
00568   return tmp;
00569 }
00570 
00571 double HcalNumberingFromDDD::deltaEta(int det, int etaR, int depth) const {
00572 
00573   double tmp = 0;
00574   if (det == static_cast<int>(HcalForward)) {
00575     int ir = nR + etaMin[2] - etaR - 1;
00576     if (ir > 0 && ir < nR) {
00577       double z = zVcal;
00578       if (depth != 1) z += dlShort;
00579       tmp = 0.5*(getEta(rTable[ir-1],z)-getEta(rTable[ir],z));
00580     }
00581   } else {
00582     if (etaR > 0 && etaR < nEta) {
00583       if (etaR == nOff[1]-1 && depth > 2) {
00584         tmp = 0.5*(etaTable[etaR+1]-etaTable[etaR-1]);
00585       } else if (det == static_cast<int>(HcalOuter) && nOff.size() > 13) {
00586         if (etaR == nOff[13]) {
00587           tmp = 0.5*(etaHO[0]-etaTable[etaR-1]);
00588         } else if (etaR == nOff[13]+1) {
00589           tmp = 0.5*(etaTable[etaR]-etaHO[1]);
00590         } else if (etaR == nOff[14]) {
00591           tmp = 0.5*(etaHO[2]-etaTable[etaR-1]);
00592         } else if (etaR == nOff[14]+1) {
00593           tmp = 0.5*(etaTable[etaR]-etaHO[3]);
00594         } else {
00595           tmp = 0.5*(etaTable[etaR]-etaTable[etaR-1]);
00596         }
00597       } else {
00598         tmp = 0.5*(etaTable[etaR]-etaTable[etaR-1]);
00599       }
00600     } 
00601   }
00602 #ifdef DebugLog
00603   LogDebug("HCalGeom") << "HcalNumberingFromDDD::deltaEta " << etaR << " " 
00604                        << depth << " ==> " << tmp;
00605 #endif
00606   return tmp;
00607 }
00608 
00609 void HcalNumberingFromDDD::initialize(std::string & name, 
00610                                       const DDCompactView & cpv) {
00611 
00612   std::string attribute = "ReadOutName";
00613   edm::LogInfo("HCalGeom") << "HcalNumberingFromDDD: Initailise for " << name 
00614                            << " as " << attribute;
00615 
00616   DDSpecificsFilter filter;
00617   DDValue           ddv(attribute,name,0);
00618   filter.setCriteria(ddv,DDSpecificsFilter::equals);
00619   DDFilteredView fv(cpv);
00620   fv.addFilter(filter);
00621   bool ok = fv.firstChild();
00622 
00623   if (ok) {
00624     //Load the SpecPars
00625     loadSpecPars(fv);
00626 
00627     //Load the Geometry parameters
00628     loadGeometry(fv);
00629   } else {
00630     edm::LogError("HCalGeom") << "HcalNumberingFromDDD: cannot get filtered "
00631                               << " view for " << attribute << " matching "
00632                               << name;
00633     throw DDException("HcalNumberingFromDDD: cannot match "+attribute+" to "+name);
00634   }
00635 
00636 #ifdef DebugLog
00637   std::vector<HcalCellType> cellTypes = HcalCellTypes();
00638   LogDebug ("HCalGeom") << "HcalNumberingFromDDD: " << cellTypes.size()
00639                         << " cells of type HCal (All)";
00640   for (unsigned int i=0; i<cellTypes.size(); i++)
00641     LogDebug ("HCalGeom") << "Cell " << i << " " << cellTypes[i];
00642 #endif
00643 }
00644 
00645 void HcalNumberingFromDDD::loadSpecPars(DDFilteredView fv) {
00646 
00647   DDsvalues_type sv(fv.mergedSpecifics());
00648 
00649   // Phi Offset
00650   int i, nphi=4;
00651   std::vector<double> tmp1 = getDDDArray("phioff",sv,nphi);
00652   phioff.resize(tmp1.size());
00653   for (i=0; i<nphi; i++) {
00654     phioff[i] = tmp1[i];
00655 #ifdef DebugLog
00656     LogDebug("HCalGeom") << "HcalNumberingFromDDD: phioff[" << i << "] = "
00657                          << phioff[i]/CLHEP::deg;
00658 #endif
00659   }
00660 
00661   //Eta table
00662   nEta     = -1;
00663   std::vector<double> tmp2 = getDDDArray("etaTable",sv,nEta);
00664   etaTable.resize(tmp2.size());
00665   for (i=0; i<nEta; i++) {
00666     etaTable[i] = tmp2[i];
00667 #ifdef DebugLog
00668     LogDebug("HCalGeom") << "HcalNumberingFromDDD: etaTable[" << i << "] = "
00669                          << etaTable[i];
00670 #endif
00671   }
00672 
00673   //R table
00674   nR     = -1;
00675   std::vector<double> tmp3 = getDDDArray("rTable",sv,nR);
00676   rTable.resize(tmp3.size());
00677   for (i=0; i<nR; i++) {
00678     rTable[i] = tmp3[i];
00679 #ifdef DebugLog
00680     LogDebug("HCalGeom") << "HcalNumberingFromDDD: rTable[" << i << "] = "
00681                          << rTable[i]/CLHEP::cm;
00682 #endif
00683   }
00684 
00685   //Phi bins
00686   nPhi   = nEta + nR - 2;
00687   std::vector<double> tmp4 = getDDDArray("phibin",sv,nPhi);
00688   phibin.resize(tmp4.size());
00689   for (i=0; i<nPhi; i++) {
00690     phibin[i] = tmp4[i];
00691 #ifdef DebugLog
00692     LogDebug("HCalGeom") << "HcalNumberingFromDDD: phibin[" << i << "] = "
00693                          << phibin[i]/CLHEP::deg;
00694 #endif
00695   }
00696 
00697   //Layer boundaries for depths 1, 2, 3, 4
00698   nDepth            = nEta - 1;
00699   std::vector<double> d1 = getDDDArray("depth1",sv,nDepth);
00700   nDepth            = nEta - 1;
00701   std::vector<double> d2 = getDDDArray("depth2",sv,nDepth);
00702   nDepth            = nEta - 1;
00703   std::vector<double> d3 = getDDDArray("depth3",sv,nDepth);
00704 #ifdef DebugLog
00705   LogDebug("HCalGeom") << "HcalNumberingFromDDD: " << nDepth << " Depths";
00706 #endif
00707   depth1.resize(nDepth);
00708   depth2.resize(nDepth);
00709   depth3.resize(nDepth);
00710   for (i=0; i<nDepth; i++) {
00711     depth1[i] = static_cast<int>(d1[i]);
00712     depth2[i] = static_cast<int>(d2[i]);
00713     depth3[i] = static_cast<int>(d3[i]);
00714 #ifdef DebugLog
00715     LogDebug("HCalGeom") << "HcalNumberingFromDDD: depth1[" << i << "] = " 
00716                          << depth1[i] << " depth2[" << i << "]  = "<< depth2[i]
00717                          << " depth3[" << i << "] = " << depth3[i];
00718 #endif
00719   }
00720 
00721   // Minimum and maximum eta boundaries
00722   int ndx  = 3;
00723   std::vector<double>  tmp5 = getDDDArray("etaMin",sv,ndx);
00724   std::vector<double>  tmp6 = getDDDArray("etaMax",sv,ndx);
00725   etaMin.resize(ndx);
00726   etaMax.resize(ndx);
00727   for (i=0; i<ndx; i++) {
00728     etaMin[i] = static_cast<int>(tmp5[i]);
00729     etaMax[i] = static_cast<int>(tmp6[i]);
00730   }
00731   etaMin[0] = 1;
00732   etaMax[1] = nEta-1;
00733   etaMax[2] = etaMin[2]+nR-2;
00734 #ifdef DebugLog
00735   for (i=0; i<ndx; i++) 
00736     LogDebug("HCalGeom") << "HcalNumberingFromDDD: etaMin[" << i << "] = "
00737                          << etaMin[i] << " etaMax[" << i << "] = "<< etaMax[i];
00738 #endif
00739 
00740   // Geometry parameters for HF
00741   int ngpar = 7;
00742   std::vector<double> gpar = getDDDArray("gparHF",sv,ngpar);
00743   dlShort = gpar[0];
00744   zVcal   = gpar[4];
00745 #ifdef DebugLog
00746   LogDebug("HCalGeom") << "HcalNumberingFromDDD: dlShort " << dlShort
00747                        << " zVcal " << zVcal;
00748 #endif
00749 
00750   // nOff
00751   int noff = 3;
00752   std::vector<double>  nvec = getDDDArray("noff",sv,noff);
00753   nOff.resize(noff);
00754   for (i=0; i<noff; i++) {
00755     nOff[i] = static_cast<int>(nvec[i]);
00756 #ifdef DebugLog
00757     LogDebug("HCalGeom") << "HcalNumberingFromDDD: nOff[" << i << "] = " 
00758                          << nOff[i];
00759 #endif
00760   }
00761 
00762   //Gains and Shifts for HB depths
00763   ndx                  = 4;
00764   gainHB               = getDDDArray("HBGains",sv,ndx);
00765   std::vector<double>  tmp7 = getDDDArray("HBShift",sv,ndx);
00766   shiftHB.resize(ndx);
00767 #ifdef DebugLog
00768   LogDebug("HCalGeom") << "HcalNumberingFromDDD:: Gain factor and Shift for "
00769                        << "HB depth layers:";
00770 #endif
00771   for (i=0; i<ndx; i++) {
00772     shiftHB[i] = static_cast<int>(tmp7[i]);
00773 #ifdef DebugLog
00774     LogDebug("HCalGeom") <<"HcalNumberingFromDDD:: gainHB[" <<  i << "] = " 
00775                          << gainHB[i] << " shiftHB[" << i << "] = " 
00776                          << shiftHB[i];
00777 #endif
00778   }
00779 
00780   //Gains and Shifts for HB depths
00781   ndx                  = 4;
00782   gainHE               = getDDDArray("HEGains",sv,ndx);
00783   std::vector<double>  tmp8 = getDDDArray("HEShift",sv,ndx);
00784   shiftHE.resize(ndx);
00785 #ifdef DebugLog
00786   LogDebug("HCalGeom") << "HcalNumberingFromDDD:: Gain factor and Shift for "
00787                        << "HE depth layers:";
00788 #endif
00789   for (i=0; i<ndx; i++) {
00790     shiftHE[i] = static_cast<int>(tmp8[i]);
00791 #ifdef DebugLog
00792     LogDebug("HCalGeom") <<"HcalNumberingFromDDD:: gainHE[" <<  i << "] = " 
00793                          << gainHE[i] << " shiftHE[" << i << "] = " 
00794                          << shiftHE[i];
00795 #endif
00796   }
00797 
00798   //Gains and Shifts for HF depths
00799   ndx                  = 4;
00800   gainHF               = getDDDArray("HFGains",sv,ndx);
00801   std::vector<double>  tmp9 = getDDDArray("HFShift",sv,ndx);
00802   shiftHF.resize(ndx);
00803 #ifdef DebugLog
00804   LogDebug("HCalGeom") << "HcalNumberingFromDDD:: Gain factor and Shift for "
00805                        << "HF depth layers:";
00806 #endif
00807   for (i=0; i<ndx; i++) {
00808     shiftHF[i] = static_cast<int>(tmp9[i]);
00809 #ifdef DebugLog
00810     LogDebug("HCalGeom") <<"HcalNumberingFromDDD:: gainHF[" <<  i << "] = " 
00811                          << gainHF[i] << " shiftHF[" << i << "] = " 
00812                          << shiftHF[i];
00813 #endif
00814   }
00815 }
00816 
00817 void HcalNumberingFromDDD::loadGeometry(DDFilteredView fv) {
00818 
00819   bool dodet=true, hf=false;
00820   std::vector<double> rb(20,0.0), ze(20,0.0), thkb(20,-1.0), thke(20,-1.0);
00821   std::vector<double> zho;
00822   std::vector<int>    ib(20,0),   ie(20,0);
00823   std::vector<int>    izb, phib, ize, phie, izf, phif;
00824   std::vector<double> rxb;
00825   rhoxb.clear(); zxb.clear(); dyxb.clear(); dzxb.clear();
00826   layb.clear(); laye.clear();
00827   zxe.clear(); rhoxe.clear(); dyxe.clear(); dx1e.clear(); dx2e.clear();
00828   double zf = 0;
00829   dzVcal = -1.;
00830 
00831   while (dodet) {
00832     DDTranslation    t    = fv.translation();
00833     std::vector<int> copy = fv.copyNumbers();
00834     const DDSolid & sol  = fv.logicalPart().solid();
00835     int idet = 0, lay = -1;
00836     int nsiz = (int)(copy.size());
00837     if (nsiz>0) lay  = copy[nsiz-1]/10;
00838     if (nsiz>1) idet = copy[nsiz-2]/1000;
00839     double dx=0, dy=0, dz=0, dx1=0, dx2=0;
00840     if (sol.shape() == 1) {
00841       const DDBox & box = static_cast<DDBox>(fv.logicalPart().solid());
00842       dx = box.halfX();
00843       dy = box.halfY();
00844       dz = box.halfZ();
00845     } else if (sol.shape() == 3) {
00846       const DDTrap & trp = static_cast<DDTrap>(fv.logicalPart().solid());
00847       dx1= trp.x1();
00848       dx2= trp.x2();
00849       dx = 0.25*(trp.x1()+trp.x2()+trp.x3()+trp.x4());
00850       dy = 0.5*(trp.y1()+trp.y2());
00851       dz = trp.halfZ();
00852     } else if (sol.shape() == 2) {
00853       const DDTubs & tub = static_cast<DDTubs>(fv.logicalPart().solid());
00854       dx = tub.rIn();
00855       dy = tub.rOut();
00856       dz = tub.zhalf();
00857     }
00858     if (idet == 3) {
00859       // HB
00860 #ifdef DebugLog
00861       LogDebug("HCalGeom") << "HB " << sol.name() << " Shape " << sol.shape()
00862                            << " Layer " << lay << " R " << t.Rho();
00863 #endif
00864       if (lay >=0 && lay < 20) {
00865         ib[lay]++;
00866         rb[lay] += t.Rho();
00867         if (thkb[lay] <= 0) {
00868           if (lay < 17) thkb[lay] = dx;
00869           else          thkb[lay] = std::min(dx,dy);
00870         }
00871         if (lay < 17) {
00872           bool found = false;
00873           for (unsigned int k=0; k<rxb.size(); k++) {
00874             if (std::abs(rxb[k]-t.Rho()) < 0.01) {
00875               found = true;
00876               break;
00877             }
00878           }
00879           if (!found) {
00880             rxb.push_back(t.Rho());
00881             rhoxb.push_back(t.Rho()*std::cos(t.phi()));
00882             zxb.push_back(std::abs(t.z()));
00883             dyxb.push_back(2.*dy);
00884             dzxb.push_back(2.*dz);
00885             layb.push_back(lay);
00886           }
00887         }
00888       }
00889       if (lay == 2) {
00890         int iz = copy[nsiz-5];
00891         int fi = copy[nsiz-4];
00892         unsigned int it1 = find(iz, izb);
00893         if (it1 == izb.size())  izb.push_back(iz);
00894         unsigned int it2 = find(fi, phib);
00895         if (it2 == phib.size()) phib.push_back(fi);
00896       }
00897       if (lay == 18) {
00898         int ifi=-1, ich=-1;
00899         if (nsiz>2) ifi = copy[nsiz-3];
00900         if (nsiz>3) ich = copy[nsiz-4];
00901         double z1 = std::abs((t.z()) + dz);
00902         double z2 = std::abs((t.z()) - dz);
00903         if (std::abs(z1-z2) < 0.01) z1 = 0;
00904         if (ifi == 1 && ich == 4) {
00905           if (z1 > z2) {
00906             double tmp = z1;
00907             z1 = z2;
00908             z2 = tmp;
00909           }
00910           bool sok = true;
00911           for (unsigned int kk=0; kk<zho.size(); kk++) {
00912             if (std::abs(z2-zho[kk]) < 0.01) {
00913               sok = false;
00914               break;
00915             }   else if (z2 < zho[kk]) {
00916               zho.resize(zho.size()+2);
00917               for (unsigned int kz=zho.size()-1; kz>kk+1; kz=kz-2) {
00918                 zho[kz]   = zho[kz-2];
00919                 zho[kz-1] = zho[kz-3];
00920               }
00921               zho[kk+1] = z2;
00922               zho[kk]   = z1;
00923               sok = false;
00924               break;
00925             }
00926           }
00927           if (sok) {
00928             zho.push_back(z1);
00929             zho.push_back(z2);
00930           }
00931           LogDebug("HCalGeom") << "Detector " << idet << " Lay " << lay << " fi " << ifi << " " << ich << " z " << z1 << " " << z2;
00932         }
00933       }
00934     } else if (idet == 4) {
00935       // HE
00936 #ifdef DebugLog
00937       LogDebug("HCalGeom") << "HE " << sol.name() << " Shape " << sol.shape()
00938                            << " Layer " << lay << " Z " << t.z();
00939 #endif
00940       if (lay >=0 && lay < 20) {
00941         ie[lay]++;
00942         ze[lay] += std::abs(t.z());
00943         if (thke[lay] <= 0) thke[lay] = dz;
00944         bool found = false;
00945         for (unsigned int k=0; k<zxe.size(); k++) {
00946           if (std::abs(zxe[k]-std::abs(t.z())) < 0.01) {
00947             found = true;
00948             break;
00949           }
00950         }
00951         if (!found) {
00952           zxe.push_back(std::abs(t.z()));
00953           rhoxe.push_back(t.Rho()*std::cos(t.phi()));
00954           dyxe.push_back(dy*std::cos(t.phi()));
00955           dx1 -= 0.5*(t.rho()-dy)*std::cos(t.phi())*std::tan(10*CLHEP::deg);
00956           dx2 -= 0.5*(t.rho()+dy)*std::cos(t.phi())*std::tan(10*CLHEP::deg);
00957           dx1e.push_back(-dx1);
00958           dx2e.push_back(-dx2);
00959           laye.push_back(lay);
00960         }
00961       }
00962       if (copy[nsiz-1] == 21) {
00963         int iz = copy[nsiz-7];
00964         int fi = copy[nsiz-5];
00965         unsigned int it1 = find(iz, ize);
00966         if (it1 == ize.size())  ize.push_back(iz);
00967         unsigned int it2 = find(fi, phie);
00968         if (it2 == phie.size()) phie.push_back(fi);
00969       }
00970     } else if (idet == 5) {
00971       // HF
00972       if (!hf) {
00973         const std::vector<double> & paras = sol.parameters();
00974 #ifdef DebugLog
00975         LogDebug("HCalGeom") << "HF " << sol.name() << " Shape " << sol.shape()
00976                              << " Z " << t.z() << " with " << paras.size()
00977                              << " Parameters";
00978         for (unsigned j=0; j<paras.size(); j++)
00979           LogDebug("HCalGeom") << "HF Parameter[" << j << "] = " << paras[j];
00980 #endif
00981         zf  = fabs(t.z());
00982         if (sol.shape() == ddpolycone_rrz) {
00983           int nz  = (int)(paras.size())-3;
00984           zf     += paras[3];
00985           dzVcal  = 0.5*(paras[nz]-paras[3]);
00986           hf      = true;
00987         } else if (sol.shape() == ddtubs || sol.shape() == ddcons) {
00988           dzVcal  = paras[0];
00989           zf     -= paras[0];
00990           hf      = true;
00991         }
00992       }
00993 #ifdef DebugLog
00994     } else {
00995       LogDebug("HCalGeom") << "Unknown Detector " << idet << " for " 
00996                            << sol.name() << " Shape " << sol.shape() << " R " 
00997                            << t.Rho() << " Z " << t.z();
00998 #endif
00999     }
01000     dodet = fv.next();
01001   }
01002 
01003   int ibmx = 0, iemx = 0;
01004   for (int i = 0; i < 20; i++) {
01005     if (ib[i]>0) {
01006       rb[i] /= (double)(ib[i]);
01007       ibmx   = i+1;
01008     }
01009     if (ie[i]>0) {
01010       ze[i] /= (double)(ie[i]);
01011       iemx   = i+1;
01012     }
01013 #ifdef DebugLog
01014     LogDebug("HCalGeom") << "Index " << i << " Barrel " << ib[i] << " "
01015                          << rb[i] << " Endcap " << ie[i] << " " << ze[i];
01016 #endif
01017   }
01018   for (int i = 4; i >= 0; i--) {
01019     if (ib[i] == 0) {rb[i] = rb[i+1]; thkb[i] = thkb[i+1];}
01020     if (ie[i] == 0) {ze[i] = ze[i+1]; thke[i] = thke[i+1];}
01021 #ifdef DebugLog
01022     if (ib[i] == 0 || ie[i] == 0)
01023       LogDebug("HCalGeom") << "Index " << i << " Barrel " << ib[i] << " "
01024                            << rb[i] << " Endcap " << ie[i] << " " << ze[i];
01025 #endif
01026   }
01027 
01028 #ifdef DebugLog
01029   for (unsigned int k=0; k<layb.size(); ++k)
01030     std::cout << "HB: " << layb[k] << " R " << rxb[k] << " " << rhoxb[k] << " Z " << zxb[k] << " DY " << dyxb[k] << " DZ " << dzxb[k] << "\n";
01031   for (unsigned int k=0; k<laye.size(); ++k) 
01032     std::cout << "HE: " << laye[k] << " R " << rhoxe[k] << " Z " << zxe[k] << " X1|X2 " << dx1e[k] << "|" << dx2e[k] << " DY " << dyxe[k] << "\n";
01033 
01034   printTile();
01035   LogDebug("HCalGeom") << "HcalNumberingFromDDD: Maximum Layer for HB " 
01036                        << ibmx << " for HE " << iemx << " Z for HF " << zf 
01037                        << " extent " << dzVcal;
01038 #endif
01039 
01040   if (ibmx > 0) {
01041     rHB.resize(ibmx);
01042     drHB.resize(ibmx);
01043     for (int i=0; i<ibmx; i++) {
01044       rHB[i]  = rb[i];
01045       drHB[i] = thkb[i];
01046 #ifdef DebugLog
01047       LogDebug("HCalGeom") << "HcalNumberingFromDDD: rHB[" << i << "] = "
01048                            << rHB[i] << " drHB[" << i << "] = " << drHB[i];
01049 #endif
01050     }
01051   }
01052   if (iemx > 0) {
01053     zHE.resize(iemx);
01054     dzHE.resize(iemx);
01055     for (int i=0; i<iemx; i++) {
01056       zHE[i]  = ze[i];
01057       dzHE[i] = thke[i];
01058 #ifdef DebugLog
01059       LogDebug("HCalGeom") << "HcalNumberingFromDDD: zHE[" << i << "] = "
01060                            << zHE[i] << " dzHE[" << i << "] = " << dzHE[i];
01061 #endif
01062     }
01063   }
01064 
01065   nzHB   = (int)(izb.size());
01066   nmodHB = (int)(phib.size());
01067 #ifdef DebugLog
01068   LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nzHB
01069                        << " barrel half-sectors";
01070   for (int i=0; i<nzHB; i++)
01071     LogDebug("HCalGeom") << "Section " << i << " Copy number " << izb[i];
01072   LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nmodHB
01073                        << " barrel modules";
01074   for (int i=0; i<nmodHB; i++)
01075     LogDebug("HCalGeom") << "Module " << i << " Copy number " << phib[i];
01076 #endif
01077 
01078   nzHE   = (int)(ize.size());
01079   nmodHE = (int)(phie.size());
01080 #ifdef DebugLog
01081   LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nzHE
01082                        << " endcap half-sectors";
01083   for (int i=0; i<nzHE; i++)
01084     LogDebug("HCalGeom") << "Section " << i << " Copy number " << ize[i];
01085   LogDebug("HCalGeom") << "HcalNumberingFromDDD::loadGeometry: " << nmodHE
01086                        << " endcap modules";
01087   for (int i=0; i<nmodHE; i++)
01088     LogDebug("HCalGeom") << "Module " << i << " Copy number " << phie[i];
01089 #endif
01090 
01091   LogDebug("HCalGeom") << "HO has Z of size " << zho.size();
01092   for (unsigned int kk=0; kk<zho.size(); kk++)
01093     LogDebug("HCalGeom") << "ZHO[" << kk << "] = " << zho[kk];
01094   if (ibmx > 17 && zho.size() > 2) {
01095     etaHO[0] = getEta(0.5*(rHB[17]+rHB[18]), zho[1]);
01096     etaHO[1] = getEta(rHB[18]+drHB[18], zho[2]);
01097     etaHO[2] = getEta(rHB[18]-drHB[18], zho[3]);
01098     etaHO[3] = getEta(rHB[18]+drHB[18], zho[4]);
01099   } else {
01100     etaHO[0] = etaTable[4];
01101     etaHO[1] = etaTable[4];
01102     etaHO[2] = etaTable[10];
01103     etaHO[3] = etaTable[10];
01104   }
01105   LogDebug("HCalGeom") << "HO Eta boundaries " << etaHO[0] << " " << etaHO[1]
01106                        << " " << etaHO[2] << " " << etaHO[3];
01107 }
01108 
01109 std::vector<double> HcalNumberingFromDDD::getDDDArray(const std::string & str, 
01110                                                       const DDsvalues_type & sv,
01111                                                       int & nmin) const {
01112 #ifdef DebugLog
01113   LogDebug("HCalGeom") << "HcalNumberingFromDDD:getDDDArray called for " 
01114                        << str << " with nMin "  << nmin;
01115 #endif
01116   DDValue value(str);
01117   if (DDfetch(&sv,value)) {
01118 #ifdef DebugLog
01119     LogDebug("HCalGeom") << "HcalNumberingFromDDD: " << value;
01120 #endif
01121     const std::vector<double> & fvec = value.doubles();
01122     int nval = fvec.size();
01123     if (nmin > 0) {
01124       if (nval < nmin) {
01125         edm::LogError("HCalGeom") << "HcalNumberingFromDDD : # of " << str 
01126                                   << " bins " << nval << " < " << nmin 
01127                                   << " ==> illegal";
01128         throw DDException("HcalNumberingFromDDD: cannot get array "+str);
01129       }
01130     } else {
01131       if (nval < 2) {
01132         edm::LogError("HCalGeom") << "HcalNumberingFromDDD : # of " << str
01133                                   << " bins " << nval << " < 2 ==> illegal"
01134                                   << " (nmin=" << nmin << ")";
01135         throw DDException("HcalNumberingFromDDD: cannot get array "+str);
01136       }
01137     }
01138     nmin = nval;
01139     return fvec;
01140   } else {
01141     edm::LogError("HCalGeom") << "HcalNumberingFromDDD: cannot get array "
01142                               << str;
01143     throw DDException("HcalNumberingFromDDD: cannot get array "+str);
01144   }
01145 }
01146 
01147 int HcalNumberingFromDDD::getShift(HcalSubdetector subdet, int depth) const {
01148 
01149   int shift;
01150   switch(subdet) {
01151   case HcalEndcap:
01152     shift = shiftHE[depth-1];
01153     break;
01154   case HcalForward:
01155     shift = shiftHF[depth-1];
01156     break;
01157   default:
01158     shift = shiftHB[depth-1];
01159     break;
01160   }
01161   return shift;
01162 }
01163 
01164 double HcalNumberingFromDDD::getGain(HcalSubdetector subdet, int depth) const {
01165 
01166   double gain;
01167   switch(subdet) {
01168   case HcalEndcap:
01169     gain = gainHE[depth-1];
01170     break;
01171   case HcalForward:
01172     gain = gainHF[depth-1];
01173     break;
01174   default:
01175     gain = gainHB[depth-1];
01176     break;
01177   }
01178   return gain;
01179 }
01180 
01181 unsigned int HcalNumberingFromDDD::find(int element, 
01182                                         std::vector<int> array) const {
01183 
01184   unsigned int id = array.size();
01185   for (unsigned int i = 0; i < array.size(); i++) {
01186     if (element == array[i]) {
01187       id = i;
01188       break;
01189     }
01190   }
01191   return id;
01192 }
01193 
01194 int HcalNumberingFromDDD::unitPhi(int det, int etaR) const {
01195 
01196   const double fiveDegInRad = 2*M_PI/72;
01197   int units=0;
01198   if (det == static_cast<int>(HcalForward))
01199     units=int(phibin[nEta+etaR-etaMin[2]-1]/fiveDegInRad+0.5);
01200   else 
01201     units=int(phibin[etaR-1]/fiveDegInRad+0.5);
01202 
01203   return units;
01204 }
01205 
01206 void HcalNumberingFromDDD::tileHB(int eta, int depth) {
01207 
01208   double etaL   = etaTable[eta-1];
01209   double thetaL = 2.*atan(exp(-etaL));
01210   double etaH   = etaTable[eta];
01211   double thetaH = 2.*atan(exp(-etaH));
01212   int    layL=0, layH=0;
01213   if (depth == 1) {
01214     layH = depth1[eta-1];
01215   } else {
01216     layL = depth1[eta-1];
01217     layH = depth2[eta-1];
01218   }
01219   std::cout << "\ntileHB:: eta|depth " << eta << "|" << depth << " theta " << thetaH/CLHEP::deg << ":" << thetaL/CLHEP::deg << " Layer " << layL << ":" << layH-1 << "\n";
01220   for (int lay=layL; lay<layH; ++lay) {
01221     std::vector<double> area(2,0);
01222     int kk=0;
01223     for (unsigned int k=0; k<layb.size(); ++k) {
01224       if (lay == layb[k]) {
01225         double zmin = rhoxb[k]*std::cos(thetaL)/std::sin(thetaL);
01226         double zmax = rhoxb[k]*std::cos(thetaH)/std::sin(thetaH);
01227         double dz   = (std::min(zmax,dzxb[k]) - zmin);
01228         if (dz > 0) {
01229           area[kk] = dz*dyxb[k];
01230           kk++;
01231         }
01232       }
01233     }
01234     if (area[0] > 0) std::cout << std::setw(2) << lay << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << "\n";
01235   }
01236 }
01237 
01238 void HcalNumberingFromDDD::tileHE(int eta, int depth) {
01239 
01240   double etaL   = etaTable[eta-1];
01241   double thetaL = 2.*atan(exp(-etaL));
01242   double etaH   = etaTable[eta];
01243   double thetaH = 2.*atan(exp(-etaH));
01244   int    layL=0, layH=0;
01245   if (eta == 16) {
01246     layH = depth3[eta-1];
01247   } else if (depth == 1) {
01248     layH = depth1[eta-1];
01249   } else if (depth == 2) {
01250     layL = depth1[eta-1];
01251     layH = depth2[eta-1];
01252   } else {
01253     layL = depth2[eta-1];
01254     layH = depth3[eta-1];
01255   }
01256   double phib  = phibin[eta-1];
01257   int nphi = 2;
01258   if (phib > 6*CLHEP::deg) nphi = 1;
01259   std::cout << "\ntileHE:: Eta/depth " << eta << "|" << depth << " theta " << thetaH/CLHEP::deg << ":" << thetaL/CLHEP::deg << " Layer " << layL << ":" << layH-1 << " phi " << nphi << "\n";
01260   for (int lay=layL; lay<layH; ++lay) {
01261     std::vector<double> area(4,0);
01262     int kk=0;
01263     for (unsigned int k=0; k<laye.size(); ++k) {
01264       if (lay == laye[k]) {
01265         double rmin = zxe[k]*std::tan(thetaH);
01266         double rmax = zxe[k]*std::tan(thetaL);
01267         if ((lay != 0 || eta == 18) && 
01268             (lay != 1 || (eta == 18 && rhoxe[k]-dyxe[k] > 1000) || (eta != 18 && rhoxe[k]-dyxe[k] < 1000)) &&
01269             rmin+30 < rhoxe[k]+dyxe[k] && rmax > rhoxe[k]-dyxe[k]) {
01270           rmin = std::max(rmin,rhoxe[k]-dyxe[k]);
01271           rmax = std::min(rmax,rhoxe[k]+dyxe[k]);
01272           double dx1 = rmin*std::tan(phib);
01273           double dx2 = rmax*std::tan(phib);
01274           double ar1=0, ar2=0;
01275           if (nphi == 1) {
01276             ar1 = 0.5*(rmax-rmin)*(dx1+dx2-4.*dx1e[k]);
01277           } else {
01278             ar1 = 0.5*(rmax-rmin)*(dx1+dx2-2.*dx1e[k]);
01279             ar2 = 0.5*(rmax-rmin)*((rmax+rmin)*tan(10.*CLHEP::deg)-4*dx1e[k])-ar1;
01280           }
01281           area[kk] = ar1;
01282           area[kk+2] = ar2;
01283           kk++;
01284         }
01285       }
01286     }
01287     if (area[0] > 0 && area[1] > 0) {
01288       int lay0 = lay-1;
01289       if (eta == 18) lay0++;
01290       if (nphi == 1) {
01291         std::cout << std::setw(2) << lay0 << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << "\n";
01292       } else {
01293         std::cout << std::setw(2) << lay0 << " Area " << std::setw(8) << area[0] << " " << std::setw(8) << area[1] << ":" << std::setw(8) << area[2] << " " << std::setw(8) << area[3] << "\n";
01294       }
01295     }
01296   }
01297 }
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/Geometry/HcalCommonData/src/HcalNumberingFromDDD.cc
