NuclearSeedGenerator/src/NuclearTester.cc

0001 #include "RecoTracker/NuclearSeedGenerator/interface/NuclearTester.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 //----------------------------------------------------------------------
0005 NuclearTester::NuclearTester(unsigned int max_hits, const MeasurementEstimator* est, const TrackerGeometry* track_geom)
0006     : maxHits(max_hits), theEstimator(est), trackerGeom(track_geom) {
0007   NuclearIndex = 0;
0008 }
0009
0010 //----------------------------------------------------------------------
0011 bool NuclearTester::isNuclearInteraction() {
0012   // TODO : if energy of primary track is below a threshold don't use checkWithMultiplicity but only checkwith compatible_hits.front
0013
0014   // 1. if momentum of the primary track is below 5 GeV and if the number of compatible hits >0
0015   //    assume that a nuclear interaction occured at the end of the track
0016   if (allTM.front().first.updatedState().globalMomentum().mag() < 5.0 && compatible_hits.front() > 0) {
0017     NuclearIndex = 1;
0018     return true;
0019   }
0020
0021   // 2. else to use multiplicity we require at least 3 TM vectors to check if nuclear interactions occurs
0022   if (nHitsChecked() < 3)
0023     return false;
0024
0025   // 2. check with multiplicity :
0026   if (checkWithMultiplicity() == true)
0027     return true;
0028   else {
0029     // 3. last case : uncompleted track with at least 1 compatible hits in the last layer
0030     if (nHitsChecked() >= maxHits && compatible_hits.front() > 0) {
0031       NuclearIndex = 1;
0032       return true;
0033     }
0034   }
0035
0036   return false;
0037 }
0038 //----------------------------------------------------------------------
0039 bool NuclearTester::checkWithMultiplicity() {
0040   //RQ: assume that the input vector of compatible hits has been filled from Outside to Inside the tracker !
0041
0042   // find the first min nb of compatible TM :
0043   std::vector<int>::iterator min_it = min_element(compatible_hits.begin(), compatible_hits.end());
0044
0045   // if the outermost hit has no compatible TM, min_it has to be recalculated :
0046   if (min_it == compatible_hits.begin() && *min_it != 0)
0047     return false;
0048   if (min_it == compatible_hits.begin() && *min_it == 0)
0049     min_it = min_element(compatible_hits.begin() + 1, compatible_hits.end());
0050
0051   // this first min cannot be the innermost TM :
0052   if (min_it == compatible_hits.end() - 1)
0053     return false;
0054
0055   // if the previous nb of compatible TM is > min+2 and if the next compatible TM is min+-1 -> NUCLEAR
0056   // example : Nhits = 5, 8, 2, 2, ...
0057   if ((*(min_it - 1) - *min_it) > 2 && (*(min_it + 1) - *min_it) < 2) {
0058     NuclearIndex = min_it - compatible_hits.begin();
0059     return true;
0060   }
0061
0062   // case of : Nhits = 5, 8, 3, 2, 2, ...
0063   if (min_it - 1 != compatible_hits.begin())  //because min_it must be at least at the third position
0064   {
0065     if (min_it - 1 != compatible_hits.begin() && (*(min_it - 1) - *min_it) < 2 && (*(min_it - 2) - *(min_it - 1)) > 2) {
0066       NuclearIndex = min_it - 1 - compatible_hits.begin();
0067       return true;
0068     }
0069   }
0070
0071   return false;
0072 }
0073
0074 //----------------------------------------------------------------------
0075 double NuclearTester::meanHitDistance(const std::vector<TrajectoryMeasurement>& vecTM) const {
0076   std::vector<GlobalPoint> vgp = this->HitPositions(vecTM);
0077   double mean_dist = 0;
0078   int ncomb = 0;
0079   if (vgp.size() < 2)
0080     return 0;
0081   for (std::vector<GlobalPoint>::iterator itp = vgp.begin(); itp != vgp.end() - 1; itp++) {
0082     for (std::vector<GlobalPoint>::iterator itq = itp + 1; itq != vgp.end(); itq++) {
0083       double dist = ((*itp) - (*itq)).mag();
0084       // to calculate mean distance between particles and not hits (to not take into account twice stereo hits)
0085       if (dist > 1E-12) {
0086         mean_dist += dist;
0087         ncomb++;
0088       }
0089     }
0090   }
0091   return mean_dist / ncomb;
0092 }
0093 //----------------------------------------------------------------------
0094 std::vector<GlobalPoint> NuclearTester::HitPositions(const std::vector<TrajectoryMeasurement>& vecTM) const {
0095   std::vector<GlobalPoint> gp;
0096
0097   std::vector<TM>::const_iterator last = this->lastValidTM(vecTM);
0098
0099   for (std::vector<TrajectoryMeasurement>::const_iterator itm = vecTM.begin(); itm != last; itm++) {
0100     ConstRecHitPointer trh = itm->recHit();
0101     if (trh->isValid())
0102       gp.push_back(trackerGeom->idToDet(trh->geographicalId())->surface().toGlobal(trh->localPosition()));
0103   }
0104   return gp;
0105 }
0106 //----------------------------------------------------------------------
0107 double NuclearTester::fwdEstimate(const std::vector<TrajectoryMeasurement>& vecTM) const {
0108   if (vecTM.empty())
0109     return 0;
0110
0111   auto hit = vecTM.front().recHit().get();
0112   if (hit->isValid())
0113     return theEstimator->estimate(vecTM.front().forwardPredictedState(), *hit).second;
0114   else
0115     return -1;
0116   /*
0117        double meanEst=0;
0118        int    goodTM=0;
0119        std::vector<TM>::const_iterator last;
0120        //std::vector<TM>::const_iterator last = this->lastValidTM(vecTM);
0121        if(vecTM.size() > 2) last = vecTM.begin()+2;
0122        else last = vecTM.end();
0123
0124        for(std::vector<TrajectoryMeasurement>::const_iterator itm = vecTM.begin(); itm!=last; itm++) {
0125              meanEst += itm->estimate();
0126              goodTM++;
0127        }
0128        return meanEst/goodTM;
0129 */
0130 }
0131 //----------------------------------------------------------------------
0132 std::vector<TrajectoryMeasurement>::const_iterator NuclearTester::lastValidTM(const std::vector<TM>& vecTM) const {
0133   if (vecTM.empty())
0134     return vecTM.end();
0135   if (vecTM.front().recHit()->isValid())
0136     return std::find_if(vecTM.begin(), vecTM.end(), [](auto const& meas) { return !meas.recHit()->isValid(); });
0137   else
0138     return vecTM.end();
0139 }
0140 //----------------------------------------------------------------------