TkMSParametrization/src/MSLayersKeeperX0DetLayer.cc

0001 #include "MSLayersKeeperX0DetLayer.h"
0002 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0003 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0004 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0005 #include "DataFormats/GeometrySurface/interface/BoundSurface.h"
0006 #include "DataFormats/GeometrySurface/interface/MediumProperties.h"
0007 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0008
0009 #include <vector>
0010
0011 using namespace std;
0012
0013 MSLayersKeeperX0DetLayer::MSLayersKeeperX0DetLayer() {
0014   //  vector<MSLayer> allLayers = MSLayersKeeperX0DetLayerGeom().detLayers();
0015   //MP
0016   vector<MSLayer> allLayers;
0017   theLayersData = MSLayersAtAngle(allLayers);
0018
0019   vector<MSLayer>::iterator it;
0020   PixelRecoPointRZ zero(0., 0.);
0021   for (it = allLayers.begin(); it != allLayers.end(); it++) {
0022     PixelRecoPointRZ middle = it->face() == GeomDetEnumerators::barrel
0023                                   ? PixelRecoPointRZ(it->position(), it->range().mean())
0024                                   : PixelRecoPointRZ(it->range().mean(), it->position());
0025
0026     float cotTheta = PixelRecoLineRZ(zero, middle).cotLine();
0027     float x0 = getDataX0(*it).x0;
0028
0029     DataX0 dataX0;
0030     if (it->face() == GeomDetEnumerators::barrel) {
0031       float sumX0D = theLayersData.sumX0D(zero, middle);
0032       dataX0 = DataX0(x0, sumX0D, cotTheta);
0033     } else {
0034       float hrange = (it->range().max() - it->range().min()) / 2.;
0035       float cot1 = it->position() / (it->range().mean() - hrange / 2);
0036       float cot2 = it->position() / (it->range().mean() + hrange / 2);
0037       PixelRecoLineRZ line1(zero, cot1);
0038       PixelRecoLineRZ line2(zero, cot2);
0039       float sum1 = theLayersData.sumX0D(zero, it->crossing(line1).first);
0040       float sum2 = theLayersData.sumX0D(zero, it->crossing(line2).first);
0041       float slope = (sum2 - sum1) / (1 / cot2 - 1 / cot1);
0042       float sumX0D = sum1 + slope * (1 / cotTheta - 1 / cot1);
0043       dataX0 = DataX0(x0, sumX0D, cotTheta);
0044       dataX0.setForwardSumX0DSlope(slope);
0045     }
0046     setDataX0(*it, dataX0);
0047     theLayersData.update(*it);
0048   }
0049   //cout << "MSLayersKeeperX0DetLayer LAYERS: "<<endl;
0050   //theLayersData.print();
0051 }
0052
0053 // vector<MSLayer>
0054 // MSLayersKeeperX0DetLayer::MSLayersKeeperX0DetLayerGeom::detLayers() const
0055 // {
0056
0057 //   vector<MSLayer> result;
0058
0059 //   vector<const DetLayer*>::const_iterator it;
0060 //   for (it = theLayers.begin(); it != theLayers.end(); it++) {
0061
0062 //     //    const DetUnit * du = (*it)->detUnits().front();
0063 //     const GeomDetUnit * du;
0064 //     //MP how access the geomdetunit??
0065 //    const MediumProperties * mp = du->surface().mediumProperties();
0066 //     float x0 = (mp) ? mp->radLen() : 0.03;
0067 //     cout << "MediumProperties: "<<mp<<" x0="<<x0<<endl;
0068 //     MSLayer layer(*it, DataX0(x0,0,0));
0069 //     result.push_back( layer);
0070 //   }
0071 //   return result;
0072 // }