TkMSParametrization/src/MSLayersKeeperX0Averaged.cc

0001 #include "MSLayersKeeperX0Averaged.h"
0002 #include "MSLayersKeeperX0AtEta.h"
0003 #include "MultipleScatteringGeometry.h"
0004
0005 using namespace std;
0006
0007 MSLayersKeeperX0Averaged::MSLayersKeeperX0Averaged(const GeometricSearchTracker &tracker, const MagneticField &bfield) {
0008   //  cout << "HERE INITIALISATION! MSLayersKeeperX0Averaged"<<endl;
0009   MSLayersKeeperX0AtEta layersX0Eta(tracker, bfield);
0010   MultipleScatteringGeometry geom(tracker);
0011   vector<MSLayer> allLayers = geom.detLayers();
0012   vector<MSLayer>::iterator it;
0013   for (int i = -1; i <= 1; i++) {
0014     float eta = i * (-1.8);
0015     vector<MSLayer> tmpLayers = geom.otherLayers(eta);
0016     vector<MSLayer>::const_iterator ic;
0017     for (ic = tmpLayers.begin(); ic != tmpLayers.end(); ic++) {
0018       it = find(allLayers.begin(), allLayers.end(), *ic);
0019       if (it == allLayers.end())
0020         allLayers.push_back(*ic);
0021     }
0022   }
0023
0024   for (it = allLayers.begin(); it != allLayers.end(); it++) {
0025     float cotTheta = (it->face() == GeomDetEnumerators::barrel) ? it->range().mean() / it->position()
0026                                                                 : it->position() / it->range().mean();
0027
0028     int nbins = 0;
0029     float sumX0 = 0.;
0030     for (int ibin = 0; ibin < 2 * layersX0Eta.theHalfNBins; ibin++) {
0031       const MSLayersAtAngle &layers = layersX0Eta.theLayersData[ibin];
0032       const MSLayer *aLayer = layers.findLayer(*it);
0033       if (aLayer) {
0034         nbins++;
0035         sumX0 += getDataX0(*aLayer).x0;
0036       }
0037     }
0038     if (nbins == 0)
0039       nbins = 1;
0040
0041     float hrange = (it->range().max() - it->range().min()) / 2.;
0042     DataX0 dataX0;
0043     if (it->face() == GeomDetEnumerators::endcap) {
0044       float cot1 = it->position() / (it->range().mean() - hrange / 2);
0045       float cot2 = it->position() / (it->range().mean() + hrange / 2);
0046       const MSLayer *aLayer1 = layersX0Eta.layers(cot1).findLayer(*it);
0047       const MSLayer *aLayer2 = layersX0Eta.layers(cot1).findLayer(*it);
0048       float sum1 = aLayer1 ? aLayer1->sumX0D(cot1) : 0.;
0049       float sum2 = aLayer2 ? aLayer2->sumX0D(cot2) : 0.;
0050       float slope = (sum2 - sum1) / (1 / cot2 - 1 / cot1);
0051       float sumX0D = sum1 + slope * (1 / cotTheta - 1 / cot1);
0052       dataX0 = DataX0(sumX0 / nbins, sumX0D, cotTheta);
0053       dataX0.setForwardSumX0DSlope(slope);
0054     } else {
0055       float sumX0D = 0;
0056       int nb = 10;
0057       for (int i = 0; i < nb; i++) {
0058         float cot = (it->range().mean() + (2 * i + 1 - nb) * hrange / nb) / it->position();
0059         float sin = 1 / sqrt(1 + cot * cot);
0060         const MSLayer *aLayer = layersX0Eta.layers(cot).findLayer(*it);
0061         if (aLayer)
0062           sumX0D += aLayer->sumX0D(cot) * sqrt(sin);
0063       }
0064       dataX0 = DataX0(sumX0 / nbins, sumX0D / nb, 0);
0065     }
0066     setDataX0(*it, dataX0);
0067   }
0068   theLayersData = MSLayersAtAngle(allLayers);
0069   //  cout << "MSLayersKeeperX0Averaged - LAYERS:"<<endl;
0070   //  theLayersData.print();
0071   //  cout << "END OF LAYERS"<<endl;
0072 }