SingleTrackPattern/src/CosmicTrajectoryBuilder.cc

0001 //
0002 // Package:         RecoTracker/SingleTrackPattern
0003 // Class:           CosmicTrajectoryBuilder
0004 // Original Author:  Michele Pioppi-INFN perugia
0005 #include <vector>
0006 #include <iostream>
0007 #include <cmath>
0008
0009 #include "RecoTracker/SingleTrackPattern/interface/CosmicTrajectoryBuilder.h"
0010 #include "DataFormats/Common/interface/OwnVector.h"
0011 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0012 #include "TrackingTools/TrajectoryState/interface/BasicSingleTrajectoryState.h"
0013 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0018 #include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h"
0019 using namespace std;
0020 CosmicTrajectoryBuilder::CosmicTrajectoryBuilder(const edm::ParameterSet &conf, edm::ConsumesCollector iC)
0021     : magfieldToken_(iC.esConsumes()),
0022       trackerToken_(iC.esConsumes()),
0023       builderToken_(iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("TTRHBuilder")))) {
0024   //minimum number of hits per tracks
0025
0026   theMinHits = conf.getParameter<int>("MinHits");
0027   //cut on chi2
0028   chi2cut = conf.getParameter<double>("Chi2Cut");
0029   edm::LogInfo("CosmicTrackFinder") << "Minimum number of hits " << theMinHits << " Cut on Chi2= " << chi2cut;
0030
0031   geometry = conf.getUntrackedParameter<std::string>("GeometricStructure", "STANDARD");
0032 }
0033
0034 CosmicTrajectoryBuilder::~CosmicTrajectoryBuilder() {}
0035
0036 void CosmicTrajectoryBuilder::init(const edm::EventSetup &es, bool seedplus) {
0037   // FIXME: this is a memory leak generator
0038
0039   //services
0040   magfield = &es.getData(magfieldToken_);
0041   tracker = &es.getData(trackerToken_);
0042
0043   if (seedplus) {
0044     seed_plus = true;
0045     thePropagator = new PropagatorWithMaterial(alongMomentum, 0.1057, &(*magfield));
0046     thePropagatorOp = new PropagatorWithMaterial(oppositeToMomentum, 0.1057, &(*magfield));
0047   } else {
0048     seed_plus = false;
0049     thePropagator = new PropagatorWithMaterial(oppositeToMomentum, 0.1057, &(*magfield));
0050     thePropagatorOp = new PropagatorWithMaterial(alongMomentum, 0.1057, &(*magfield));
0051   }
0052
0053   theUpdator = new KFUpdator();
0054   theEstimator = new Chi2MeasurementEstimator(chi2cut);
0055
0056   RHBuilder = &es.getData(builderToken_);
0057   hitCloner = static_cast<TkTransientTrackingRecHitBuilder const *>(RHBuilder)->cloner();
0058
0059   theFitter = new KFTrajectoryFitter(*thePropagator, *theUpdator, *theEstimator);
0060   theFitter->setHitCloner(&hitCloner);
0061
0062   theSmoother = new KFTrajectorySmoother(*thePropagatorOp, *theUpdator, *theEstimator);
0063   theSmoother->setHitCloner(&hitCloner);
0064 }
0065
0066 void CosmicTrajectoryBuilder::run(const TrajectorySeedCollection &collseed,
0067                                   const SiStripRecHit2DCollection &collstereo,
0068                                   const SiStripRecHit2DCollection &collrphi,
0069                                   const SiStripMatchedRecHit2DCollection &collmatched,
0070                                   const SiPixelRecHitCollection &collpixel,
0071                                   const edm::EventSetup &es,
0072                                   edm::Event &e,
0073                                   vector<Trajectory> &trajoutput) {
0074   std::vector<Trajectory> trajSmooth;
0075   std::vector<Trajectory>::iterator trajIter;
0076
0077   TrajectorySeedCollection::const_iterator iseed;
0078   unsigned int IS = 0;
0079   for (iseed = collseed.begin(); iseed != collseed.end(); iseed++) {
0080     bool seedplus = ((*iseed).direction() == alongMomentum);
0081     init(es, seedplus);
0082     hits.clear();
0083     trajFit.clear();
0084     trajSmooth.clear();
0085     //order all the hits
0086     vector<const TrackingRecHit *> allHits = SortHits(collstereo, collrphi, collmatched, collpixel, *iseed);
0087     Trajectory startingTraj = createStartingTrajectory(*iseed);
0088     AddHit(startingTraj, allHits);
0089     for (trajIter = trajFit.begin(); trajIter != trajFit.end(); trajIter++) {
0090       trajSmooth = theSmoother->trajectories((*trajIter));
0091     }
0092     for (trajIter = trajSmooth.begin(); trajIter != trajSmooth.end(); trajIter++) {
0093       if ((*trajIter).isValid()) {
0094         trajoutput.push_back((*trajIter));
0095       }
0096     }
0097     delete theUpdator;
0098     delete theEstimator;
0099     delete thePropagator;
0100     delete thePropagatorOp;
0101     delete theFitter;
0102     delete theSmoother;
0103     //Only the first 30 seeds are considered
0104     if (IS > 30)
0105       return;
0106     IS++;
0107   }
0108 }
0109
0110 Trajectory CosmicTrajectoryBuilder::createStartingTrajectory(const TrajectorySeed &seed) const {
0111   Trajectory result(seed, seed.direction());
0112   std::vector<TM> &&seedMeas = seedMeasurements(seed);
0113   for (auto &i : seedMeas)
0114     result.push(std::move(i));
0115   return result;
0116 }
0117
0118 std::vector<TrajectoryMeasurement> CosmicTrajectoryBuilder::seedMeasurements(const TrajectorySeed &seed) const {
0119   std::vector<TrajectoryMeasurement> result;
0120   auto const &hitRange = seed.recHits();
0121   for (auto ihit = hitRange.begin(); ihit != hitRange.end(); ihit++) {
0122     //RC TransientTrackingRecHit* recHit = RHBuilder->build(&(*ihit));
0123     TransientTrackingRecHit::RecHitPointer recHit = RHBuilder->build(&(*ihit));
0124     const GeomDet *hitGeomDet = (&(*tracker))->idToDet(ihit->geographicalId());
0125     TSOS invalidState(new BasicSingleTrajectoryState(hitGeomDet->surface()));
0126
0127     if (ihit == hitRange.end() - 1) {
0128       TSOS updatedState = startingTSOS(seed);
0129       result.emplace_back(invalidState, updatedState, recHit);
0130     } else {
0131       result.emplace_back(invalidState, recHit);
0132     }
0133   }
0134
0135   return result;
0136 }
0137
0138 vector<const TrackingRecHit *> CosmicTrajectoryBuilder::SortHits(const SiStripRecHit2DCollection &collstereo,
0139                                                                  const SiStripRecHit2DCollection &collrphi,
0140                                                                  const SiStripMatchedRecHit2DCollection &collmatched,
0141                                                                  const SiPixelRecHitCollection &collpixel,
0142                                                                  const TrajectorySeed &seed) {
0143   //The Hits with global y more than the seed are discarded
0144   //The Hits correspondign to the seed are discarded
0145   //At the end all the hits are sorted in y
0146   vector<const TrackingRecHit *> allHits;
0147
0148   SiStripRecHit2DCollection::DataContainer::const_iterator istrip;
0149   float yref = 0.;
0150   for (auto const &recHit : seed.recHits()) {
0151     yref = RHBuilder->build(&recHit)->globalPosition().y();
0152     hits.push_back((RHBuilder->build(&recHit)));
0153     LogDebug("CosmicTrackFinder") << "SEED HITS" << RHBuilder->build(&recHit)->globalPosition();
0154   }
0155
0156   SiPixelRecHitCollection::DataContainer::const_iterator ipix;
0157   for (ipix = collpixel.data().begin(); ipix != collpixel.data().end(); ipix++) {
0158     float ych = RHBuilder->build(&(*ipix))->globalPosition().y();
0159     if ((seed_plus && (ych < yref)) || (!(seed_plus) && (ych > yref)))
0160       allHits.push_back(&(*ipix));
0161   }
0162
0163   for (istrip = collrphi.data().begin(); istrip != collrphi.data().end(); istrip++) {
0164     float ych = RHBuilder->build(&(*istrip))->globalPosition().y();
0165     if ((seed_plus && (ych < yref)) || (!(seed_plus) && (ych > yref)))
0166       allHits.push_back(&(*istrip));
0167   }
0168
0169   for (istrip = collstereo.data().begin(); istrip != collstereo.data().end(); istrip++) {
0170     float ych = RHBuilder->build(&(*istrip))->globalPosition().y();
0171     if ((seed_plus && (ych < yref)) || (!(seed_plus) && (ych > yref)))
0172       allHits.push_back(&(*istrip));
0173   }
0174
0175   if (seed_plus)
0176     stable_sort(allHits.begin(), allHits.end(), CompareHitY_plus(*tracker));
0177   else
0178     stable_sort(allHits.begin(), allHits.end(), CompareHitY(*tracker));
0179
0180   return allHits;
0181 }
0182
0183 TrajectoryStateOnSurface CosmicTrajectoryBuilder::startingTSOS(const TrajectorySeed &seed) const {
0184   PTrajectoryStateOnDet pState(seed.startingState());
0185   const GeomDet *gdet = (&(*tracker))->idToDet(DetId(pState.detId()));
0186   TSOS State = trajectoryStateTransform::transientState(pState, &(gdet->surface()), &(*magfield));
0187   return State;
0188 }
0189
0190 void CosmicTrajectoryBuilder::AddHit(Trajectory &traj, const vector<const TrackingRecHit *> &Hits) {
0191   unsigned int icosm2;
0192   unsigned int ibestdet;
0193   float chi2min;
0194   for (unsigned int icosmhit = 0; icosmhit < Hits.size(); icosmhit++) {
0195     GlobalPoint gphit = RHBuilder->build(Hits[icosmhit])->globalPosition();
0196     unsigned int iraw = Hits[icosmhit]->geographicalId().rawId();
0197     LogDebug("CosmicTrackFinder") << " HIT POSITION " << gphit;
0198     //RC TransientTrackingRecHit* tmphit=RHBuilder->build(Hits[icosmhit]);
0199     TransientTrackingRecHit::RecHitPointer tmphit = RHBuilder->build(Hits[icosmhit]);
0200     TSOS prSt = thePropagator->propagate(traj.lastMeasurement().updatedState(),
0201                                          tracker->idToDet(Hits[icosmhit]->geographicalId())->surface());
0202
0203     //After propagating the trajectory to a detector,
0204     //find the most compatible hit in the det
0205     chi2min = 20000000;
0206     ibestdet = 1000;
0207     if (prSt.isValid()) {
0208       LocalPoint prLoc = prSt.localPosition();
0209       LogDebug("CosmicTrackFinder") << "STATE PROPAGATED AT DET " << iraw << " " << prSt;
0210       for (icosm2 = icosmhit; icosm2 < Hits.size(); icosm2++) {
0211         if (iraw == Hits[icosm2]->geographicalId().rawId()) {
0212           TransientTrackingRecHit::RecHitPointer tmphit = RHBuilder->build(Hits[icosm2]);
0213           float contr = theEstimator->estimate(prSt, *tmphit).second;
0214           if (contr < chi2min) {
0215             chi2min = contr;
0216             ibestdet = icosm2;
0217           }
0218           if (icosm2 != icosmhit)
0219             icosmhit++;
0220
0221         } else
0222           icosm2 = Hits.size();
0223       }
0224
0225       if (chi2min < chi2cut)
0226         LogDebug("CosmicTrackFinder") << "Chi2 contribution for hit at "
0227                                       << RHBuilder->build(Hits[ibestdet])->globalPosition() << " is " << chi2min;
0228
0229       if (traj.foundHits() < 3 && (chi2min < chi2cut)) {
0230         //check on the first hit after the seed
0231         GlobalVector ck = RHBuilder->build(Hits[ibestdet])->globalPosition() -
0232                           traj.firstMeasurement().updatedState().globalPosition();
0233         if ((abs(ck.x() / ck.y()) > 2) || (abs(ck.z() / ck.y()) > 2))
0234           chi2min = 300;
0235       }
0236       if (chi2min < chi2cut) {
0237         if (abs(prLoc.x()) < 25 && abs(prLoc.y()) < 25) {
0238           TransientTrackingRecHit::RecHitPointer tmphitbestdet = RHBuilder->build(Hits[ibestdet]);
0239           TSOS UpdatedState = theUpdator->update(prSt, *tmphitbestdet);
0240           if (UpdatedState.isValid()) {
0241             traj.push(TM(prSt, UpdatedState, RHBuilder->build(Hits[ibestdet]), chi2min));
0242             LogDebug("CosmicTrackFinder") << "STATE UPDATED WITH HIT AT POSITION " << tmphitbestdet->globalPosition()
0243                                           << UpdatedState << " " << traj.chiSquared();
0244
0245             hits.push_back(tmphitbestdet);
0246           }
0247         } else
0248           LogDebug("CosmicTrackFinder") << " Hits outside module surface " << prLoc;
0249       } else
0250         LogDebug("CosmicTrackFinder") << " State can not be updated with hit at position " << gphit;
0251     } else
0252       LogDebug("CosmicTrackFinder") << " State can not be propagated at det " << iraw;
0253   }
0254
0255   if (qualityFilter(traj)) {
0256     const TrajectorySeed &tmpseed = traj.seed();
0257     if (thePropagatorOp
0258             ->propagate(traj.lastMeasurement().updatedState(),
0259                         tracker->idToDet((*hits.begin())->geographicalId())->surface())
0260             .isValid()) {
0261       TSOS startingState = thePropagatorOp->propagate(traj.lastMeasurement().updatedState(),
0262                                                       tracker->idToDet((*hits.begin())->geographicalId())->surface());
0263
0264       trajFit = theFitter->fit(tmpseed, hits, startingState);
0265     }
0266   }
0267 }
0268
0269 bool CosmicTrajectoryBuilder::qualityFilter(const Trajectory &traj) {
0270   int ngoodhits = 0;
0271   if (geometry == "MTCC") {
0272     auto hits = traj.recHits();
0273     for (auto hit = hits.begin(); hit != hits.end(); hit++) {
0274       unsigned int iid = (*hit)->hit()->geographicalId().rawId();
0275       //CHECK FOR 3 hits r-phi
0276       if (((iid >> 0) & 0x3) != 1)
0277         ngoodhits++;
0278     }
0279   } else
0280     ngoodhits = traj.foundHits();
0281
0282   if (ngoodhits >= theMinHits) {
0283     return true;
0284   } else {
0285     return false;
0286   }
0287 }