CSCGeometryBuilder/src/CSCGeometryBuilder.cc

0001 #include "CSCGeometryBuilder.h"
0002
0003 #include <Geometry/CSCGeometry/interface/CSCGeometry.h>
0004
0005 #include <DataFormats/DetId/interface/DetId.h>
0006 #include <DataFormats/MuonDetId/interface/CSCDetId.h>
0007 #include <Geometry/CSCGeometry/interface/CSCWireGroupPackage.h>
0008
0009 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0010
0011 #include <vector>
0012
0013 CSCGeometryBuilder::CSCGeometryBuilder() : myName("CSCGeometryBuilder") {}
0014
0015 CSCGeometryBuilder::~CSCGeometryBuilder() {}
0016
0017 void CSCGeometryBuilder::build(CSCGeometry& theGeometry,
0018                                const RecoIdealGeometry& rig,
0019                                const CSCRecoDigiParameters& cscpars) {
0020   std::vector<float> fpar;
0021   std::vector<float> gtran;
0022   std::vector<float> grmat;
0023   std::vector<float> fupar;
0024   std::vector<double>::const_iterator it, endIt;
0025   const std::vector<DetId>& detids(rig.detIds());
0026
0027   for (size_t idt = 0; idt < detids.size(); ++idt) {
0028     CSCDetId detid = CSCDetId(detids[idt]);
0029     int jstation = detid.station();
0030     int jring = detid.ring();
0031
0032     endIt = rig.shapeEnd(idt);
0033     fpar.clear();
0034     for (it = rig.shapeStart(idt); it != endIt; ++it) {
0035       fpar.emplace_back((float)(*it));
0036     }
0037
0038     gtran.clear();
0039     endIt = rig.tranEnd(idt);
0040     for (it = rig.tranStart(idt); it != endIt; ++it) {
0041       gtran.emplace_back((float)(*it));
0042     }
0043     grmat.clear();
0044     endIt = rig.rotEnd(idt);
0045     for (it = rig.rotStart(idt); it != endIt; ++it) {
0046       grmat.emplace_back((float)(*it));
0047     }
0048
0049     // get the chamber type from existing info
0050     int chamberType = CSCChamberSpecs::whatChamberType(jstation, jring);
0051     size_t cs = 0;
0052     assert(!cscpars.pChamberType.empty());
0053     while (cs < cscpars.pChamberType.size() && chamberType != cscpars.pChamberType[cs]) {
0054       ++cs;
0055     }
0056     assert(cs != cscpars.pChamberType.size());
0057
0058     // check the existence of the specs for this type WHY? Remove it...
0059     size_t fu, numfuPars;
0060     CSCWireGroupPackage wg;
0061     fu = cscpars.pUserParOffset[cs];
0062     numfuPars = fu + 1 + size_t(cscpars.pfupars[fu]);
0063
0064     // upars from db are now uparvals + wg info so we need to unwrap only part here first...
0065     LogTrace(myName) << myName << ": I think I have " << cscpars.pUserParSize[cs] << " values in pfupars (uparvals)."
0066                      << std::endl;
0067     LogTrace(myName) << myName << ": For fupar I will start at " << cscpars.pUserParOffset[cs] + 1
0068                      << " in pfupars and go to " << numfuPars << "." << std::endl;
0069     for (++fu; fu < numfuPars; ++fu) {
0070       LogTrace(myName) << myName << ": pfupars[" << fu << "]=" << cscpars.pfupars[fu] << std::endl;
0071       fupar.emplace_back(cscpars.pfupars[fu]);
0072     }
0073     // now, we need to start from "here" at fu to go on and build wg...
0074     wg.wireSpacing = cscpars.pfupars[fu++];
0075     wg.alignmentPinToFirstWire = cscpars.pfupars[fu++];
0076     wg.numberOfGroups = int(cscpars.pfupars[fu++]);
0077     wg.narrowWidthOfWirePlane = cscpars.pfupars[fu++];
0078     wg.wideWidthOfWirePlane = cscpars.pfupars[fu++];
0079     wg.lengthOfWirePlane = cscpars.pfupars[fu++];
0080     size_t numgrp = static_cast<size_t>(cscpars.pfupars[fu]);
0081     size_t maxFu = fu + 1 + numgrp;
0082     fu++;
0083     for (; fu < maxFu; ++fu) {
0084       wg.wiresInEachGroup.emplace_back(int(cscpars.pfupars[fu]));
0085     }
0086     maxFu = fu + numgrp;
0087     for (; fu < maxFu; ++fu) {
0088       wg.consecutiveGroups.emplace_back(int(cscpars.pfupars[fu]));
0089     }
0090
0091     if (wg.numberOfGroups != 0) {
0092       LogTrace(myName) << myName << ": TotNumWireGroups     = " << wg.numberOfGroups;
0093       LogTrace(myName) << myName << ": WireSpacing          = " << wg.wireSpacing;
0094       LogTrace(myName) << myName << ": AlignmentPinToFirstWire = " << wg.alignmentPinToFirstWire;
0095       LogTrace(myName) << myName << ": Narrow width of wire plane = " << wg.narrowWidthOfWirePlane;
0096       LogTrace(myName) << myName << ": Wide width of wire plane = " << wg.wideWidthOfWirePlane;
0097       LogTrace(myName) << myName << ": Length in y of wire plane = " << wg.lengthOfWirePlane;
0098       LogTrace(myName) << myName << ": wg.consecutiveGroups.size() = " << wg.consecutiveGroups.size();
0099       LogTrace(myName) << myName << ": wg.wiresInEachGroup.size() = " << wg.wiresInEachGroup.size();
0100       LogTrace(myName) << myName << ": \tNumGroups\tWiresInGroup";
0101       for (size_t i = 0; i < wg.consecutiveGroups.size(); i++) {
0102         LogTrace(myName) << myName << " \t" << wg.consecutiveGroups[i] << "\t\t" << wg.wiresInEachGroup[i];
0103       }
0104     } else {
0105       LogTrace(myName) << myName << ": DDD is MISSING SpecPars for wire groups";
0106     }
0107     LogTrace(myName) << myName << ": end of wire group info. ";
0108
0109     // Are we going to apply centre-to-intersection offsets, even if values exist in the specs file?
0110     if (!theGeometry.centreTIOffsets())
0111       fupar[30] = 0.;  // reset to zero if flagged 'off'
0112
0113     buildChamber(theGeometry, detid, fpar, fupar, gtran, grmat, wg);  //, cscpars.pWGPs[cs] );
0114     fupar.clear();
0115   }
0116 }
0117
0118 void CSCGeometryBuilder::buildChamber(CSCGeometry& theGeometry,         // the geometry container
0119                                       CSCDetId chamberId,               // the DetId for this chamber
0120                                       const std::vector<float>& fpar,   // volume parameters hB, hT. hD, hH
0121                                       const std::vector<float>& fupar,  // user parameters
0122                                       const std::vector<float>& gtran,  // translation vector
0123                                       const std::vector<float>& grmat,  // rotation matrix
0124                                       const CSCWireGroupPackage& wg     // wire group info
0125 ) {
0126   LogTrace(myName) << myName << ": entering buildChamber";
0127
0128   int jend = chamberId.endcap();
0129   int jstat = chamberId.station();
0130   int jring = chamberId.ring();
0131   int jch = chamberId.chamber();
0132   int jlay = chamberId.layer();
0133
0134   if (jlay != 0)
0135     edm::LogWarning(myName) << "Error! CSCGeometryBuilderFromDDD was fed layer id = " << jlay << "\n";
0136
0137   const size_t kNpar = 4;
0138   if (fpar.size() != kNpar)
0139     edm::LogError(myName) << "Error, expected npar=" << kNpar << ", found npar=" << fpar.size() << std::endl;
0140
0141   LogTrace(myName) << myName << ":  E" << jend << " S" << jstat << " R" << jring << " C" << jch << " L" << jlay;
0142   LogTrace(myName) << myName << ": npar=" << fpar.size() << " hB=" << fpar[0] << " hT=" << fpar[1] << " hD=" << fpar[2]
0143                    << " hH=" << fpar[3];
0144   LogTrace(myName) << myName << ": gtran[0,1,2]=" << gtran[0] << " " << gtran[1] << " " << gtran[2];
0145   LogTrace(myName) << myName << ": grmat[0-8]=" << grmat[0] << " " << grmat[1] << " " << grmat[2] << " " << grmat[3]
0146                    << " " << grmat[4] << " " << grmat[5] << " " << grmat[6] << " " << grmat[7] << " " << grmat[8];
0147   LogTrace(myName) << myName << ": nupar=" << fupar.size() << " upar[0]=" << fupar[0] << " upar[" << fupar.size() - 1
0148                    << "]=" << fupar[fupar.size() - 1];
0149
0150   const CSCChamber* chamber = theGeometry.chamber(chamberId);
0151   if (chamber) {
0152   } else {  // this chamber not yet built/stored
0153
0154     LogTrace(myName) << myName << ": CSCChamberSpecs::build requested for ME" << jstat << jring;
0155     int chamberType = CSCChamberSpecs::whatChamberType(jstat, jring);
0156     const CSCChamberSpecs* aSpecs = theGeometry.findSpecs(chamberType);
0157     if (!fupar.empty() && aSpecs == nullptr) {
0158       // make new one:
0159       aSpecs = theGeometry.buildSpecs(chamberType, fpar, fupar, wg);
0160     } else if (fupar.empty() && aSpecs == nullptr) {
0161       edm::LogError(myName)
0162           << "SHOULD BE THROW? Error, wg and/or fupar size are 0 BUT this Chamber Spec has not been built!";
0163     }
0164
0165     // Build a Transformation out of GEANT gtran and grmat...
0166     // These are used to transform a point in the local reference frame
0167     // of a subdetector to the global frame of CMS by
0168     //         (grmat)^(-1)*local + (gtran)
0169     // The corresponding transformation from global to local is
0170     //         (grmat)*(global - gtran)
0171
0172     Surface::RotationType aRot(
0173         grmat[0], grmat[1], grmat[2], grmat[3], grmat[4], grmat[5], grmat[6], grmat[7], grmat[8]);
0174
0175     // This rotation from GEANT considers the detector face as the x-z plane.
0176     // We want this to be the local x-y plane.
0177     // Furthermore, the -z_global endcap has LH local coordinates, since it is built
0178     // in GEANT as a *reflection* of the +z_global endcap.
0179     // So we need to rotate, and in -z flip local x.
0180
0181     // aRot.rotateAxes will transform aRot in place so that it becomes
0182     // applicable to the new local coordinates: detector face in x-y plane
0183     // looking out along z, in either endcap.
0184
0185     // The interface for rotateAxes specifies 'new' X,Y,Z but the
0186     // implementation deals with them as the 'old'.
0187
0188     Basic3DVector<float> oldX(1., 0., 0.);
0189     Basic3DVector<float> oldY(0., 0., -1.);
0190     Basic3DVector<float> oldZ(0., 1., 0.);
0191
0192     if (gtran[2] < 0.)
0193       oldX *= -1;
0194
0195     aRot.rotateAxes(oldX, oldY, oldZ);
0196
0197     // Need to know z of layers w.r.t to z of centre of chamber.
0198
0199     float frameThickness = fupar[31] / 10.;   // mm -> cm
0200     float gapThickness = fupar[32] / 10.;     // mm -> cm
0201     float panelThickness = fupar[33] / 10.;   // mm -> cm
0202     float zAverageAGVtoAF = fupar[34] / 10.;  // mm -> cm
0203
0204     float layerThickness = gapThickness;                    // consider the layer to be the gas gap
0205     float layerSeparation = gapThickness + panelThickness;  // centre-to-centre of neighbouring layers
0206
0207     float chamberThickness = 7. * panelThickness + 6. * gapThickness + 2. * frameThickness;  // chamber frame thickness
0208     float hChamberThickness = chamberThickness / 2.;  // @@ should match value returned from DDD directly
0209
0210     // distAverageAGVtoAF is offset between centre of chamber (AF) and (L1+L6)/2 (average AGVs)
0211     // where AF = AluminumFrame and AGV=ActiveGasVolume (volume names in DDD).
0212     // It is signed based on global z values: zc - (zl1+zl6)/2
0213
0214     // Local z values w.r.t. AF...
0215     //     z of wires in layer 1 = z_w1 = +/- zAverageAGVtoAF + 2.5*layerSeparation; // layer 1 is at most +ve local z
0216     // The sign in '+/-' depends on relative directions of local and global z.
0217     // It is '-' if they are the same direction, and '+' if opposite directions.
0218     //     z of wires in layer N   = z_wN = z_w1 - (N-1)*layerSeparation;
0219     //     z of strips in layer N  = z_sN = z_wN + gapThickness/2.; @@ BEWARE: need to check if it should be '-gapThickness/2' !
0220
0221     // Set dimensions of trapezoidal chamber volume
0222     // N.B. apothem is 4th in fpar but 3rd in ctor
0223
0224     // hChamberThickness and fpar[2] should be the same - but using the above value at least shows
0225     // how chamber structure works
0226
0227     TrapezoidalPlaneBounds* bounds = new TrapezoidalPlaneBounds(fpar[0], fpar[1], fpar[3], hChamberThickness);
0228
0229     // Centre of chamber in z is specified in DDD
0230     Surface::PositionType aVec(gtran[0], gtran[1], gtran[2]);
0231
0232     Plane::PlanePointer plane = Plane::build(aVec, aRot, bounds);
0233
0234     CSCChamber* chamber = new CSCChamber(plane, chamberId, aSpecs);
0235     theGeometry.addChamber(chamber);
0236
0237     LogTrace(myName) << myName << ": Create chamber E" << jend << " S" << jstat << " R" << jring << " C" << jch
0238                      << " z=" << gtran[2] << " t/2=" << fpar[2] << " (DDD) or " << hChamberThickness
0239                      << " (specs) adr=" << chamber;
0240
0241     // Create the component layers of this chamber
0242     // We're taking the z as the z of the wire plane within the layer (middle of gas gap)
0243
0244     // Specify global z of layer by offsetting from centre of chamber: since layer 1
0245     // is nearest to IP in stations 1/2 but layer 6 is nearest in stations 3/4,
0246     // we need to adjust sign of offset appropriately...
0247     int localZwrtGlobalZ = +1;
0248     if ((jend == 1 && jstat < 3) || (jend == 2 && jstat > 2))
0249       localZwrtGlobalZ = -1;
0250     int globalZ = +1;
0251     if (jend == 2)
0252       globalZ = -1;
0253
0254     LogTrace(myName) << myName << ": layerSeparation=" << layerSeparation << ", zAF-zAverageAGV=" << zAverageAGVtoAF
0255                      << ", localZwrtGlobalZ=" << localZwrtGlobalZ << ", gtran[2]=" << gtran[2];
0256
0257     for (short j = 1; j <= 6; ++j) {
0258       CSCDetId layerId = CSCDetId(jend, jstat, jring, jch, j);
0259
0260       // extra-careful check that we haven't already built this layer
0261       const CSCLayer* cLayer = dynamic_cast<const CSCLayer*>(theGeometry.idToDet(layerId));
0262
0263       if (cLayer == nullptr) {
0264         // build the layer - need the chamber's specs and an appropriate layer-geometry
0265         const CSCChamberSpecs* aSpecs = chamber->specs();
0266         const CSCLayerGeometry* geom = (j % 2 != 0) ? aSpecs->oddLayerGeometry(jend) : aSpecs->evenLayerGeometry(jend);
0267
0268         // Build appropriate BoundPlane, based on parent chamber, with gas gap as thickness
0269
0270         // centre of chamber is at global z = gtran[2]
0271         float zlayer = gtran[2] - globalZ * zAverageAGVtoAF + localZwrtGlobalZ * (3.5 - j) * layerSeparation;
0272
0273         Surface::RotationType chamberRotation = chamber->surface().rotation();
0274         Surface::PositionType layerPosition(gtran[0], gtran[1], zlayer);
0275         std::array<const float, 4> const& dims = geom->parameters();  // returns hb, ht, d, a
0276         // dims[2] = layerThickness/2.; // half-thickness required and note it is 3rd value in vector
0277         TrapezoidalPlaneBounds* bounds = new TrapezoidalPlaneBounds(dims[0], dims[1], dims[3], layerThickness / 2.);
0278         Plane::PlanePointer plane = Plane::build(layerPosition, chamberRotation, bounds);
0279
0280         CSCLayer* layer = new CSCLayer(plane, layerId, chamber, geom);
0281
0282         LogTrace(myName) << myName << ": Create layer E" << jend << " S" << jstat << " R" << jring << " C" << jch
0283                          << " L" << j << " z=" << zlayer << " t=" << layerThickness << " or "
0284                          << layer->surface().bounds().thickness() << " adr=" << layer << " layerGeom adr=" << geom;
0285
0286         chamber->addComponent(j, layer);
0287         theGeometry.addLayer(layer);
0288       } else {
0289         edm::LogError(myName) << ": ERROR, layer " << j << " for chamber = " << (chamber->id())
0290                               << " already exists: layer address=" << cLayer << " chamber address=" << chamber << "\n";
0291       }
0292
0293     }  // layer construction within chamber
0294   }  // chamber construction
0295 }