GeomBuilder/src/bLayer.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author N. Amapane - INFN Torino
0005  */
0006
0007 #include "bLayer.h"
0008 #include "printUniqueNames.h"
0009 #include "MagneticField/VolumeGeometry/interface/MagVolume6Faces.h"
0010 #include "MagneticField/Layers/interface/MagBLayer.h"
0011
0012 #include "Utilities/General/interface/precomputed_value_sort.h"
0013
0014 #include <iostream>
0015
0016 using namespace SurfaceOrientation;
0017 using namespace magneticfield;
0018
0019 //The ctor is in charge of finding sectors inside the layer.
0020 bLayer::bLayer(handles::const_iterator begin, handles::const_iterator end, bool debugFlag)
0021     : size(end - begin), theVolumes(begin, end), mlayer(nullptr), debug(debugFlag) {
0022   // Sort in phi
0023   precomputed_value_sort(theVolumes.begin(), theVolumes.end(), ExtractPhi());
0024
0025   if (debug) {
0026     std::cout << " elements: " << theVolumes.size() << " unique volumes: ";
0027     printUniqueNames(theVolumes.begin(), theVolumes.end());
0028   }
0029
0030   // Find sectors in phi
0031   handles::iterator secBegin = theVolumes.begin();
0032   handles::iterator secEnd = secBegin;
0033   int binOffset = 0;
0034
0035   const Surface& refSurf = (*secBegin)->surface(outer);
0036
0037   int newbegin = 0;
0038   int newend = 0;
0039
0040   // A sector is made of several volumes in R, and, for planar layers
0041   // (box and traps) also in phi, so it might cross the -phi boundary.
0042   // For those, we have to look for the end of first sector and rotate the
0043   // vector of volumes.
0044   // ASSUMPTION: all volumes in a layer must be of compatible type.
0045
0046   if (size == 1) {  // Only one volume; this is the case for barrel
0047     // cylinders.
0048     // FIXME sectors.push_back(bSector(theVolumes.begin(),theVolumes.end());
0049     if (debug)
0050       std::cout << "      Sector is just one volume." << std::endl;
0051
0052   } else if (size == 12 ||  // In this case, each volume is a sector.
0053              (((*secBegin)->shape() != DDSolidShape::ddtrap) && (*secBegin)->shape() != DDSolidShape::ddbox)) {
0054     secEnd = secBegin + size / 12;
0055
0056   } else {                    // there are more than one volume per sector.
0057     float tolerance = 0.025;  // 250 micron
0058     do {
0059       if (debug)
0060         std::cout << (*secBegin)->name << " " << (*secBegin)->copyno << std::endl;
0061       ++secBegin;
0062     } while (
0063         (secBegin != theVolumes.end()) &&
0064         (*secBegin)->sameSurface(
0065             refSurf,
0066             outer,
0067             tolerance));  // This works only if outer surface is a plane, otherwise sameSurface returns always true!
0068
0069     secEnd = secBegin;
0070     secBegin = theVolumes.begin() + bin((secEnd - theVolumes.begin()) - size / 12);
0071     ;
0072     newend = secEnd - theVolumes.begin();
0073     newbegin = secBegin - theVolumes.begin();
0074
0075     // Rotate the begin of the first sector to the vector beginning.
0076     rotate(theVolumes.begin(), secBegin, theVolumes.end());
0077     secBegin = theVolumes.begin();
0078     secEnd = secBegin + size / 12;
0079
0080     // Test it is correct...
0081     if (!((*secBegin)->sameSurface((*(secEnd - 1))->surface(outer), outer, tolerance))) {
0082       std::cout << "*** ERROR: Big mess while looking for sectors " << (*secBegin)->name << " " << (*secBegin)->copyno
0083                 << " " << (*(secEnd - 1))->name << " " << (*(secEnd - 1))->copyno << std::endl;
0084     }
0085   }
0086
0087   if (debug) {
0088     std::cout << "      First sector: volumes " << secEnd - theVolumes.begin() << " from " << newbegin
0089               << " (phi = " << (*secBegin)->center().phi() << ") "
0090               << " to " << newend << " (phi = " << (*secEnd)->center().phi() << ") "
0091               << " # " << (*secBegin)->copyno << " ";
0092     std::cout << GlobalVector(refSurf.rotation().zx(), refSurf.rotation().zy(), refSurf.rotation().zz()) << std::endl;
0093   }
0094
0095   if (size != 1) {  // Build the 12 sectors
0096     int offset = size / 12;
0097     sectors.resize(12);
0098     for (int i = 0; i < 12; ++i) {
0099       int isec = (i + binOffset) % 12;
0100       sectors[isec >= 0 ? isec : isec + 12] =
0101           bSector(theVolumes.begin() + ((i)*offset), theVolumes.begin() + ((i + 1) * offset), debug);
0102     }
0103   }
0104
0105   if (debug)
0106     std::cout << "-----------------------" << std::endl;
0107 }
0108
0109 int bLayer::bin(int i) const {
0110   i = i % size;
0111   return (i >= 0 ? i : i + size);
0112 }
0113
0114 // const bSector &
0115 // bLayer::sector(int i) const {
0116 //   i = i%12;
0117 //   return sectors[i>=0?i:i+12];
0118 // }
0119
0120 double bLayer::minR() const {
0121   // ASSUMPTION: a layer is only 1 volume thick (by construction).
0122   return theVolumes.front()->minR();
0123 }
0124
0125 // double bLayer::maxR() const {
0126 //   // ASSUMPTION: a layer is only 1 volume thick (by construction).
0127 //   return theVolumes.front()->maxR();
0128 // }
0129
0130 MagBLayer* bLayer::buildMagBLayer() const {
0131   if (mlayer == nullptr) {
0132     // If we have only one volume, do not build any MagBSector.
0133     if (sectors.empty()) {
0134       if (debug && size != 0) {
0135         std::cout << "ERROR: bLayer::buildMagBLayer, 0 sectors but " << size << " volumes" << std::endl;
0136       }
0137       // Technically we might have only one bSector built and we would
0138       // not need a separate MagBLayer constructor...
0139       mlayer = new MagBLayer(theVolumes.front()->magVolume, minR());
0140     }
0141
0142     // If we have several sectors, create the MagBSector
0143     std::vector<MagBSector*> mSectors;
0144     for (unsigned int i = 0; i < sectors.size(); ++i) {
0145       mSectors.push_back(sectors[i].buildMagBSector());
0146     }
0147     mlayer = new MagBLayer(mSectors, minR());
0148   }
0149   return mlayer;
0150 }