CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_1/src/Geometry/HcalTowerAlgo/src/HcalHardcodeGeometryLoader.cc

Go to the documentation of this file.
00001 #include "Geometry/HcalTowerAlgo/interface/HcalHardcodeGeometryLoader.h"
00002 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00003 #include "Geometry/CaloGeometry/interface/IdealObliquePrism.h"
00004 #include "Geometry/CaloGeometry/interface/IdealZPrism.h"
00005 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
00006 #include "HcalHardcodeGeometryData.h"
00007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00008 #include <algorithm>
00009 
00010 typedef CaloCellGeometry::CCGFloat CCGFloat ;
00011 
00012 HcalHardcodeGeometryLoader::HcalHardcodeGeometryLoader(const HcalTopology& ht) : 
00013    theTopology ( 0   ) ,
00014    extTopology ( &ht )
00015 {
00016    init();
00017 }
00018 
00019 void 
00020 HcalHardcodeGeometryLoader::init() 
00021 {
00022   theBarrelRadius = 190.;
00023   theHBThickness = 93.6; // just from drawings.  All thicknesses needs to be done right
00024 
00025   theHB15aThickness = theHBThickness * 12.0/17.0; // relative weight from layer count!
00026   theHB15bThickness = theHBThickness * 4.0/17.0;  // relative weight from layer count!
00027   theHB16aThickness = theHBThickness * 1.0/17.0;  // relative weight from layer count!
00028   theHB16bThickness = theHBThickness * 7.0/17.0;  // relative weight from layer count!
00029 
00030   theOuterRadius  = 406;
00031   theHOThickness = 1.;
00032 
00033   theHEZPos[0] = 388.0;
00034   theHEZPos[1] = 397.0;
00035   theHEZPos[2] = 450.0;
00036   theHEZPos[3] = 568.0;
00037 
00038   theHFZPos[0] = 1100.0;
00039   theHFZPos[1] = 1120.0;
00040   theHFThickness = 165;
00041 
00042 }
00043 
00044 
00045 HcalHardcodeGeometryLoader::ReturnType
00046 HcalHardcodeGeometryLoader::load( DetId::Detector det, 
00047                                   int             subdet ) 
00048 {
00049    HcalSubdetector hsub=static_cast<HcalSubdetector>( subdet );
00050    ReturnType hg( new HcalGeometry( *extTopology) );
00051 
00052    if( hg->cornersMgr() == 0 ) hg->allocateCorners( extTopology->ncells() );
00053    if( hg->parMgr()     == 0 ) hg->allocatePar( 
00054        HcalGeometry::k_NumberOfParametersPerShape*hg->numberOfShapes(),
00055       HcalGeometry::k_NumberOfParametersPerShape ) ;
00056 
00057    switch (hsub) 
00058    {
00059       case (HcalBarrel) :  fill(hsub, extTopology->firstHBRing(), extTopology->lastHBRing(), hg ); break;
00060       case (HcalEndcap) :  fill(hsub, extTopology->firstHERing(), extTopology->lastHERing(), hg ); break;
00061       case (HcalForward) : fill(hsub, extTopology->firstHFRing(), extTopology->lastHFRing(), hg ); break;
00062       case (HcalOuter) :   fill(hsub, extTopology->firstHORing(), extTopology->lastHORing(), hg ); break;
00063       default:  break;
00064    }
00065    return hg;
00066 }
00067 
00068 HcalHardcodeGeometryLoader::ReturnType
00069 HcalHardcodeGeometryLoader::load() 
00070 {
00071    ReturnType hg( new HcalGeometry( *extTopology ) ) ;
00072 
00073    if( hg->cornersMgr() == 0 ) hg->allocateCorners(  extTopology->ncells() );
00074    if( hg->parMgr()     == 0 ) hg->allocatePar( 500, 5 ) ;
00075 
00076    fill(HcalBarrel,  extTopology->firstHBRing(), extTopology->lastHBRing(), hg); 
00077    fill(HcalEndcap,  extTopology->firstHERing(), extTopology->lastHERing(), hg); 
00078    fill(HcalForward, extTopology->firstHFRing(), extTopology->lastHFRing(), hg); 
00079    fill(HcalOuter,   extTopology->firstHORing(), extTopology->lastHORing(), hg);
00080    return hg;
00081 }
00082 
00083 void 
00084 HcalHardcodeGeometryLoader::fill( HcalSubdetector subdet, 
00085                                   int             firstEtaRing, 
00086                                   int             lastEtaRing, 
00087                                   ReturnType      geom           ) 
00088 {
00089    // start by making the new HcalDetIds
00090    std::vector<HcalDetId> hcalIds;
00091    int nDepthSegments, startingDepthSegment;
00092    for(int etaRing = firstEtaRing; etaRing <= lastEtaRing; ++etaRing) {
00093       extTopology->depthBinInformation(subdet, etaRing, nDepthSegments, startingDepthSegment);
00094       unsigned int nPhiBins = extTopology->nPhiBins(etaRing);
00095       unsigned int phiInc=72/std::max(36u,nPhiBins);
00096       for(int idepth = 0; idepth < nDepthSegments; ++idepth) {
00097          int depthBin = startingDepthSegment + idepth;
00098          
00099          for(unsigned iphi = 1; iphi <= 72; iphi+=phiInc) {
00100             for(int zsign = -1; zsign <= 1; zsign += 2) {
00101                HcalDetId id( subdet, zsign * etaRing, iphi, depthBin);
00102                if (extTopology->valid(id)) hcalIds.push_back(id);
00103             }
00104          } 
00105       }
00106    }
00107 
00108    edm::LogInfo("HcalHardcodeGeometry") << "Number of HCAL DetIds made: " 
00109                                         << subdet 
00110                                         << " " << hcalIds.size();
00111    // for each new HcalDetId, make a CaloCellGeometry
00112    for(std::vector<HcalDetId>::const_iterator hcalIdItr = hcalIds.begin();
00113        hcalIdItr != hcalIds.end(); ++hcalIdItr)
00114    {
00115       makeCell(*hcalIdItr,geom) ;
00116    }
00117 }
00118      
00119 
00120 inline double theta_from_eta(double eta){return (2.0*atan(exp(-eta)));}
00121 
00122 void
00123 HcalHardcodeGeometryLoader::makeCell( const HcalDetId& detId ,
00124                                       ReturnType       geom    ) const 
00125 {
00126    // the two eta boundaries of the cell
00127    double eta1, eta2;
00128    HcalSubdetector subdet = detId.subdet();
00129    int etaRing = detId.ietaAbs();
00130    if(subdet == HcalForward) 
00131    {
00132       eta1 = theHFEtaBounds[etaRing-extTopology->firstHFRing()];
00133       eta2 = theHFEtaBounds[etaRing-extTopology->firstHFRing()+1];
00134    } 
00135    else if (etaRing==28 && detId.depth()==3) 
00136    {
00137       // double-width in eta at depth 3 in ieta=28
00138       eta1 = theHBHEEtaBounds[etaRing-1];
00139       eta2 = theHBHEEtaBounds[etaRing+1];
00140    } 
00141    else 
00142    {
00143       eta1 = theHBHEEtaBounds[etaRing-1];
00144       eta2 = theHBHEEtaBounds[etaRing];
00145    }
00146    double eta = 0.5*(eta1+eta2) * detId.zside();
00147    double deta = 0.5*(eta2-eta1);
00148    double theta = theta_from_eta(eta);
00149 
00150    int nDepthBins, startingBin;
00151    extTopology->depthBinInformation(subdet,etaRing,nDepthBins,startingBin);
00152 
00153    // in radians
00154    double dphi_nominal = 2.0*M_PI / extTopology->nPhiBins(1); // always the same
00155    double dphi_half = M_PI / extTopology->nPhiBins(etaRing); // half-width
00156   
00157    double phi_low = dphi_nominal*(detId.iphi()-1); // low-edge boundaries are constant...
00158    double phi = phi_low+dphi_half;
00159 
00160    bool isBarrel = (subdet == HcalBarrel || subdet == HcalOuter);
00161 
00162    double x,y,z,r;
00163    double thickness;
00164 
00165    if(isBarrel) 
00166    {
00167       if(subdet == HcalBarrel) 
00168       {
00169          if (detId.ietaAbs()==15) 
00170          {
00171             if (detId.depth()==1) 
00172             {
00173                r = theBarrelRadius;
00174                thickness = theHB15aThickness;
00175             } 
00176             else 
00177             {
00178                r = theBarrelRadius+theHB15aThickness;
00179                thickness = theHB15bThickness;
00180             }
00181          } 
00182          else if (detId.ietaAbs()==16) 
00183          {
00184             if (detId.depth()==1) 
00185             {
00186                r = theBarrelRadius;
00187                thickness = theHB16aThickness;
00188             } 
00189             else 
00190             {
00191                r = theBarrelRadius+theHB16aThickness;
00192                thickness = theHB16bThickness;
00193             }
00194          } 
00195          else 
00196          {
00197             r = theBarrelRadius;
00198             thickness = theHBThickness;
00199          }
00200       } 
00201       else 
00202       { // HO
00203          r = theOuterRadius;
00204          thickness = theHOThickness;
00205       } 
00206       z = r * sinh(eta); // sinh(eta) == 1/tan(theta)
00207       thickness *= cosh(eta); // cosh(eta) == 1/sin(theta)
00208    } 
00209    else 
00210    {
00211       int depth = detId.depth();
00212       if(subdet == HcalEndcap) 
00213       { // Sadly, Z must be made a function of ieta.
00214          if (nDepthBins==1) 
00215          {
00216             z = theHEZPos[depth - 1];
00217             thickness = theHEZPos[3] - z;
00218          } 
00219          else if (nDepthBins==2 && depth==2) 
00220          {
00221             z = theHEZPos[depth - 1];
00222             if (etaRing==29) // special case for tower 29
00223                thickness = theHEZPos[depth] - z;
00224             else
00225                thickness = theHEZPos[3] - z;
00226          } 
00227          else 
00228          {
00229             z = theHEZPos[depth - 1];
00230             thickness = theHEZPos[depth] - z;
00231          }
00232          thickness /= fabs(tanh(eta));
00233       } 
00234       else 
00235       {
00236          z = theHFZPos[depth - 1];
00237          thickness = theHFThickness-(theHFZPos[depth-1]-theHFZPos[0]);
00238       }
00239       z*=detId.zside(); // get the sign right.
00240       r = z * tan(theta);
00241       assert(r>0.);
00242    }
00243 
00244    x = r * cos(phi);
00245    y = r * sin(phi);
00246    GlobalPoint point(x,y,z);
00247 
00248    std::vector<CCGFloat> hp ;
00249    hp.reserve(5) ;
00250 
00251    if (subdet==HcalForward) 
00252    {
00253       hp.push_back( deta ) ;
00254       hp.push_back( dphi_half ) ;
00255       hp.push_back( thickness/2 ) ;
00256       hp.push_back( fabs( point.eta() ) ) ;
00257       hp.push_back( fabs( point.z() ) ) ;
00258    } 
00259    else 
00260    { 
00261       const double mysign ( isBarrel ? 1 : -1 ) ;
00262       hp.push_back( deta ) ;
00263       hp.push_back( dphi_half ) ;
00264       hp.push_back( mysign*thickness/2. ) ;
00265       hp.push_back( fabs( point.eta() ) ) ;
00266       hp.push_back( fabs( point.z() ) ) ;
00267    }
00268    geom->newCell( point, point, point,
00269                   CaloCellGeometry::getParmPtr( hp, 
00270                                                 geom->parMgr(), 
00271                                                 geom->parVecVec() ),
00272                   detId );
00273 }
00274 
00275