File indexing completed on 2024-04-06 12:22:29
0001
0002
0003
0004
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
0022 bSector::bSector() : debug(false) {}
0023
0024
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
0038
0039
0040
0041
0042
0043
0044
0045 precomputed_value_sort(volumes.begin(), volumes.end(), ExtractPhiMax(), LessDPhi());
0046
0047 const float resolution(0.01);
0048 Geom::Phi<float> phi0 = volumes.front()->maxPhi();
0049 float phiMin = -resolution;
0050
0051
0052
0053 float phiTmp = static_cast<float>(volumes.back()->maxPhi()) - static_cast<float>(phi0) + resolution;
0054 const float phiMax = angle0to2pi::make0To2pi(phiTmp);
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
0100
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) {
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) {
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());
0140
0141 }
0142 return msector;
0143 }