Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:14:45

0001 #include "ShallowGainCalibration.h"
0002 
0003 using namespace edm;
0004 using namespace reco;
0005 using namespace std;
0006 
0007 ShallowGainCalibration::ShallowGainCalibration(const edm::ParameterSet& iConfig)
0008     : tracks_token_(consumes<edm::View<reco::Track>>(iConfig.getParameter<edm::InputTag>("Tracks"))),
0009       association_token_(consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("Tracks"))),
0010       trackerGeometry_token_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0011       gain_token_(esConsumes<SiStripGain, SiStripGainRcd>()),
0012       tkGeom_token_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0013       Suffix(iConfig.getParameter<std::string>("Suffix")),
0014       Prefix(iConfig.getParameter<std::string>("Prefix")) {
0015   produces<std::vector<int>>(Prefix + "trackindex" + Suffix);
0016   produces<std::vector<unsigned int>>(Prefix + "rawid" + Suffix);
0017   produces<std::vector<double>>(Prefix + "localdirx" + Suffix);
0018   produces<std::vector<double>>(Prefix + "localdiry" + Suffix);
0019   produces<std::vector<double>>(Prefix + "localdirz" + Suffix);
0020   produces<std::vector<unsigned short>>(Prefix + "firststrip" + Suffix);
0021   produces<std::vector<unsigned short>>(Prefix + "nstrips" + Suffix);
0022   produces<std::vector<bool>>(Prefix + "saturation" + Suffix);
0023   produces<std::vector<bool>>(Prefix + "overlapping" + Suffix);
0024   produces<std::vector<bool>>(Prefix + "farfromedge" + Suffix);
0025   produces<std::vector<unsigned int>>(Prefix + "charge" + Suffix);
0026   produces<std::vector<double>>(Prefix + "path" + Suffix);
0027 #ifdef ExtendedCALIBTree
0028   produces<std::vector<double>>(Prefix + "chargeoverpath" + Suffix);
0029 #endif
0030   produces<std::vector<unsigned char>>(Prefix + "amplitude" + Suffix);
0031   produces<std::vector<double>>(Prefix + "gainused" + Suffix);
0032   produces<std::vector<double>>(Prefix + "gainusedTick" + Suffix);
0033 }
0034 
0035 void ShallowGainCalibration::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0036   auto trackindex = std::make_unique<std::vector<int>>();
0037   auto rawid = std::make_unique<std::vector<unsigned int>>();
0038   auto localdirx = std::make_unique<std::vector<double>>();
0039   auto localdiry = std::make_unique<std::vector<double>>();
0040   auto localdirz = std::make_unique<std::vector<double>>();
0041   auto firststrip = std::make_unique<std::vector<unsigned short>>();
0042   auto nstrips = std::make_unique<std::vector<unsigned short>>();
0043   auto saturation = std::make_unique<std::vector<bool>>();
0044   auto overlapping = std::make_unique<std::vector<bool>>();
0045   auto farfromedge = std::make_unique<std::vector<bool>>();
0046   auto charge = std::make_unique<std::vector<unsigned int>>();
0047   auto path = std::make_unique<std::vector<double>>();
0048 #ifdef ExtendedCALIBTree
0049   auto chargeoverpath = std::make_unique<std::vector<double>>();
0050 #endif
0051   auto amplitude = std::make_unique<std::vector<unsigned char>>();
0052   auto gainused = std::make_unique<std::vector<double>>();
0053   auto gainusedTick = std::make_unique<std::vector<double>>();
0054 
0055   m_tracker = &iSetup.getData(trackerGeometry_token_);
0056   edm::ESHandle<SiStripGain> gainHandle = iSetup.getHandle(gain_token_);
0057   edm::Handle<edm::View<reco::Track>> tracks;
0058   iEvent.getByToken(tracks_token_, tracks);
0059   edm::Handle<TrajTrackAssociationCollection> associations;
0060   iEvent.getByToken(association_token_, associations);
0061 
0062   for (TrajTrackAssociationCollection::const_iterator association = associations->begin();
0063        association != associations->end();
0064        association++) {
0065     const Trajectory* traj = association->key.get();
0066     const reco::Track* track = association->val.get();
0067 
0068     vector<TrajectoryMeasurement> measurements = traj->measurements();
0069     for (vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin();
0070          measurement_it != measurements.end();
0071          measurement_it++) {
0072       TrajectoryStateOnSurface trajState = measurement_it->updatedState();
0073       if (!trajState.isValid())
0074         continue;
0075 
0076       const TrackingRecHit* hit = (*measurement_it->recHit()).hit();
0077       const SiStripRecHit1D* sistripsimple1dhit = dynamic_cast<const SiStripRecHit1D*>(hit);
0078       const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(hit);
0079       const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
0080       const SiPixelRecHit* sipixelhit = dynamic_cast<const SiPixelRecHit*>(hit);
0081 
0082       const SiPixelCluster* PixelCluster = nullptr;
0083       const SiStripCluster* StripCluster = nullptr;
0084       uint32_t DetId = 0;
0085 
0086       for (unsigned int h = 0; h < 2; h++) {
0087         if (!sistripmatchedhit && h == 1) {
0088           continue;
0089         } else if (sistripmatchedhit && h == 0) {
0090           StripCluster = &sistripmatchedhit->monoCluster();
0091           DetId = sistripmatchedhit->monoId();
0092         } else if (sistripmatchedhit && h == 1) {
0093           StripCluster = &sistripmatchedhit->stereoCluster();
0094           ;
0095           DetId = sistripmatchedhit->stereoId();
0096         } else if (sistripsimplehit) {
0097           StripCluster = (sistripsimplehit->cluster()).get();
0098           DetId = sistripsimplehit->geographicalId().rawId();
0099         } else if (sistripsimple1dhit) {
0100           StripCluster = (sistripsimple1dhit->cluster()).get();
0101           DetId = sistripsimple1dhit->geographicalId().rawId();
0102         } else if (sipixelhit) {
0103           PixelCluster = (sipixelhit->cluster()).get();
0104           DetId = sipixelhit->geographicalId().rawId();
0105         } else {
0106           continue;
0107         }
0108 
0109         LocalVector trackDirection = trajState.localDirection();
0110         double cosine = trackDirection.z() / trackDirection.mag();
0111         bool Saturation = false;
0112         bool Overlapping = false;
0113         unsigned int Charge = 0;
0114         double Path = (10.0 * thickness(DetId)) / fabs(cosine);
0115         double PrevGain = -1;
0116         double PrevGainTick = -1;
0117         int FirstStrip = 0;
0118         int NStrips = 0;
0119 
0120         if (StripCluster) {
0121           const auto& Ampls = StripCluster->amplitudes();
0122           FirstStrip = StripCluster->firstStrip();
0123           NStrips = Ampls.size();
0124           int APVId = FirstStrip / 128;
0125 
0126           if (gainHandle.isValid()) {
0127             PrevGain = gainHandle->getApvGain(APVId, gainHandle->getRange(DetId, 1), 1);
0128             PrevGainTick = gainHandle->getApvGain(APVId, gainHandle->getRange(DetId, 0), 1);
0129           }
0130 
0131           for (unsigned int a = 0; a < Ampls.size(); a++) {
0132             Charge += Ampls[a];
0133             if (Ampls[a] >= 254)
0134               Saturation = true;
0135             amplitude->push_back(Ampls[a]);
0136           }
0137 
0138           if (FirstStrip == 0)
0139             Overlapping = true;
0140           if (FirstStrip == 128)
0141             Overlapping = true;
0142           if (FirstStrip == 256)
0143             Overlapping = true;
0144           if (FirstStrip == 384)
0145             Overlapping = true;
0146           if (FirstStrip == 512)
0147             Overlapping = true;
0148           if (FirstStrip == 640)
0149             Overlapping = true;
0150 
0151           if (FirstStrip <= 127 && FirstStrip + Ampls.size() > 127)
0152             Overlapping = true;
0153           if (FirstStrip <= 255 && FirstStrip + Ampls.size() > 255)
0154             Overlapping = true;
0155           if (FirstStrip <= 383 && FirstStrip + Ampls.size() > 383)
0156             Overlapping = true;
0157           if (FirstStrip <= 511 && FirstStrip + Ampls.size() > 511)
0158             Overlapping = true;
0159           if (FirstStrip <= 639 && FirstStrip + Ampls.size() > 639)
0160             Overlapping = true;
0161 
0162           if (FirstStrip + Ampls.size() == 127)
0163             Overlapping = true;
0164           if (FirstStrip + Ampls.size() == 255)
0165             Overlapping = true;
0166           if (FirstStrip + Ampls.size() == 383)
0167             Overlapping = true;
0168           if (FirstStrip + Ampls.size() == 511)
0169             Overlapping = true;
0170           if (FirstStrip + Ampls.size() == 639)
0171             Overlapping = true;
0172           if (FirstStrip + Ampls.size() == 767)
0173             Overlapping = true;
0174         } else if (PixelCluster) {
0175           const auto& Ampls = PixelCluster->pixelADC();
0176           int FirstRow = PixelCluster->minPixelRow();
0177           int FirstCol = PixelCluster->minPixelCol();
0178           FirstStrip = ((FirstRow / 80) << 3 | (FirstCol / 52)) * 128;  //Hack to save the APVId
0179           NStrips = 0;
0180           Saturation = false;
0181           Overlapping = false;
0182 
0183           for (unsigned int a = 0; a < Ampls.size(); a++) {
0184             Charge += Ampls[a];
0185             if (Ampls[a] >= 254)
0186               Saturation = true;
0187           }
0188         }
0189 #ifdef ExtendedCALIBTree
0190         double ChargeOverPath = (double)Charge / Path;
0191 #endif
0192 
0193         trackindex->push_back(shallow::findTrackIndex(tracks, track));
0194         rawid->push_back(DetId);
0195         localdirx->push_back(trackDirection.x());
0196         localdiry->push_back(trackDirection.y());
0197         localdirz->push_back(trackDirection.z());
0198         firststrip->push_back(FirstStrip);
0199         nstrips->push_back(NStrips);
0200         saturation->push_back(Saturation);
0201         overlapping->push_back(Overlapping);
0202         farfromedge->push_back(StripCluster ? isFarFromBorder(&trajState, DetId, &iSetup) : true);
0203         charge->push_back(Charge);
0204         path->push_back(Path);
0205 #ifdef ExtendedCALIBTree
0206         chargeoverpath->push_back(ChargeOverPath);
0207 #endif
0208         gainused->push_back(PrevGain);
0209         gainusedTick->push_back(PrevGainTick);
0210       }
0211     }
0212   }
0213 
0214   iEvent.put(std::move(trackindex), Prefix + "trackindex" + Suffix);
0215   iEvent.put(std::move(rawid), Prefix + "rawid" + Suffix);
0216   iEvent.put(std::move(localdirx), Prefix + "localdirx" + Suffix);
0217   iEvent.put(std::move(localdiry), Prefix + "localdiry" + Suffix);
0218   iEvent.put(std::move(localdirz), Prefix + "localdirz" + Suffix);
0219   iEvent.put(std::move(firststrip), Prefix + "firststrip" + Suffix);
0220   iEvent.put(std::move(nstrips), Prefix + "nstrips" + Suffix);
0221   iEvent.put(std::move(saturation), Prefix + "saturation" + Suffix);
0222   iEvent.put(std::move(overlapping), Prefix + "overlapping" + Suffix);
0223   iEvent.put(std::move(farfromedge), Prefix + "farfromedge" + Suffix);
0224   iEvent.put(std::move(charge), Prefix + "charge" + Suffix);
0225   iEvent.put(std::move(path), Prefix + "path" + Suffix);
0226 #ifdef ExtendedCALIBTree
0227   iEvent.put(std::move(chargeoverpath), Prefix + "chargeoverpath" + Suffix);
0228 #endif
0229   iEvent.put(std::move(amplitude), Prefix + "amplitude" + Suffix);
0230   iEvent.put(std::move(gainused), Prefix + "gainused" + Suffix);
0231   iEvent.put(std::move(gainusedTick), Prefix + "gainusedTick" + Suffix);
0232 }
0233 
0234 bool ShallowGainCalibration::isFarFromBorder(TrajectoryStateOnSurface* trajState,
0235                                              const uint32_t detid,
0236                                              const edm::EventSetup* iSetup) {
0237   edm::ESHandle<TrackerGeometry> tkGeom = iSetup->getHandle(tkGeom_token_);
0238 
0239   LocalPoint HitLocalPos = trajState->localPosition();
0240   LocalError HitLocalError = trajState->localError().positionError();
0241 
0242   const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
0243   if (dynamic_cast<const StripGeomDetUnit*>(it) == nullptr && dynamic_cast<const PixelGeomDetUnit*>(it) == nullptr) {
0244     throw cms::Exception("Logic Error") << "\t\t this detID doesn't seem to belong to the Tracker";
0245   }
0246 
0247   const BoundPlane plane = it->surface();
0248   const TrapezoidalPlaneBounds* trapezoidalBounds(dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
0249   const RectangularPlaneBounds* rectangularBounds(dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
0250 
0251   double DistFromBorder = 1.0;
0252   double HalfLength = it->surface().bounds().length() / 2.0;
0253 
0254   if (trapezoidalBounds) {
0255     std::array<const float, 4> const& parameters = (*trapezoidalBounds).parameters();
0256     HalfLength = parameters[3];
0257   } else if (rectangularBounds) {
0258     HalfLength = it->surface().bounds().length() / 2.0;
0259   } else {
0260     return false;
0261   }
0262 
0263   if (fabs(HitLocalPos.y()) + HitLocalError.yy() >= (HalfLength - DistFromBorder))
0264     return false;
0265 
0266   return true;
0267 }
0268 
0269 double ShallowGainCalibration::thickness(DetId id) {
0270   map<DetId, double>::iterator th = m_thicknessMap.find(id);
0271   if (th != m_thicknessMap.end())
0272     return (*th).second;
0273   else {
0274     double detThickness = 1.;
0275     //compute thickness normalization
0276     const GeomDetUnit* it = m_tracker->idToDetUnit(DetId(id));
0277     bool isPixel = dynamic_cast<const PixelGeomDetUnit*>(it) != nullptr;
0278     bool isStrip = dynamic_cast<const StripGeomDetUnit*>(it) != nullptr;
0279     if (!isPixel && !isStrip) {
0280       throw cms::Exception("Logic Error") << "\t\t this detID doesn't seem to belong to the Tracker";
0281     } else {
0282       detThickness = it->surface().bounds().thickness();
0283     }
0284 
0285     m_thicknessMap[id] = detThickness;  //computed value
0286     return detThickness;
0287   }
0288 }