Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:28:43 

0001 //
0002 // Package:         RecoTracker/SingleTrackPattern
0003 // Class:           CRackTrajectoryBuilder
0004 // Original Author:  Michele Pioppi-INFN perugia
0005 //
0006 // Package:         RecoTracker/SingleTrackPattern
0007 // Class:           CRackTrajectoryBuilder
0008 // Original Author:  Michele Pioppi-INFN perugia
0009 
0010 #include <vector>
0011 #include <iostream>
0012 #include <cmath>
0013 
0014 #include "RecoTracker/SingleTrackPattern/interface/CRackTrajectoryBuilder.h"
0015 #include "DataFormats/Common/interface/OwnVector.h"
0016 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0017 #include "TrackingTools/TrajectoryState/interface/BasicSingleTrajectoryState.h"
0018 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0023 #include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h"
0024 
0025 //#include "RecoTracker/CkfPattern/interface/TransientInitialStateEstimator.h"
0026 
0027 using namespace std;
0028 
0029 CRackTrajectoryBuilder::CRackTrajectoryBuilder(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0030     : magfieldToken_(iC.esConsumes()),
0031       trackerToken_(iC.esConsumes()),
0032       builderToken_(iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("TTRHBuilder")))) {
0033   //minimum number of hits per tracks
0034 
0035   theMinHits = conf.getParameter<int>("MinHits");
0036   //cut on chi2
0037   chi2cut = conf.getParameter<double>("Chi2Cut");
0038   edm::LogInfo("CosmicTrackFinder") << "Minimum number of hits " << theMinHits << " Cut on Chi2= " << chi2cut;
0039 
0040   debug_info = conf.getUntrackedParameter<bool>("debug", false);
0041   fastPropagation = conf.getUntrackedParameter<bool>("fastPropagation", false);
0042   useMatchedHits = conf.getUntrackedParameter<bool>("useMatchedHits", true);
0043 
0044   geometry = conf.getUntrackedParameter<std::string>("GeometricStructure", "STANDARD");
0045 }
0046 
0047 CRackTrajectoryBuilder::~CRackTrajectoryBuilder() {
0048   //  delete theInitialState;
0049 }
0050 
0051 void CRackTrajectoryBuilder::init(const edm::EventSetup& es, bool seedplus) {
0052   //  edm::ParameterSet tise_params = conf_.getParameter<edm::ParameterSet>("TransientInitialStateEstimatorParameters") ;
0053   // theInitialState          = new TransientInitialStateEstimator( es,tise_params);
0054 
0055   //services
0056   magfield = &es.getData(magfieldToken_);
0057   tracker = &es.getData(trackerToken_);
0058 
0059   if (seedplus) {
0060     seed_plus = true;
0061     thePropagator = new PropagatorWithMaterial(alongMomentum, 0.1057, &(*magfield));
0062     thePropagatorOp = new PropagatorWithMaterial(oppositeToMomentum, 0.1057, &(*magfield));
0063   } else {
0064     seed_plus = false;
0065     thePropagator = new PropagatorWithMaterial(oppositeToMomentum, 0.1057, &(*magfield));
0066     thePropagatorOp = new PropagatorWithMaterial(alongMomentum, 0.1057, &(*magfield));
0067   }
0068 
0069   theUpdator = new KFUpdator();
0070   //  theUpdator=       new KFStripUpdator();
0071   theEstimator = new Chi2MeasurementEstimator(chi2cut);
0072 
0073   RHBuilder = &es.getData(builderToken_);
0074 
0075   theFitter = new KFTrajectoryFitter(*thePropagator, *theUpdator, *theEstimator);
0076 
0077   theSmoother = new KFTrajectorySmoother(*thePropagatorOp, *theUpdator, *theEstimator);
0078 }
0079 
0080 void CRackTrajectoryBuilder::run(const TrajectorySeedCollection& collseed,
0081                                  const SiStripRecHit2DCollection& collstereo,
0082                                  const SiStripRecHit2DCollection& collrphi,
0083                                  const SiStripMatchedRecHit2DCollection& collmatched,
0084                                  const SiPixelRecHitCollection& collpixel,
0085                                  const edm::EventSetup& es,
0086                                  edm::Event& e,
0087                                  vector<Trajectory>& trajoutput) {
0088   std::vector<Trajectory> trajSmooth;
0089   std::vector<Trajectory>::iterator trajIter;
0090 
0091   TrajectorySeedCollection::const_iterator iseed;
0092   unsigned int IS = 0;
0093   for (iseed = collseed.begin(); iseed != collseed.end(); iseed++) {
0094     bool seedplus = ((*iseed).direction() == alongMomentum);
0095     init(es, seedplus);
0096     hits.clear();
0097     trajFit.clear();
0098     trajSmooth.clear();
0099 
0100     Trajectory startingTraj = createStartingTrajectory(*iseed);
0101     Trajectory trajTmp;  // used for the opposite direction
0102     Trajectory traj = startingTraj;
0103 
0104     //first propagate track in opposite direction ...
0105     seed_plus = !seed_plus;
0106     vector<const TrackingRecHit*> allHitsOppsite = SortHits(collstereo, collrphi, collmatched, collpixel, *iseed, true);
0107     seed_plus = !seed_plus;
0108     if (!allHitsOppsite.empty()) {
0109       //there are hits which are above the seed,
0110       //cout << "Number of hits higher than seed " <<allHitsOppsite.size() << endl;
0111 
0112       AddHit(traj, allHitsOppsite, thePropagatorOp);
0113 
0114       if (debug_info) {
0115         cout << "Hits in opposite direction..." << endl;
0116         TransientTrackingRecHit::RecHitContainer::const_iterator iHit;
0117         for (iHit = hits.begin(); iHit != hits.end(); iHit++) {
0118           cout << (**iHit).globalPosition() << endl;
0119         }
0120       }
0121 
0122       // now add hist opposite to seed
0123       //now crate new starting trajectory
0124       reverse(hits.begin(), hits.end());
0125       if (thePropagatorOp
0126               ->propagate(traj.firstMeasurement().updatedState(),
0127                           tracker->idToDet((hits.front())->geographicalId())->surface())
0128               .isValid()) {
0129         TSOS startingStateNew =  //TrajectoryStateWithArbitraryError()//
0130             (thePropagatorOp->propagate(traj.firstMeasurement().updatedState(),
0131                                         tracker->idToDet((hits.front())->geographicalId())->surface()));
0132 
0133         if (debug_info) {
0134           cout << "Hits in opposite direction reversed..." << endl;
0135           TransientTrackingRecHit::RecHitContainer::const_iterator iHit;
0136           for (iHit = hits.begin(); iHit != hits.end(); iHit++) {
0137             cout << (**iHit).globalPosition() << endl;
0138           }
0139         }
0140 
0141         const TrajectorySeed& tmpseed = traj.seed();
0142         trajTmp = theFitter->fit(tmpseed, hits, startingStateNew).front();
0143 
0144         if (debug_info) {
0145           cout << "Debugging show fitted hits" << endl;
0146           auto hitsFit = trajTmp.recHits();
0147           for (auto hit = hitsFit.begin(); hit != hitsFit.end(); hit++) {
0148             cout << RHBuilder->build(&(*(*hit)->hit()))->globalPosition() << endl;
0149           }
0150         }
0151       }
0152     } else {
0153       if (debug_info)
0154         cout << "There are no hits in opposite direction ..." << endl;
0155     }
0156 
0157     vector<const TrackingRecHit*> allHits;
0158     if (trajTmp.foundHits()) {
0159       traj = trajTmp;
0160       allHits = SortHits(collstereo, collrphi, collmatched, collpixel, *iseed, false);
0161     } else {
0162       traj = startingTraj;
0163       hits.clear();
0164       allHits = SortHits(collstereo, collrphi, collmatched, collpixel, *iseed, true);
0165     }
0166 
0167     AddHit(traj, allHits, thePropagator);
0168 
0169     if (debug_info) {
0170       cout << "Debugging show All fitted hits" << endl;
0171       auto hits = traj.recHits();
0172       for (auto hit = hits.begin(); hit != hits.end(); hit++) {
0173         cout << (*hit)->globalPosition() << endl;
0174       }
0175 
0176       cout << qualityFilter(traj) << " <- quality filter good?" << endl;
0177     }
0178 
0179     if (debug_info)
0180       cout << "now do quality checks" << endl;
0181     if (qualityFilter(traj)) {
0182       const TrajectorySeed& tmpseed = traj.seed();
0183       std::pair<TrajectoryStateOnSurface, const GeomDet*> initState =
0184           innerState(traj);  //theInitialState->innerState(traj);
0185       if (initState.first.isValid())
0186         trajFit = theFitter->fit(tmpseed, hits, initState.first);
0187     }
0188 
0189     for (trajIter = trajFit.begin(); trajIter != trajFit.end(); trajIter++) {
0190       trajSmooth = theSmoother->trajectories((*trajIter));
0191     }
0192     for (trajIter = trajSmooth.begin(); trajIter != trajSmooth.end(); trajIter++) {
0193       if ((*trajIter).isValid()) {
0194         if (debug_info)
0195           cout << "adding track ... " << endl;
0196         trajoutput.push_back((*trajIter));
0197       }
0198     }
0199     delete theUpdator;
0200     delete theEstimator;
0201     delete thePropagator;
0202     delete thePropagatorOp;
0203     delete theFitter;
0204     delete theSmoother;
0205     //Only the first 30 seeds are considered
0206     if (IS > 30)
0207       return;
0208     IS++;
0209   }
0210 }
0211 
0212 Trajectory CRackTrajectoryBuilder::createStartingTrajectory(const TrajectorySeed& seed) const {
0213   Trajectory result(seed, seed.direction());
0214   std::vector<TM>&& seedMeas = seedMeasurements(seed);
0215   for (auto i : seedMeas)
0216     result.push(std::move(i));
0217   return result;
0218 }
0219 
0220 std::vector<TrajectoryMeasurement> CRackTrajectoryBuilder::seedMeasurements(const TrajectorySeed& seed) const {
0221   std::vector<TrajectoryMeasurement> result;
0222   auto const& hitRange = seed.recHits();
0223   for (auto ihit = hitRange.begin(); ihit != hitRange.end(); ihit++) {
0224     //RC TransientTrackingRecHit* recHit = RHBuilder->build(&(*ihit));
0225     TransientTrackingRecHit::RecHitPointer recHit = RHBuilder->build(&(*ihit));
0226     const GeomDet* hitGeomDet = (&(*tracker))->idToDet(ihit->geographicalId());
0227     TSOS invalidState(new BasicSingleTrajectoryState(hitGeomDet->surface()));
0228 
0229     if (ihit == hitRange.end() - 1) {
0230       TSOS updatedState = startingTSOS(seed);
0231       result.emplace_back(invalidState, updatedState, recHit);
0232 
0233     } else {
0234       result.emplace_back(invalidState, recHit);
0235     }
0236   }
0237 
0238   return result;
0239 }
0240 
0241 vector<const TrackingRecHit*> CRackTrajectoryBuilder::SortHits(const SiStripRecHit2DCollection& collstereo,
0242                                                                const SiStripRecHit2DCollection& collrphi,
0243                                                                const SiStripMatchedRecHit2DCollection& collmatched,
0244                                                                const SiPixelRecHitCollection& collpixel,
0245                                                                const TrajectorySeed& seed,
0246                                                                const bool bAddSeedHits) {
0247   //The Hits with global y more than the seed are discarded
0248   //The Hits correspondign to the seed are discarded
0249   //At the end all the hits are sorted in y
0250   vector<const TrackingRecHit*> allHits;
0251 
0252   SiStripRecHit2DCollection::DataContainer::const_iterator istrip;
0253   float yref = 0.;
0254 
0255   if (debug_info)
0256     cout << "SEED " << startingTSOS(seed) << endl;
0257   if (debug_info)
0258     cout << "seed hits size " << seed.nHits() << endl;
0259 
0260   //seed debugging:
0261   // GeomDet *seedDet =  tracker->idToDet(seed.startingState().detId());
0262 
0263   //  edm::LogInfo("CRackTrajectoryBuilder::SortHits") << "SEED " << seed.startingState();
0264 
0265   edm::LogInfo("CRackTrajectoryBuilder::SortHits") << "SEED " << startingTSOS(seed);
0266   //  edm::LogInfo("CRackTrajectoryBuilder::SortHits" << "seed hits size " << seed.nHits();
0267 
0268   //  if (seed.nHits()<2)
0269   //  return allHits;
0270 
0271   float_t yMin = 0.;
0272   float_t yMax = 0.;
0273 
0274   int seedHitSize = seed.nHits();
0275 
0276   vector<int> detIDSeedMatched(seedHitSize);
0277   vector<int> detIDSeedRphi(seedHitSize);
0278   vector<int> detIDSeedStereo(seedHitSize);
0279 
0280   auto const& hRange = seed.recHits();
0281   for (auto ihit = hRange.begin(); ihit != hRange.end(); ihit++) {
0282     // need to find track with lowest (seed_plus)/ highest y (seed_minus)
0283     // split matched hits ...
0284     const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(&(*ihit));
0285 
0286     yref = RHBuilder->build(&(*ihit))->globalPosition().y();
0287     if (ihit == hRange.begin()) {
0288       yMin = yref;
0289       yMax = yref;
0290     }
0291 
0292     if (matchedhit) {
0293       auto m = matchedhit->monoHit();
0294       auto s = matchedhit->stereoHit();
0295       float_t yGlobRPhi = RHBuilder->build(&m)->globalPosition().y();
0296       float_t yGlobStereo = RHBuilder->build(&s)->globalPosition().y();
0297 
0298       if (debug_info)
0299         cout << "Rphi ..." << yGlobRPhi << endl;
0300       if (debug_info)
0301         cout << "Stereo ..." << yGlobStereo << endl;
0302 
0303       if (yGlobStereo < yMin)
0304         yMin = yGlobStereo;
0305       if (yGlobRPhi < yMin)
0306         yMin = yGlobRPhi;
0307 
0308       if (yGlobStereo > yMax)
0309         yMax = yGlobStereo;
0310       if (yGlobRPhi > yMax)
0311         yMax = yGlobRPhi;
0312 
0313       detIDSeedMatched.push_back(matchedhit->geographicalId().rawId());
0314       detIDSeedRphi.push_back(m.geographicalId().rawId());
0315       detIDSeedStereo.push_back(s.geographicalId().rawId());
0316 
0317       if (bAddSeedHits) {
0318         if (useMatchedHits) {
0319           hits.push_back((RHBuilder->build(&(*ihit))));
0320         } else {
0321           if (((yGlobRPhi > yGlobStereo) && seed_plus) || ((yGlobRPhi < yGlobStereo) && !seed_plus)) {
0322             hits.push_back((RHBuilder->build(&m)));
0323             hits.push_back((RHBuilder->build(&s)));
0324           } else {
0325             hits.push_back((RHBuilder->build(&s)));
0326             hits.push_back((RHBuilder->build(&m)));
0327           }
0328         }
0329       }
0330     } else if (bAddSeedHits) {
0331       hits.push_back((RHBuilder->build(&(*ihit))));
0332       detIDSeedRphi.push_back(ihit->geographicalId().rawId());
0333       detIDSeedMatched.push_back(-1);
0334       detIDSeedStereo.push_back(-1);
0335     }
0336 
0337     if (yref < yMin)
0338       yMin = yref;
0339     if (yref > yMax)
0340       yMax = yref;
0341 
0342     //    if (bAddSeedHits)
0343     //      hits.push_back((RHBuilder->build(&(*ihit))));
0344 
0345     LogDebug("CosmicTrackFinder") << "SEED HITS" << RHBuilder->build(&(*ihit))->globalPosition();
0346     if (debug_info)
0347       cout << "SEED HITS" << RHBuilder->build(&(*ihit))->globalPosition() << endl;
0348     //    if (debug_info) cout <<"SEED HITS"<< seed.startingState().parameters() << endl;
0349   }
0350 
0351   yref = (seed_plus) ? yMin : yMax;
0352 
0353   SiPixelRecHitCollection::DataContainer::const_iterator ipix;
0354   for (ipix = collpixel.data().begin(); ipix != collpixel.data().end(); ipix++) {
0355     float ych = RHBuilder->build(&(*ipix))->globalPosition().y();
0356     if ((seed_plus && (ych < yref)) || (!(seed_plus) && (ych > yref)))
0357       allHits.push_back(&(*ipix));
0358   }
0359 
0360   if (useMatchedHits)  // use matched
0361   {
0362     //add the matched hits ...
0363     SiStripMatchedRecHit2DCollection::DataContainer::const_iterator istripm;
0364 
0365     for (istripm = collmatched.data().begin(); istripm != collmatched.data().end(); istripm++) {
0366       float ych = RHBuilder->build(&(*istripm))->globalPosition().y();
0367 
0368       int cDetId = istripm->geographicalId().rawId();
0369       bool noSeedDet = (detIDSeedMatched.end() == find(detIDSeedMatched.begin(), detIDSeedMatched.end(), cDetId));
0370 
0371       if (noSeedDet)
0372         if ((seed_plus && (ych < yref)) || (!(seed_plus) && (ych > yref)))
0373           allHits.push_back(&(*istripm));
0374     }
0375 
0376     //add the rpi hits, but only accept hits that are not matched hits
0377     for (istrip = collrphi.data().begin(); istrip != collrphi.data().end(); istrip++) {
0378       float ych = RHBuilder->build(&(*istrip))->globalPosition().y();
0379       StripSubdetector monoDetId(istrip->geographicalId());
0380       if (monoDetId.partnerDetId()) {
0381         edm::LogInfo("CRackTrajectoryBuilder::SortHits") << "this det belongs to a glued det " << ych << endl;
0382         continue;
0383       }
0384       int cDetId = istrip->geographicalId().rawId();
0385       bool noSeedDet = (detIDSeedRphi.end() == find(detIDSeedRphi.begin(), detIDSeedRphi.end(), cDetId));
0386       if (noSeedDet)
0387         if ((seed_plus && (ych < yref)) || (!(seed_plus) && (ych > yref))) {
0388           bool hitIsUnique = true;
0389           //now
0390           for (istripm = collmatched.data().begin(); istripm != collmatched.data().end(); istripm++) {
0391             //       if ( isDifferentStripReHit2D ( *istrip, (istripm->stereoHit() ) ) == false)
0392             if (isDifferentStripReHit2D(*istrip, (istripm->monoHit())) == false) {
0393               hitIsUnique = false;
0394               edm::LogInfo("CRackTrajectoryBuilder::SortHits") << "rphi hit is in matched hits; y: " << ych << endl;
0395               break;
0396             }
0397           }  //end loop over all matched
0398           if (hitIsUnique) {
0399             //      if (debug_info) cout << "adding rphi hit " << &(*istrip) << endl;
0400             allHits.push_back(&(*istrip));
0401           }
0402         }
0403     }
0404 
0405   } else  // dont use matched ...
0406   {
0407     for (istrip = collrphi.data().begin(); istrip != collrphi.data().end(); istrip++) {
0408       float ych = RHBuilder->build(&(*istrip))->globalPosition().y();
0409       if ((seed_plus && (ych < yref)) || (!(seed_plus) && (ych > yref)))
0410         allHits.push_back(&(*istrip));
0411     }
0412 
0413     for (istrip = collstereo.data().begin(); istrip != collstereo.data().end(); istrip++) {
0414       float ych = RHBuilder->build(&(*istrip))->globalPosition().y();
0415       if ((seed_plus && (ych < yref)) || (!(seed_plus) && (ych > yref)))
0416         allHits.push_back(&(*istrip));
0417     }
0418   }
0419 
0420   if (seed_plus)
0421     stable_sort(allHits.begin(), allHits.end(), CompareHitY_plus(*tracker));
0422   else
0423     stable_sort(allHits.begin(), allHits.end(), CompareHitY(*tracker));
0424 
0425   if (debug_info) {
0426     if (debug_info)
0427       cout << "all hits" << endl;
0428 
0429     //starting trajectory
0430     Trajectory startingTraj = createStartingTrajectory(seed);
0431 
0432     if (debug_info)
0433       cout << "START " << startingTraj.lastMeasurement().updatedState() << endl;
0434     if (debug_info)
0435       cout << "START  Err" << startingTraj.lastMeasurement().updatedState().localError().matrix() << endl;
0436 
0437     vector<const TrackingRecHit*>::iterator iHit;
0438     for (iHit = allHits.begin(); iHit < allHits.end(); iHit++) {
0439       GlobalPoint gphit = RHBuilder->build(*iHit)->globalPosition();
0440       if (debug_info)
0441         cout << "GH " << gphit << endl;
0442 
0443       //      tracker->idToDet((*iHit)->geographicalId())->surface();
0444 
0445       TSOS prSt = thePropagator->propagate(startingTraj.lastMeasurement().updatedState(),
0446                                            tracker->idToDet((*iHit)->geographicalId())->surface());
0447 
0448       if (prSt.isValid()) {
0449         if (debug_info)
0450           cout << "PR " << prSt.globalPosition() << endl;
0451         //if (debug_info) cout << "PR  Err" << prSt.localError().matrix() << endl;
0452 
0453       } else {
0454         if (debug_info)
0455           cout << "not valid" << endl;
0456       }
0457     }
0458     if (debug_info)
0459       cout << "all hits end" << endl;
0460   }
0461 
0462   return allHits;
0463 }
0464 
0465 TrajectoryStateOnSurface CRackTrajectoryBuilder::startingTSOS(const TrajectorySeed& seed) const {
0466   PTrajectoryStateOnDet pState(seed.startingState());
0467   const GeomDet* gdet = (&(*tracker))->idToDet(DetId(pState.detId()));
0468   TSOS State = trajectoryStateTransform::transientState(pState, &(gdet->surface()), &(*magfield));
0469   return State;
0470 }
0471 
0472 void CRackTrajectoryBuilder::AddHit(Trajectory& traj,
0473                                     const vector<const TrackingRecHit*>& _Hits,
0474                                     Propagator* currPropagator) {
0475   vector<const TrackingRecHit*> Hits = _Hits;
0476   if (Hits.empty())
0477     return;
0478 
0479   if (debug_info)
0480     cout << "CRackTrajectoryBuilder::AddHit" << endl;
0481   if (debug_info)
0482     cout << "START " << traj.lastMeasurement().updatedState() << endl;
0483 
0484   vector<TrackingRecHitRange> hitRangeByDet;
0485   TrackingRecHitIterator prevDet;
0486 
0487   prevDet = Hits.begin();
0488   for (TrackingRecHitIterator iHit = Hits.begin(); iHit != Hits.end(); iHit++) {
0489     if ((*prevDet)->geographicalId() == (*iHit)->geographicalId())
0490       continue;
0491 
0492     hitRangeByDet.push_back(make_pair(prevDet, iHit));
0493     prevDet = iHit;
0494   }
0495   hitRangeByDet.push_back(make_pair(prevDet, Hits.end()));
0496 
0497   /// do the old version ....
0498 
0499   if (fastPropagation) {
0500     for (TrackingRecHitRangeIterator iHitRange = hitRangeByDet.begin(); iHitRange != hitRangeByDet.end(); iHitRange++) {
0501       const TrackingRecHit* currHit = *(iHitRange->first);
0502       DetId currDet = currHit->geographicalId();
0503 
0504       TSOS prSt =
0505           currPropagator->propagate(traj.lastMeasurement().updatedState(), tracker->idToDet(currDet)->surface());
0506 
0507       if (!prSt.isValid()) {
0508         if (debug_info)
0509           cout << "Not Valid: PRST" << prSt.globalPosition();
0510         //        if (debug_info) cout << "Not Valid: HIT" << *currHit;
0511 
0512         continue;
0513       }
0514 
0515       TransientTrackingRecHit::RecHitPointer bestHit = RHBuilder->build(currHit);
0516       double chi2min = theEstimator->estimate(prSt, *bestHit).second;
0517 
0518       if (debug_info)
0519         cout << "Size " << iHitRange->first - (*iHitRange).second << endl;
0520       for (TrackingRecHitIterator iHit = (*iHitRange).first + 1; iHit != iHitRange->second; iHit++) {
0521         if (debug_info)
0522           cout << "loop3 "
0523                << " " << Hits.end() - iHit << endl;
0524 
0525         TransientTrackingRecHit::RecHitPointer tmpHit = RHBuilder->build(*iHit);
0526         double currChi2 = theEstimator->estimate(prSt, *tmpHit).second;
0527         if (currChi2 < chi2min) {
0528           chi2min = currChi2;
0529           bestHit = tmpHit;
0530         }
0531       }
0532       //now we have check if the track can be added to the trajectory
0533       if (debug_info)
0534         cout << chi2min << endl;
0535       if (chi2min < chi2cut) {
0536         if (debug_info)
0537           cout << "chi2 fine : " << chi2min << endl;
0538         TSOS UpdatedState = theUpdator->update(prSt, *bestHit);
0539         if (UpdatedState.isValid()) {
0540           hits.push_back(bestHit);
0541           traj.push(TM(prSt, UpdatedState, bestHit, chi2min));
0542           if (debug_info)
0543             edm::LogInfo("CosmicTrackFinder") << "STATE UPDATED WITH HIT AT POSITION "
0544 
0545                                               << bestHit->globalPosition() << UpdatedState << " " << traj.chiSquared();
0546           if (debug_info)
0547             cout << "STATE UPDATED WITH HIT AT POSITION "
0548 
0549                  << bestHit->globalPosition() << UpdatedState << " " << traj.chiSquared();
0550           if (debug_info)
0551             cout << "State is valid ..." << endl;
0552           break;  // now we need to
0553         } else {
0554           edm::LogWarning("CosmicTrackFinder") << " State can not be updated with hit at position " << endl;
0555           TSOS UpdatedState = theUpdator->update(prSt, *bestHit);
0556           if (UpdatedState.isValid()) {
0557             cout << "NOT! UPDATED WITH HIT AT POSITION "
0558 
0559                  << bestHit->globalPosition() << UpdatedState << " " << traj.chiSquared();
0560           }
0561         }
0562       }
0563     }
0564   }  //simple version end
0565   else {
0566     //first sort the dets in the order they are traversed by the trajectory
0567     // we need three loops:
0568     // 1: loop as long as there can be an new hit added to the trajectory
0569     // 2: loop over all dets that might be hit
0570     // 3: loop over all hits on a certain det
0571 
0572     std::vector<std::pair<TrackingRecHitRangeIterator, TSOS> > trackHitCandidates;
0573     std::vector<std::pair<TrackingRecHitRangeIterator, TSOS> >::iterator iHitRange;
0574     std::vector<uint32_t> processedDets;
0575     do {
0576       //create vector of possibly hit detectors...
0577       trackHitCandidates.clear();
0578       DetId currDet;
0579       for (TrackingRecHitRangeIterator iHit = hitRangeByDet.begin(); iHit != hitRangeByDet.end(); iHit++) {
0580         const TrackingRecHit* currHit = *(iHit->first);
0581         currDet = currHit->geographicalId();
0582 
0583         if (find(processedDets.begin(), processedDets.end(), currDet.rawId()) != processedDets.end())
0584           continue;
0585 
0586         TSOS prSt =
0587             currPropagator->propagate(traj.lastMeasurement().updatedState(), tracker->idToDet(currDet)->surface());
0588         if ((!prSt.isValid()) ||
0589             (theEstimator->Chi2MeasurementEstimatorBase::estimate(prSt, tracker->idToDet(currDet)->surface()) == false))
0590           //          if ( ( !prSt.isValid() ) ||  (theEstimator->estimate(prSt,tracker->idToDet(currDet)->surface() ) == false) )
0591           continue;
0592 
0593         trackHitCandidates.push_back(make_pair(iHit, prSt));
0594       }
0595 
0596       if (trackHitCandidates.empty())
0597         break;
0598 
0599       if (debug_info)
0600         cout << Hits.size() << " (int) trackHitCandidates.begin() " << trackHitCandidates.size() << endl;
0601       if (debug_info)
0602         cout << "Before sorting ... " << endl;
0603 
0604       if (debug_info)
0605         for (iHitRange = trackHitCandidates.begin(); iHitRange != trackHitCandidates.end(); iHitRange++) {
0606           if (debug_info)
0607             cout << (tracker->idToDet((*(iHitRange->first->first))->geographicalId()))->position();
0608         }
0609       if (debug_info)
0610         cout << endl;
0611 
0612       stable_sort(trackHitCandidates.begin(),
0613                   trackHitCandidates.end(),
0614                   CompareDetByTraj(traj.lastMeasurement().updatedState()));
0615 
0616       if (debug_info)
0617         cout << "After sorting ... " << endl;
0618       if (debug_info) {
0619         for (iHitRange = trackHitCandidates.begin(); iHitRange != trackHitCandidates.end(); iHitRange++) {
0620           if (debug_info)
0621             cout << (tracker->idToDet((*(iHitRange->first->first))->geographicalId()))->position();
0622         }
0623         cout << endl;
0624       }
0625 
0626       for (iHitRange = trackHitCandidates.begin(); iHitRange != trackHitCandidates.end(); iHitRange++)  //loop over dets
0627       {
0628         //now find the best hit of the detector
0629         if (debug_info)
0630           cout << "loop2 " << trackHitCandidates.size() << " " << trackHitCandidates.end() - iHitRange << endl;
0631         const TrackingRecHit* currHit = *(iHitRange->first->first);
0632 
0633         TransientTrackingRecHit::RecHitPointer bestHit = RHBuilder->build(currHit);
0634         TSOS currPrSt = (*iHitRange).second;
0635 
0636         if (debug_info)
0637           cout << "curr position" << bestHit->globalPosition();
0638         for (TrackingRecHitIterator iHit = (*iHitRange).first->first + 1; iHit != iHitRange->first->second; iHit++) {
0639           TransientTrackingRecHit::RecHitPointer tmpHit = RHBuilder->build(*iHit);
0640           if (debug_info)
0641             cout << "curr position" << tmpHit->globalPosition();
0642         }
0643       }
0644       if (debug_info)
0645         cout << "Cross check end ..." << endl;
0646 
0647       //just a simple test if the same hit can be added twice ...
0648       //      for( iHitRange = trackHitCandidates.begin(); iHitRange != trackHitCandidates.end(); iHitRange++ )      //loop over all hits
0649 
0650       //      break;
0651 
0652       for (iHitRange = trackHitCandidates.begin(); iHitRange != trackHitCandidates.end();
0653            iHitRange++)  //loop over detsall hits
0654       {
0655         //now find the best hit of the detector
0656         if (debug_info)
0657           cout << "loop2 " << trackHitCandidates.size() << " " << trackHitCandidates.end() - iHitRange << endl;
0658 
0659         const TrackingRecHit* currHit = *(iHitRange->first->first);
0660 
0661         processedDets.push_back(currHit->geographicalId().rawId());
0662 
0663         TransientTrackingRecHit::RecHitPointer bestHit = RHBuilder->build(currHit);
0664 
0665         if (debug_info)
0666           cout << "curr position A" << bestHit->globalPosition() << endl;
0667         TSOS currPrSt = (*iHitRange).second;
0668         double chi2min = theEstimator->estimate(currPrSt, *bestHit).second;
0669 
0670         if (debug_info)
0671           cout << "Size " << iHitRange->first->second - (*iHitRange).first->first << endl;
0672         for (TrackingRecHitIterator iHit = (*iHitRange).first->first + 1; iHit != iHitRange->first->second; iHit++) {
0673           if (debug_info)
0674             cout << "loop3 "
0675                  << " " << Hits.end() - iHit << endl;
0676 
0677           TransientTrackingRecHit::RecHitPointer tmpHit = RHBuilder->build(*iHit);
0678           if (debug_info)
0679             cout << "curr position B" << tmpHit->globalPosition() << endl;
0680           double currChi2 = theEstimator->estimate(currPrSt, *tmpHit).second;
0681           if (currChi2 < chi2min) {
0682             if (debug_info)
0683               cout << "Is best hit" << endl;
0684             chi2min = currChi2;
0685             bestHit = tmpHit;
0686           }
0687         }
0688         //now we have checked the det and can remove the entry from the vector...
0689 
0690         //if (debug_info) cout << "before erase ..." << endl;
0691         //this is to slow ...
0692         // hitRangeByDet.erase( (*iHitRange).first,(*iHitRange).first+1 );
0693         //if (debug_info) cout << "after erase ..." << endl;
0694 
0695         if (debug_info)
0696           cout << chi2min << endl;
0697         //if the track can be added to the trajectory
0698         if (chi2min < chi2cut) {
0699           if (debug_info)
0700             cout << "chi2 fine : " << chi2min << endl;
0701 
0702           //  if (debug_info) cout << "previaous state " << traj.lastMeasurement().updatedState() <<endl;
0703           TSOS UpdatedState = theUpdator->update(currPrSt, *bestHit);
0704           if (UpdatedState.isValid()) {
0705             hits.push_back(bestHit);
0706             traj.push(TM(currPrSt, UpdatedState, bestHit, chi2min));
0707             if (debug_info)
0708               edm::LogInfo("CosmicTrackFinder") << "STATE UPDATED WITH HIT AT POSITION "
0709                                                 //   <<tmphitbestdet->globalPosition()
0710                                                 << UpdatedState << " " << traj.chiSquared();
0711             if (debug_info)
0712               cout << "Added Hit" << bestHit->globalPosition() << endl;
0713             if (debug_info)
0714               cout << "State is valid ..." << UpdatedState << endl;
0715             //cout << "updated state " << traj.lastMeasurement().updatedState() <<endl;
0716 
0717             //        return; //break;
0718             //
0719             //        TSOS prSt= currPropagator->propagate(traj.lastMeasurement().updatedState(),
0720             //                        tracker->idToDet( bestHit->geographicalId() )->surface());
0721             //
0722             //        if ( prSt.isValid())
0723             //      cout << "the same hit can be added twice ..." << endl;
0724             //
0725 
0726             break;
0727           } else {
0728             if (debug_info)
0729               edm::LogWarning("CosmicTrackFinder")
0730                   << " State can not be updated with hit at position " << bestHit->globalPosition();
0731             //      cout << "State can not be updated with hit at " << bestHit->globalPosition() << endl;
0732           }
0733           continue;
0734         } else {
0735           //      cout << "chi2 to big : " << chi2min << endl;
0736         }
0737         if (debug_info)
0738           cout << " continue 1 " << endl;
0739       }
0740       //now we remove all already processed dets from the list ...
0741       // hitRangeByDet.erase( (*trackHitCandidates.begin()).first,(*iHitRange).first+1 );
0742 
0743       if (debug_info)
0744         cout << " continue 2 " << endl;
0745     }
0746     //if this was the last exit
0747     while (iHitRange != trackHitCandidates.end());
0748   }
0749 }
0750 
0751 bool CRackTrajectoryBuilder::qualityFilter(const Trajectory& traj) {
0752   int ngoodhits = 0;
0753   if (geometry == "MTCC") {
0754     auto hits = traj.recHits();
0755     for (auto hit = hits.begin(); hit != hits.end(); hit++) {
0756       unsigned int iid = (*hit)->hit()->geographicalId().rawId();
0757       //CHECK FOR 3 hits r-phi
0758       if (((iid >> 0) & 0x3) != 1)
0759         ngoodhits++;
0760     }
0761   } else
0762     ngoodhits = traj.foundHits();
0763 
0764   if (ngoodhits >= theMinHits) {
0765     return true;
0766   } else {
0767     return false;
0768   }
0769 }
0770 
0771 //----------------------------------------------------------------------------------------------------------
0772 // little helper function that returns false if both hits conatin the same information
0773 
0774 bool CRackTrajectoryBuilder::isDifferentStripReHit2D(const SiStripRecHit2D& hitA, const SiStripRecHit2D& hitB) {
0775   if (hitA.geographicalId() != hitB.geographicalId())
0776     return true;
0777   if (hitA.localPosition().x() != hitB.localPosition().x())
0778     return true;
0779   if (hitA.localPosition().y() != hitB.localPosition().y())
0780     return true;
0781   if (hitA.localPosition().z() != hitB.localPosition().z())
0782     return true;
0783 
0784   //  if (debug_info) cout << hitA.localPosition() << endl;
0785   //  if (debug_info) cout << hitB << endl;
0786 
0787   return false;
0788 }
0789 
0790 //----backfitting
0791 std::pair<TrajectoryStateOnSurface, const GeomDet*> CRackTrajectoryBuilder::innerState(const Trajectory& traj) const {
0792   int lastFitted = 999;
0793   int nhits = traj.foundHits();
0794   if (nhits < lastFitted + 1)
0795     lastFitted = nhits - 1;
0796 
0797   std::vector<TrajectoryMeasurement> measvec = traj.measurements();
0798   TransientTrackingRecHit::ConstRecHitContainer firstHits;
0799 
0800   bool foundLast = false;
0801   int actualLast = -99;
0802   for (int i = lastFitted; i >= 0; i--) {
0803     if (measvec[i].recHit()->isValid()) {
0804       if (!foundLast) {
0805         actualLast = i;
0806         foundLast = true;
0807       }
0808       firstHits.push_back(measvec[i].recHit());
0809     }
0810   }
0811   TSOS unscaledState = measvec[actualLast].updatedState();
0812   AlgebraicSymMatrix55 C = ROOT::Math::SMatrixIdentity();
0813   // C *= 100.;
0814 
0815   TSOS startingState(unscaledState.localParameters(),
0816                      LocalTrajectoryError(C),
0817                      unscaledState.surface(),
0818                      thePropagator->magneticField());
0819 
0820   // cout << endl << "FitTester starts with state " << startingState << endl;
0821 
0822   KFTrajectoryFitter backFitter(*thePropagator, KFUpdator(), Chi2MeasurementEstimator(100., 3));
0823 
0824   PropagationDirection backFitDirection = traj.direction() == alongMomentum ? oppositeToMomentum : alongMomentum;
0825 
0826   // only direction matters in this contest
0827   TrajectorySeed fakeSeed = TrajectorySeed(PTrajectoryStateOnDet(), edm::OwnVector<TrackingRecHit>(), backFitDirection);
0828 
0829   vector<Trajectory> fitres = backFitter.fit(fakeSeed, firstHits, startingState);
0830 
0831   if (fitres.size() != 1) {
0832     // cout << "FitTester: first hits fit failed!" << endl;
0833     return std::pair<TrajectoryStateOnSurface, const GeomDet*>();
0834   }
0835 
0836   TrajectoryMeasurement firstMeas = fitres[0].lastMeasurement();
0837   const TSOS& firstState = firstMeas.updatedState();
0838 
0839   //  cout << "FitTester: Fitted first state " << firstState << endl;
0840   //cout << "FitTester: chi2 = " << fitres[0].chiSquared() << endl;
0841 
0842   TSOS initialState(
0843       firstState.localParameters(), LocalTrajectoryError(C), firstState.surface(), thePropagator->magneticField());
0844 
0845   return std::pair<TrajectoryStateOnSurface, const GeomDet*>(initialState, firstMeas.recHit()->det());
0846 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!