SiTrackerMRHTools/src/GroupedDAFHitCollector.cc

0001 #include "RecoTracker/SiTrackerMRHTools/interface/GroupedDAFHitCollector.h"
0002 #include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h"
0003 #include "RecoTracker/SiTrackerMRHTools/interface/MeasurementByLayerGrouper.h"
0004 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0005 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0006 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0007 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0008 #include "TrackingTools/MeasurementDet/interface/TrajectoryMeasurementGroup.h"
0009 #include "TrackingTools/TransientTrackingRecHit/interface/InvalidTransientRecHit.h"
0010 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0011 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0012
0013 #include <vector>
0014 #include <map>
0015
0016 using namespace std;
0017
0018 vector<TrajectoryMeasurement> GroupedDAFHitCollector::recHits(const Trajectory& traj,
0019                                                               const MeasurementTrackerEvent* theMTE) const {
0020   LayerMeasurements theLM(theMTE->measurementTracker(), *theMTE);
0021
0022   //WARNING: At the moment the trajectories has the measurements with reversed sorting after the track smoothing
0023   const vector<TrajectoryMeasurement>& meas = traj.measurements();
0024   const Propagator* forwardPropagator = getPropagator();
0025   const Propagator* backwardPropagator = getReversePropagator();
0026   if (traj.direction() == alongMomentum) {
0027     forwardPropagator = getReversePropagator();
0028     backwardPropagator = getPropagator();
0029   }
0030   if (meas.empty())  //return TransientTrackingRecHit::ConstRecHitContainer();
0031     return vector<TrajectoryMeasurement>();
0032
0033   //groups the TrajectoryMeasurements on a layer by layer
0034   vector<pair<const DetLayer*, vector<TrajectoryMeasurement> > > mol;
0035   mol = MeasurementByLayerGrouper(getMeasurementTracker()->geometricSearchTracker())(meas);
0036
0037   vector<TrajectoryMeasurement> result;
0038
0039   //add a protection if all the measurement are on the same layer
0040   if (mol.size() < 2)
0041     return vector<TrajectoryMeasurement>();
0042
0043   //it assumes that the measurements are sorted in the smoothing direction
0044   //TrajectoryStateOnSurface current = (*(mol.begin()+1)).second.front().updatedState();
0045   TrajectoryStateOnSurface current = (*(mol.rbegin() + 1)).second.back().updatedState();
0046   //if (current.isValid()) current.rescaleError(10);
0047
0048   //protection for layers with invalid meas with no id associated
0049   //to be fixed
0050   //for the moment no hit are lookerd for in these layers
0051   //remind that:
0052   //groupedMeasurements will return at least a measurement with an invalid hit with no detid
0053   LogDebug("MultiRecHitCollector") << "Layer " << mol.back().first << " has " << mol.back().second.size()
0054                                    << " measurements";
0055   LogTrace("MultiRecHitCollector") << "Original measurements are:";
0056   for (unsigned int iLay = 0; iLay < mol.size(); iLay++) {
0057     LogTrace("MultiRecHitCollector") << "  Layer " << mol.at(iLay).first << " has " << mol.at(iLay).second.size()
0058                                      << " measurements:";
0059     vector<TrajectoryMeasurement>::const_iterator ibeg = (mol.at(iLay)).second.begin();
0060     vector<TrajectoryMeasurement>::const_iterator iend = (mol.at(iLay)).second.end();
0061     for (vector<TrajectoryMeasurement>::const_iterator imeas = ibeg; imeas != iend; ++imeas) {
0062       if (imeas->recHit()->isValid()) {
0063         LogTrace("MultiRecHitCollector") << "   Valid Hit with DetId " << imeas->recHit()->geographicalId().rawId()
0064                                          << " local position " << imeas->recHit()->hit()->localPosition()
0065                                          << " global position " << imeas->recHit()->hit()->globalPosition();
0066       } else {
0067         LogTrace("MultiRecHitCollector") << "   Invalid Hit with DetId " << imeas->recHit()->geographicalId().rawId();
0068       }
0069     }
0070   }
0071
0072   //ERICA: I have to understand how are set the TM now. REPLACE THIS PART!!
0073   vector<TrajectoryMeasurementGroup> groupedMeas;
0074   if (mol.back().first)
0075     groupedMeas = theLM.groupedMeasurements(*(mol.back().first), current, *backwardPropagator, *(getEstimator()));
0076
0077   //Since we have passed the backwardPropagator, we have to sort the detGroups in the opposite way
0078   //(according the forward propagator, not the backward one)
0079   vector<TrajectoryMeasurementGroup> sortedgroupedMeas;
0080   for (vector<TrajectoryMeasurementGroup>::reverse_iterator iter = groupedMeas.rbegin(); iter != groupedMeas.rend();
0081        iter++) {
0082     sortedgroupedMeas.push_back(*iter);
0083   }
0084
0085   //for the first layer
0086   buildMultiRecHits(sortedgroupedMeas, result, theMTE);
0087
0088   //for other layers
0089   current = mol.back().second.front().updatedState();
0090   //if (current.isValid()) current.rescaleError(10);
0091
0092   for (vector<pair<const DetLayer*, vector<TrajectoryMeasurement> > >::reverse_iterator imol = mol.rbegin() + 1;
0093        imol != mol.rend();
0094        imol++) {
0095     const DetLayer* lay = (*imol).first;
0096     LogDebug("MultiRecHitCollector") << "Layer " << lay << " has " << (*imol).second.size() << " measurements";
0097     //debug
0098     vector<TrajectoryMeasurementGroup> currentLayerMeas;
0099     if (lay) {
0100       currentLayerMeas = theLM.groupedMeasurements(*lay, current, *forwardPropagator, *(getEstimator()));
0101     }
0102
0103     buildMultiRecHits(currentLayerMeas, result, theMTE);
0104     current = (*imol).second.front().updatedState();
0105     //if (current.isValid()) current.rescaleError(10);
0106   }
0107
0108   LogTrace("MultiRecHitCollector") << " Ending GroupedDAFHitCollector::recHits >> Original Measurement size "
0109                                    << meas.size()
0110                                    << "\n                                      >> GroupedDAFHitCollector returned "
0111                                    << result.size() << " measurements";
0112   //results are sorted in the fitting direction
0113
0114   //    adding a protection against too few hits and invalid hits (due to failed propagation on the same surface of the original hits)
0115   if (result.size() > 2) {
0116     int hitcounter = 0;
0117     //check if the vector result has more than 3 valid hits
0118     for (vector<TrajectoryMeasurement>::const_iterator iimeas = result.begin(); iimeas != result.end(); ++iimeas) {
0119       if (iimeas->recHit()->isValid())
0120         hitcounter++;
0121     }
0122
0123     if (hitcounter > 2) {
0124       return result;
0125     }
0126
0127     else
0128       return vector<TrajectoryMeasurement>();
0129   }
0130
0131   else {
0132     return vector<TrajectoryMeasurement>();
0133   }
0134 }
0135
0136 void GroupedDAFHitCollector::buildMultiRecHits(const vector<TrajectoryMeasurementGroup>& measgroup,
0137                                                vector<TrajectoryMeasurement>& result,
0138                                                const MeasurementTrackerEvent*& theMTE) const {
0139   unsigned int initial_size = result.size();
0140
0141   //TransientTrackingRecHit::ConstRecHitContainer rhits;
0142   if (measgroup.empty()) {
0143     LogTrace("MultiRecHitCollector") << "No TrajectoryMeasurementGroups found for this layer\n";
0144     //should we do something?
0145     //result.push_back(InvalidTransientRecHit::build(0,TrackingRecHit::missing));
0146     return;
0147   }
0148
0149   //we build a MultiRecHit out of each group
0150   //groups are sorted along momentum or opposite to momentum,
0151   //measurements in groups are sorted with increating chi2
0152   LogTrace("MultiRecHitCollector") << "Found " << measgroup.size() << " groups for this layer";
0153
0154   //trajectory state to store the last valid TrajectoryState (if present) to be used
0155   //to add an invalid Measurement in case no valid state or no valid hits are found in any group
0156   for (vector<TrajectoryMeasurementGroup>::const_iterator igroup = measgroup.begin(); igroup != measgroup.end();
0157        igroup++) {
0158     //the TrajectoryState is the first one
0159     TrajectoryStateOnSurface state = igroup->measurements().front().predictedState();
0160     if (!state.isValid()) {
0161       LogTrace("MultiRecHitCollector") << "Something wrong! no valid TSOS found in current group ";
0162       continue;
0163     }
0164
0165     LogTrace("MultiRecHitCollector") << "This group has " << igroup->measurements().size() << " measurements";
0166     LogTrace("MultiRecHitCollector") << "This group has the following " << igroup->detGroup().size()
0167                                      << " detector ids: " << endl;
0168     for (DetGroup::const_iterator idet = igroup->detGroup().begin(); idet != igroup->detGroup().end(); ++idet) {
0169       LogTrace("MultiRecHitCollector") << idet->det()->geographicalId().rawId();
0170     }
0171
0172     vector<const TrackingRecHit*> hits;
0173     for (vector<TrajectoryMeasurement>::const_iterator imeas = igroup->measurements().begin();
0174          imeas != igroup->measurements().end();
0175          imeas++) {
0176       //should be fixed!!
0177       //DetId id = imeas->recHit()->geographicalId();
0178       //MeasurementDetWithData measDet = theMTE->idToDet(id);
0179
0180       //collect the non missing hits to build the MultiRecHits
0181       //we use the recHits method; anyway only simple hits, not MultiHits should be present
0182       if (imeas->recHit()->getType() != TrackingRecHit::missing) {
0183         LogTrace("MultiRecHitCollector") << "This hit is valid ";
0184         hits.push_back(imeas->recHit()->hit());
0185       } else {
0186         LogTrace("MultiRecHitCollector") << "     This hit is not valid and will not enter in the MRH. ";
0187       }
0188     }
0189
0190     if (hits.empty()) {
0191       LogTrace("MultiRecHitCollector") << "No valid hits found in current group ";
0192       continue;
0193     }
0194
0195     LogTrace("MultiRecHitCollector") << "The best TSOS in this group is " << state << " it lays on surface located at "
0196                                      << state.surface().position();
0197
0198     LogTrace("MultiRecHitCollector") << "For the MRH on this group the following hits will be used";
0199     for (vector<const TrackingRecHit*>::iterator iter = hits.begin(); iter != hits.end(); iter++) {
0200       string validity = "valid";
0201       if ((*iter)->getType() == TrackingRecHit::missing)
0202         validity = "missing !should not happen!";
0203       else if ((*iter)->getType() == TrackingRecHit::inactive)
0204         validity = "inactive";
0205       else if ((*iter)->getType() == TrackingRecHit::bad)
0206         validity = "bad";
0207       LogTrace("MultiRecHitCollector")
0208           << "DetId " << (*iter)->geographicalId().rawId() << " validity: " << validity << " surface position "
0209           << getMeasurementTracker()->geomTracker()->idToDet((*iter)->geographicalId())->position()
0210           << " hit local position " << (*iter)->localPosition();
0211     }
0212     //should be fixed!!
0213     //result.push_back(TrajectoryMeasurement(state,theUpdator->buildMultiRecHit(hits, state, *MeasurementDetWithData())));
0214   }
0215   //can this happen? it means that the measgroup was not empty but no valid measurement was found inside
0216   //in this case we add an invalid measuremnt for this layer
0217   if (result.size() == initial_size) {
0218     LogTrace("MultiRecHitCollector") << "no valid measuremnt or no valid TSOS in none of the groups";
0219     //measgroup has been already checked for size != 0
0220     if (!measgroup.back().measurements().empty()) {
0221       result.push_back(measgroup.back().measurements().back());
0222     }
0223   }
0224 }