GeomBuilder/src/bSector.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author N. Amapane - INFN Torino
0005  */
0006
0007 #include "bSector.h"
0008 #include "printUniqueNames.h"
0009 #include "DataFormats/GeometrySurface/interface/Surface.h"
0010 #include "Utilities/BinningTools/interface/ClusterizingHistogram.h"
0011 #include "MagneticField/Layers/interface/MagBSector.h"
0012 #include "Utilities/General/interface/precomputed_value_sort.h"
0013
0014 #include <algorithm>
0015 #include <iostream>
0016
0017 using namespace SurfaceOrientation;
0018 using namespace std;
0019 using namespace magneticfield;
0020
0021 // Default ctor needed to have arrays.
0022 bSector::bSector() : debug(false) {}
0023
0024 // The ctor is in charge of finding rods inside the sector.
0025 bSector::bSector(handles::const_iterator begin, handles::const_iterator end, bool debugVal)
0026     : volumes(begin, end), msector(nullptr), debug(debugVal) {
0027   if (debug)
0028     cout << "   Sector at Phi  " << volumes.front()->center().phi() << " " << volumes.back()->center().phi() << endl;
0029
0030   if (volumes.size() == 1) {
0031     if (debug) {
0032       cout << "   Rod at: 0 elements: " << end - begin << " unique volumes: ";
0033       printUniqueNames(begin, end);
0034     }
0035     rods.push_back(bRod(begin, end, debug));
0036   } else {
0037     // Clusterize in phi. Use bin edge so that complete clusters can be
0038     // easily found (not trivial using bin centers!)
0039     // Unfortunately this makes the result more sensitive to the
0040     // "resolution" parameter...
0041     // To avoid +-pi boundary problems, take phi distance from min phi.
0042     // Caveat of implicit conversions of Geom::Phi!!!
0043
0044     // Sort volumes in DELTA phi - i.e. phi(j)-phi(i) > 0 if j>1.
0045     precomputed_value_sort(volumes.begin(), volumes.end(), ExtractPhiMax(), LessDPhi());
0046
0047     const float resolution(0.01);  // rad
0048     Geom::Phi<float> phi0 = volumes.front()->maxPhi();
0049     float phiMin = -resolution;  ///FIXME: (float) resolution; ??
0050
0051     // Careful -- Phi overloads arithmetic operators and will wrap around each step of a calculation,
0052     // so use casts to prevent intermediate value wrap-arounds.
0053     float phiTmp = static_cast<float>(volumes.back()->maxPhi()) - static_cast<float>(phi0) + resolution;
0054     const float phiMax = angle0to2pi::make0To2pi(phiTmp);  // Ensure 0-2pi
0055
0056     if (debug) {
0057       cout << "volumes size = " << volumes.size();
0058       cout << ", phi0 = " << phi0 << ", volumes.back()->maxPhi() = " << volumes.back()->maxPhi();
0059       cout << ", phiMin = " << phiMin << endl << "phiMax = " << phiMax;
0060       cout << ", int((phiMax - phiMin) / resolution) + 1 = " << int((phiMax - phiMin) / resolution) + 1 << endl;
0061     }
0062     ClusterizingHistogram hisPhi(int((phiMax - phiMin) / resolution) + 1, phiMin, phiMax);
0063
0064     handles::const_iterator first = volumes.begin();
0065     handles::const_iterator last = volumes.end();
0066
0067     for (handles::const_iterator i = first; i != last; ++i) {
0068       hisPhi.fill((*i)->maxPhi() - phi0);
0069     }
0070     vector<float> phiClust = hisPhi.clusterize(resolution);
0071
0072     if (debug)
0073       cout << "     Found " << phiClust.size() << " clusters in Phi, "
0074            << " rods: " << endl;
0075
0076     handles::const_iterator rodStart = first;
0077     handles::const_iterator separ = first;
0078
0079     float DZ = (*max_element(first, last, LessZ()))->maxZ() - (*min_element(first, last, LessZ()))->minZ();
0080
0081     float DZ1 = 0.;
0082     for (unsigned int i = 0; i < phiClust.size(); ++i) {
0083       float phiSepar;
0084       if (i < phiClust.size() - 1) {
0085         phiSepar = (phiClust[i] + phiClust[i + 1]) / 2.f;
0086       } else {
0087         phiSepar = phiMax;
0088       }
0089       if (debug)
0090         cout << "       cluster " << i << " phisepar " << phiSepar << endl;
0091       while (separ < last && (*separ)->maxPhi() - phi0 < phiSepar) {
0092         DZ1 += ((*separ)->maxZ() - (*separ)->minZ());
0093         if (debug)
0094           cout << "         " << (*separ)->name << " " << (*separ)->maxPhi() - phi0 << " " << (*separ)->maxZ() << " "
0095                << (*separ)->minZ() << " " << DZ1 << endl;
0096         ++separ;
0097       }
0098
0099       // FIXME: print warning for small discrepancies. Tolerance (below)
0100       // had to be increased since discrepancies sum to up to ~ 2 mm.
0101       if (fabs(DZ - DZ1) > 0.001 && fabs(DZ - DZ1) < 0.5) {
0102         if (debug)
0103           cout << "*** WARNING: Z lenght mismatch by " << DZ - DZ1 << " " << DZ << " " << DZ1 << endl;
0104       }
0105       if (fabs(DZ - DZ1) > 0.25) {  // FIXME hardcoded tolerance
0106         if (debug)
0107           cout << "       Incomplete, use also next cluster: " << DZ << " " << DZ1 << " " << DZ - DZ1 << endl;
0108         DZ1 = 0.;
0109         continue;
0110       } else if (DZ1 > DZ + 0.05) {  // Wrong: went past max lenght // FIXME hardcoded tolerance
0111         cout << " *** ERROR: bSector finding messed up." << endl;
0112         printUniqueNames(rodStart, separ);
0113         DZ1 = 0.;
0114       } else {
0115         if (debug) {
0116           cout << "       Rod at: " << phiClust[i] << " elements: " << separ - rodStart << " unique volumes: ";
0117           printUniqueNames(rodStart, separ);
0118         }
0119
0120         rods.push_back(bRod(rodStart, separ, debug));
0121         rodStart = separ;
0122         DZ1 = 0.;
0123       }
0124     }
0125
0126     if (rods.empty())
0127       cout << " *** ERROR: bSector has no rods " << DZ << " " << DZ1 << endl;
0128     if (debug)
0129       cout << "-----------------------" << endl;
0130   }
0131 }
0132
0133 MagBSector* bSector::buildMagBSector() const {
0134   if (msector == nullptr) {
0135     vector<MagBRod*> mRods;
0136     for (vector<bRod>::const_iterator rod = rods.begin(); rod != rods.end(); ++rod) {
0137       mRods.push_back((*rod).buildMagBRod());
0138     }
0139     msector = new MagBSector(mRods, volumes.front()->minPhi());  //FIXME
0140     // Never deleted. When is it safe to delete it?
0141   }
0142   return msector;
0143 }