Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:15

0001 #include "RecoEgamma/EgammaPhotonAlgos/interface/InOutConversionSeedFinder.h"
0002 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionBarrelEstimator.h"
0003 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionForwardEstimator.h"
0004 
0005 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionFastHelix.h"
0006 
0007 // Field
0008 //
0009 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0010 // Geometry
0011 //
0012 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0013 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0015 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0016 
0017 //
0018 //
0019 //
0020 
0021 InOutConversionSeedFinder::InOutConversionSeedFinder(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC)
0022     : ConversionSeedFinder(conf, iC), conf_(conf) {
0023   // std::cout << " InOutConversionSeedFinder CTOR " << "\n";
0024 
0025   maxNumberOfInOutSeedsPerInputTrack_ = conf_.getParameter<int>("maxNumOfSeedsInOut");
0026   //the2ndHitdphi_ = 0.008;
0027   the2ndHitdphi_ = 0.01;
0028   the2ndHitdzConst_ = 5.;
0029   the2ndHitdznSigma_ = 2.;
0030 }
0031 
0032 InOutConversionSeedFinder::~InOutConversionSeedFinder() {
0033   //std::cout << " InOutConversionSeedFinder DTOR " << "\n";
0034 }
0035 
0036 void InOutConversionSeedFinder::makeSeeds(const edm::Handle<edm::View<reco::CaloCluster> >& allBC) {
0037   //std::cout << "  InOutConversionSeedFinder::makeSeeds() " << "\n";
0038   theSeeds_.clear();
0039   //std::cout << " Check Calo cluster collection size " << allBC->size() << "\n";
0040   bcCollection_ = allBC;
0041 
0042   findLayers();
0043 
0044   fillClusterSeeds();
0045   //std::cout << "Built vector of seeds of size  " << theSeeds_.size() <<  "\n" ;
0046 
0047   theOutInTracks_.clear();
0048   inputTracks_.clear();
0049   theFirstMeasurements_.clear();
0050 }
0051 
0052 void InOutConversionSeedFinder::fillClusterSeeds() {
0053   std::vector<Trajectory>::const_iterator outInTrackItr;
0054 
0055   //std::cout << "  InOutConversionSeedFinder::fillClusterSeeds outInTracks_.size " << theOutInTracks_.size() << "\n";
0056   //Start looking for seeds for both of the 2 best tracks from the inward tracking
0057 
0058   ///// This bit is for debugging; it will go away
0059   /*
0060   for(outInTrackItr = theOutInTracks_.begin(); outInTrackItr != theOutInTracks_.end();  ++outInTrackItr) {
0061 
0062 
0063    //std::cout << " InOutConversionSeedFinder::fillClusterSeeds out in input track hits " << (*outInTrackItr).foundHits() << "\n";
0064     DetId tmpId = DetId( (*outInTrackItr).seed().startingState().detId());
0065     const GeomDet* tmpDet  = this->getMeasurementTracker()->geomTracker()->idToDet( tmpId );
0066     GlobalVector gv = tmpDet->surface().toGlobal( (*outInTrackItr).seed().startingState().parameters().momentum() );
0067 
0068 
0069    //std::cout << " InOutConversionSeedFinder::fillClusterSeed was built from seed position " <<gv   <<  " charge " << (*outInTrackItr).seed().startingState().parameters().charge() << "\n";
0070 
0071     Trajectory::DataContainer m=  outInTrackItr->measurements();
0072     int nHit=0;
0073     for (Trajectory::DataContainer::iterator itm = m.begin(); itm != m.end(); ++itm) {
0074       if ( itm->recHit()->isValid()  ) {
0075         nHit++;
0076     //std::cout << nHit << ")  Valid RecHit global position " << itm->recHit()->globalPosition() << " R " <<  itm->recHit()->globalPosition().perp() << " phi " << itm->recHit()->globalPosition().phi() << " eta " << itm->recHit()->globalPosition().eta() << "\n";
0077       } 
0078 
0079     }
0080 
0081   }
0082 
0083   */
0084 
0085   //Start looking for seeds for both of the 2 best tracks from the inward tracking
0086   for (outInTrackItr = theOutInTracks_.begin(); outInTrackItr != theOutInTracks_.end(); ++outInTrackItr) {
0087     //std::cout << " InOutConversionSeedFinder::fillClusterSeeds out in input track hits " << (*outInTrackItr).foundHits() << "\n";
0088     nSeedsPerInputTrack_ = 0;
0089 
0090     //Find the first valid hit of the track
0091     // Measurements are ordered according to the direction in which the trajectories were built
0092     std::vector<TrajectoryMeasurement> measurements = (*outInTrackItr).measurements();
0093 
0094     std::vector<const DetLayer*> allLayers = layerList();
0095 
0096     //std::cout << "  InOutConversionSeedFinder::fill clusterSeed allLayers.size " <<  allLayers.size() << "\n";
0097     for (unsigned int i = 0; i < allLayers.size(); ++i) {
0098       //std::cout <<  " allLayers " << allLayers[i] << "\n";
0099       printLayer(i);
0100     }
0101 
0102     std::vector<const DetLayer*> myLayers;
0103     myLayers.clear();
0104     std::vector<TrajectoryMeasurement>::reverse_iterator measurementItr;
0105     std::vector<TrajectoryMeasurement*> myItr;
0106     // TrajectoryMeasurement* myPointer=0;
0107     myPointer = nullptr;
0108     //std::cout << "  InOutConversionSeedFinder::fillClusterSeeds measurements.size " << measurements.size() <<"\n";
0109 
0110     for (measurementItr = measurements.rbegin(); measurementItr != measurements.rend(); ++measurementItr) {
0111       if ((*measurementItr).recHit()->isValid()) {
0112         //std::cout << "  InOutConversionSeedFinder::fillClusterSeeds measurement on  layer  " << measurementItr->layer() <<   " " <<&(*measurementItr) <<  " position " << measurementItr->recHit()->globalPosition() <<   " R " << sqrt( measurementItr->recHit()->globalPosition().x()*measurementItr->recHit()->globalPosition().x() + measurementItr->recHit()->globalPosition().y()*measurementItr->recHit()->globalPosition().y() ) << " Z " << measurementItr->recHit()->globalPosition().z() << " phi " <<  measurementItr->recHit()->globalPosition().phi() << "\n";
0113 
0114         myLayers.push_back(measurementItr->layer());
0115         myItr.push_back(&(*measurementItr));
0116       }
0117     }
0118 
0119     //std::cout << " InOutConversionSeedFinder::fillClusterSeed myLayers.size " <<  myLayers.size() << "\n";
0120     //    for( unsigned int i = 0; i < myLayers.size(); ++i) {
0121     ////std::cout <<  " myLayers " << myLayers[i] << " myItr " << myItr[i] << "\n";
0122     // }
0123 
0124     if (myItr.empty()) {
0125       //std::cout << "HORRENDOUS ERROR!  No meas on track!" << "\n";
0126     }
0127     unsigned int ilayer;
0128     for (ilayer = 0; ilayer < allLayers.size(); ++ilayer) {
0129       //std::cout <<  " allLayers in the search loop  " << allLayers[ilayer] <<  " " << myLayers[0] <<  "\n";
0130       if (allLayers[ilayer] == myLayers[0]) {
0131         myPointer = myItr[0];
0132 
0133         //std::cout <<  " allLayers in the search loop   allLayers[ilayer] == myLayers[0])  " << allLayers[ilayer] <<  " " << myLayers[0] <<  " myPointer " << myPointer << "\n";
0134 
0135         //std::cout  << "Layer " << ilayer << "  contains the first valid measurement " << "\n";
0136         printLayer(ilayer);
0137 
0138         if ((myLayers[0])->location() == GeomDetEnumerators::barrel) {
0139           //      const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(myLayers[0]);
0140           //std::cout << " InOutConversionSeedFinder::fillClusterSeeds  **** firstHit found in Barrel on layer " << ilayer  << " R= " << barrelLayer->specificSurface().radius() <<   "\n";
0141         } else {
0142           //const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(myLayers[0]);
0143           //std::cout << " InOutwardConversionSeedFinder::fillClusterSeeds  **** firstHit found in Forw on layer " << ilayer  << " Z= " << forwardLayer->specificSurface().position().z() <<  "\n";
0144         }
0145 
0146         break;
0147 
0148       } else if (allLayers[ilayer] == myLayers[1]) {
0149         myPointer = myItr[1];
0150 
0151         //std::cout <<  " allLayers in the search loop   allLayers[ilayer] == myLayers[1])  " << allLayers[ilayer] <<  " " << myLayers[1] <<  " myPointer " << myPointer << "\n";
0152 
0153         //std::cout << "Layer " << ilayer << "  contains the first innermost  valid measurement " << "\n";
0154         if ((myLayers[1])->location() == GeomDetEnumerators::barrel) {
0155           //      const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(myLayers[1]);
0156           //std::cout << " InOutConversionSeedFinder::fillClusterSeeds  **** 2ndHit found in Barrel on layer " << ilayer  << " R= " << barrelLayer->specificSurface().radius() <<   "\n";
0157         } else {
0158           //const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(myLayers[1]);
0159           //std::cout << " InOutwardConversionSeedFinder::fillClusterSeeds  ****  2ndHitfound on forw layer " << ilayer  << " Z= " << forwardLayer->specificSurface().position().z() <<  "\n";
0160         }
0161 
0162         break;
0163       }
0164     }
0165 
0166     if (ilayer == allLayers.size()) {
0167       //std::cout << "InOutConversionSeedFinder::fillClusterSeeds ERROR could not find layer on list" <<  "\n";
0168       return;
0169     }
0170 
0171     //PropagatorWithMaterial reversePropagator(oppositeToMomentum, 0.000511, &(*theMF_) );
0172     assert(myPointer);
0173     const FreeTrajectoryState* fts = myPointer->updatedState().freeTrajectoryState();
0174 
0175     //std::cout << " InOutConversionSeedFinder::fillClusterSeeds First FTS charge " << fts->charge() << " Position " << fts->position() << " momentum " << fts->momentum() << " R " << sqrt(fts->position().x()*fts->position().x() + fts->position().y()* fts->position().y() ) << " Z " << fts->position().z() << " phi " << fts->position().phi() << " fts parameters " << fts->parameters() << "\n";
0176 
0177     while (ilayer > 0) {
0178       //std::cout << " InOutConversionSeedFinder::fillClusterSeeds looking for 2nd seed from layer " << ilayer << "\n";
0179 
0180       //   if ( (allLayers[ilayer])->location() == GeomDetEnumerators::barrel ) {const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(allLayers[ilayer]);
0181       //std::cout <<  " InOutConversionSeedFinder::fillClusterSeeds  ****  Barrel on layer " << ilayer  << " R= " << barrelLayer->specificSurface().radius() <<  "\n";
0182       // } else {
0183       //const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(allLayers[ilayer]);
0184       //std::cout <<  " InOutConversionSeedFinder::fillClusterSeeds  ****  Forw on layer " << ilayer  << " Z= " << forwardLayer->specificSurface().position().z() << "\n";
0185       //      }
0186 
0187       const DetLayer* previousLayer = allLayers[ilayer];
0188       TrajectoryStateOnSurface stateAtPreviousLayer;
0189       //std::cout << " InOutConversionSeedFinder::fillClusterSeeds previousLayer->surface() position before  " <<allLayers[ilayer] << " " <<  previousLayer->surface().position() << " layer location " << previousLayer->location() << "\n";
0190       // Propagate to the previous layer
0191       // The present layer is actually included in the loop so that a partner can be searched for
0192       // Applying the propagator to the same layer does not do any harm. It simply does nothing
0193 
0194       //     const Propagator& newProp=thePropagatorOppositeToMomentum_;
0195       //std::cout << " InOutConversionSeedFinder::fillClusterSeeds reversepropagator direction " << thePropagatorOppositeToMomentum_->propagationDirection()  << "\n";
0196       if (ilayer - 1 > 0) {
0197         if (allLayers[ilayer] == myLayers[0]) {
0198           //std::cout << " innermost hit R " << myPointer->recHit()->globalPosition().perp() << " Z " << myPointer->recHit()->globalPosition().z() << " phi " <<myPointer->recHit()->globalPosition().phi() << "\n";
0199           //std::cout << " surface R " << theTrackerGeom_->idToDet(  myPointer->recHit() ->geographicalId())->surface().position().perp() <<  " Z " <<  theTrackerGeom_->idToDet(  myPointer->recHit() ->geographicalId())->surface().position().z() << " phi " << theTrackerGeom_->idToDet(  myPointer->recHit() ->geographicalId())->surface().position().phi() << "\n";
0200 
0201           stateAtPreviousLayer = thePropagatorOppositeToMomentum_->propagate(
0202               *fts, theTrackerGeom_->idToDet(myPointer->recHit()->geographicalId())->surface());
0203 
0204         } else {
0205           stateAtPreviousLayer = thePropagatorOppositeToMomentum_->propagate(*fts, previousLayer->surface());
0206           //std::cout << " InOutConversionSeedFinder::fillClusterSeeds previousLayer->surface() position after " << previousLayer->surface().position() << " layer location " << previousLayer->location() <<   "\n";
0207         }
0208 
0209       } else if (ilayer - 1 == 0) {
0210         ////std::cout << " innermost hit R " << myPointer->recHit()->globalPosition().perp() << " Z " << myPointer->recHit()->globalPosition().z() << " phi " <<myPointer->recHit()->globalPosition().phi() << "\n";
0211         ////std::cout << " surface R " << theTrackerGeom_->idToDet(  myPointer->recHit() ->geographicalId())->surface().position().perp() <<  " Z " <<  theTrackerGeom_->idToDet(  myPointer->recHit() ->geographicalId())->surface().position().z() << " phi " << theTrackerGeom_->idToDet(  myPointer->recHit() ->geographicalId())->surface().position().phi() << "\n";
0212 
0213         //stateAtPreviousLayer= thePropagatorOppositeToMomentum_->propagate(*fts,   theTrackerGeom_->idToDet(  myPointer->recHit() ->geographicalId())->surface()   );
0214         stateAtPreviousLayer = thePropagatorOppositeToMomentum_->propagate(*fts, previousLayer->surface());
0215       }
0216 
0217       if (!stateAtPreviousLayer.isValid()) {
0218         //std::cout << "InOutConversionSeedFinder::fillClusterSeeds ERROR:could not propagate back to layer "  << ilayer << "\n";
0219         ////std::cout << "  InOutConversionSeedFinder::fillClusterSeeds stateAtPreviousLayer " << stateAtPreviousLayer <<std:: endl;
0220       } else {
0221         //std::cout << "InOutConversionSeedFinder::fillClusterSeeds  stateAtPreviousLayer is valid.  Propagating back to layer "  << ilayer << "\n";
0222         //std::cout << "InOutConversionSeedFinder::fillClusterSeeds stateAtPreviousLayer R  " << stateAtPreviousLayer.globalPosition().perp()  << " Z " << stateAtPreviousLayer.globalPosition().z() << " phi " <<  stateAtPreviousLayer.globalPosition().phi() << "\n";
0223 
0224         startSeed(fts, stateAtPreviousLayer, -1, ilayer);
0225       }
0226 
0227       --ilayer;
0228     }
0229 
0230     if (ilayer == 0) {
0231       //     if ( (allLayers[ilayer])->location() == GeomDetEnumerators::barrel ) {const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(allLayers[ilayer]);
0232       // //std::cout <<  " InOutConversionSeedFinder::fillClusterSeeds  ****  Barrel on layer " << ilayer  << " R= " << barrelLayer->specificSurface().radius() <<  "\n";
0233       // } else {
0234       //const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(allLayers[ilayer]);
0235       //std::cout <<  " InOutConversionSeedFinder::fillClusterSeeds  ****  Forw on layer " << ilayer  << " Z= " << forwardLayer->specificSurface().position().z() << "\n";
0236       // }
0237       const DetLayer* previousLayer = allLayers[ilayer];
0238       TrajectoryStateOnSurface stateAtPreviousLayer;
0239       stateAtPreviousLayer = thePropagatorOppositeToMomentum_->propagate(*fts, previousLayer->surface());
0240 
0241       if (!stateAtPreviousLayer.isValid()) {
0242         //std::cout << "InOutConversionSeedFinder::fillClusterSeeds ERROR:could not propagate back to layer "  << ilayer << "\n";
0243         ////std::cout << "  InOutConversionSeedFinder::fillClusterSeeds stateAtPreviousLayer " << stateAtPreviousLayer <<std:: endl;
0244       } else {
0245         //std::cout << "InOutConversionSeedFinder::fillClusterSeeds  stateAtPreviousLayer is valid.  Propagating back to layer "  << ilayer << "\n";
0246         //std::cout << "InOutConversionSeedFinder::fillClusterSeeds stateAtPreviousLayer R  " << stateAtPreviousLayer.globalPosition().perp()  << " Z " << stateAtPreviousLayer.globalPosition().z() << " phi " <<  stateAtPreviousLayer.globalPosition().phi() << "\n";
0247 
0248         startSeed(fts, stateAtPreviousLayer, -1, ilayer);
0249       }
0250     }
0251 
0252   }  // End loop over Out In tracks
0253 }
0254 
0255 void InOutConversionSeedFinder::startSeed(const FreeTrajectoryState* fts,
0256                                           const TrajectoryStateOnSurface& stateAtPreviousLayer,
0257                                           int charge,
0258                                           int ilayer) {
0259   //std::cout << "InOutConversionSeedFinder::startSeed ilayer " << ilayer <<  "\n";
0260   // Get a list of basic clusters that are consistent with a track
0261   // starting at the assumed conversion point with opp. charge to the
0262   // inward track.  Loop over these basic clusters.
0263   track2Charge_ = charge * fts->charge();
0264   std::vector<const reco::CaloCluster*> bcVec;
0265   //std::cout << "InOutConversionSeedFinder::startSeed charge assumed for the in-out track  " << track2Charge_ <<  "\n";
0266 
0267   // Geom::Phi<float> theConvPhi( stateAtPreviousLayer.globalPosition().phi());
0268   //std::cout << "InOutConversionSeedFinder::startSeed  stateAtPreviousLayer phi " << stateAtPreviousLayer.globalPosition().phi() << " R " <<  stateAtPreviousLayer.globalPosition().perp() << " Z " << stateAtPreviousLayer.globalPosition().z() << "\n";
0269 
0270   bcVec = getSecondCaloClusters(stateAtPreviousLayer.globalPosition(), track2Charge_);
0271 
0272   std::vector<const reco::CaloCluster*>::iterator bcItr;
0273   //std::cout << "InOutConversionSeedFinder::startSeed bcVec.size " << bcVec.size() << "\n";
0274 
0275   // debug
0276   //  for(bcItr = bcVec.begin(); bcItr != bcVec.end(); ++bcItr) {
0277   //  //std::cout << "InOutConversionSeedFinder::startSeed list of  bc eta " << (*bcItr)->position().eta() << " phi " << (*bcItr)->position().phi() << " x " << (*bcItr)->position().x() << " y " << (*bcItr)->position().y() << " z " << (*bcItr)->position().z() << "\n";
0278   // }
0279 
0280   for (bcItr = bcVec.begin(); bcItr != bcVec.end(); ++bcItr) {
0281     theSecondBC_ = **bcItr;
0282     GlobalPoint bcPos((theSecondBC_.position()).x(), (theSecondBC_.position()).y(), (theSecondBC_.position()).z());
0283 
0284     //std::cout << "InOutConversionSeedFinder::startSeed for  bc position x " << bcPos.x() << " y " <<  bcPos.y() << " z  " <<  bcPos.z() << " eta " <<  bcPos.eta() << " phi " <<  bcPos.phi() << "\n";
0285     GlobalVector dir = stateAtPreviousLayer.globalDirection();
0286     GlobalPoint back1mm = stateAtPreviousLayer.globalPosition();
0287     //std::cout << "InOutConversionSeedFinder::startSeed   stateAtPreviousLayer.globalPosition() " << back1mm << "\n";
0288 
0289     back1mm -= dir.unit() * 0.1;
0290     //std::cout << " InOutConversionSeedFinder:::startSeed going to make the helix using back1mm " << back1mm <<"\n";
0291     ConversionFastHelix helix(bcPos, stateAtPreviousLayer.globalPosition(), back1mm, &(*theMF_));
0292     helix.stateAtVertex();
0293 
0294     //std::cout << " InOutConversionSeedFinder:::startSeed helix status " <<helix.isValid() << std::endl;
0295     if (!helix.isValid())
0296       continue;
0297     findSeeds(stateAtPreviousLayer, helix.stateAtVertex().transverseCurvature(), ilayer);
0298   }
0299 }
0300 
0301 std::vector<const reco::CaloCluster*> InOutConversionSeedFinder::getSecondCaloClusters(
0302     const GlobalPoint& conversionPosition, float charge) const {
0303   std::vector<const reco::CaloCluster*> result;
0304 
0305   Geom::Phi<float> convPhi(conversionPosition.phi());
0306 
0307   for (auto const& bc : *bcCollection_) {
0308     Geom::Phi<float> bcPhi(bc.position().phi());
0309 
0310     // Require phi of cluster to be consistent with the conversion position and the track charge
0311 
0312     if (std::abs(bcPhi - convPhi) < .5 &&
0313         ((charge < 0 && bcPhi - convPhi > -.5) || (charge > 0 && bcPhi - convPhi < .5)))
0314       result.push_back(&bc);
0315   }
0316 
0317   return result;
0318 }
0319 
0320 void InOutConversionSeedFinder::findSeeds(const TrajectoryStateOnSurface& startingState,
0321                                           float transverseCurvature,
0322                                           unsigned int startingLayer) {
0323   std::vector<const DetLayer*> allLayers = layerList();
0324   //std::cout << "InOutConversionSeedFinder::findSeeds starting forward propagation from  startingLayer " << startingLayer <<  "\n";
0325 
0326   // create error matrix
0327   AlgebraicSymMatrix55 m = AlgebraicMatrixID();
0328   m(0, 0) = 0.1;
0329   m(1, 1) = 0.0001;
0330   m(2, 2) = 0.0001;
0331   m(3, 3) = 0.0001;
0332   m(4, 4) = 0.001;
0333 
0334   // Make an FTS consistent with the start point, start direction and curvature
0335   FreeTrajectoryState fts(
0336       GlobalTrajectoryParameters(
0337           startingState.globalPosition(), startingState.globalDirection(), double(transverseCurvature), 0, &(*theMF_)),
0338       CurvilinearTrajectoryError(m));
0339   if (fts.momentum().mag2() == 0) {
0340     edm::LogWarning("FailedToInitiateSeeding")
0341         << " initial FTS has a zero momentum, probably because of the zero field.  ";
0342     return;
0343   }
0344   //std::cout << " InOutConversionSeedFinder::findSeeds startingState R "<< startingState.globalPosition().perp() << " Z " << startingState.globalPosition().z() << " phi " <<  startingState.globalPosition().phi() <<  " position " << startingState.globalPosition() << "\n";
0345   //std::cout << " InOutConversionSeedFinder::findSeeds Initial FTS charge " << fts.charge() << " curvature " <<  transverseCurvature << "\n";
0346   //std::cout << " InOutConversionSeedFinder::findSeeds Initial FTS parameters " << fts <<  "\n";
0347 
0348   //float dphi = 0.01;
0349   float dphi = 0.03;
0350   float zrange = 5.;
0351   for (unsigned int ilayer = startingLayer; ilayer <= startingLayer + 1 && (ilayer < allLayers.size() - 2); ++ilayer) {
0352     const DetLayer* layer = allLayers[ilayer];
0353 
0354     ///// debug
0355     //    if ( layer->location() == GeomDetEnumerators::barrel ) {const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
0356     // //std::cout << "InOutConversionSeedFinder::findSeeds  ****  Barrel on layer " << ilayer  << " R= " << barrelLayer->specificSurface().radius() <<  "\n";
0357     // } else {
0358     // const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
0359     ////std::cout << "InOutConversionSeedFinder::findSeeds  ****  Forw on layer " << ilayer  << " Z= " << forwardLayer->specificSurface().position().z() <<  "\n";
0360     // }
0361     // // end debug
0362 
0363     MeasurementEstimator* newEstimator = nullptr;
0364     if (layer->location() == GeomDetEnumerators::barrel) {
0365       //std::cout << "InOutConversionSeedFinder::findSeeds Barrel ilayer " << ilayer <<  "\n";
0366       newEstimator = new ConversionBarrelEstimator(-dphi, dphi, -zrange, zrange);
0367     } else {
0368       //std::cout << "InOutConversionSeedFinder::findSeeds Forward  ilayer " << ilayer <<  "\n";
0369       newEstimator = new ConversionForwardEstimator(-dphi, dphi, 15.);
0370     }
0371 
0372     theFirstMeasurements_.clear();
0373     // Get measurements compatible with the FTS and Estimator
0374     TSOS tsos(fts, layer->surface());
0375 
0376     //std::cout << "InOutConversionSeedFinder::findSeed propagationDirection " << int(thePropagatorAlongMomentum_->propagationDirection() ) << "\n";
0377     /// Rememeber that this alwyas give back at least one dummy-innvalid it which prevents from everything getting stopped
0378     LayerMeasurements theLayerMeasurements_(*this->getMeasurementTracker(), *theTrackerData_);
0379 
0380     theFirstMeasurements_ =
0381         theLayerMeasurements_.measurements(*layer, tsos, *thePropagatorAlongMomentum_, *newEstimator);
0382 
0383     delete newEstimator;
0384     //std::cout <<  "InOutConversionSeedFinder::findSeeds  Found " << theFirstMeasurements_.size() << " first hits" << "\n";
0385 
0386     if (theFirstMeasurements_.size() ==
0387         1) {  // only dummy hit found: start finding the seed from the innermost hit of the OutIn track
0388 
0389       GlobalPoint bcPos((theSecondBC_.position()).x(), (theSecondBC_.position()).y(), (theSecondBC_.position()).z());
0390       GlobalVector dir = startingState.globalDirection();
0391       GlobalPoint back1mm = myPointer->recHit()->globalPosition();
0392 
0393       back1mm -= dir.unit() * 0.1;
0394       //std::cout << " InOutConversionSeedFinder:::findSeeds going to make the helix using back1mm " << back1mm << "\n";
0395       ConversionFastHelix helix(bcPos, myPointer->recHit()->globalPosition(), back1mm, &(*theMF_));
0396 
0397       helix.stateAtVertex();
0398       //std::cout << " InOutConversionSeedFinder:::findSeeds helix status " <<helix.isValid() << std::endl;
0399       if (!helix.isValid())
0400         continue;
0401 
0402       track2InitialMomentum_ = helix.stateAtVertex().momentum();
0403 
0404       // Make a new FTS
0405       FreeTrajectoryState newfts(GlobalTrajectoryParameters(myPointer->recHit()->globalPosition(),
0406                                                             startingState.globalDirection(),
0407                                                             helix.stateAtVertex().transverseCurvature(),
0408                                                             0,
0409                                                             &(*theMF_)),
0410                                  CurvilinearTrajectoryError(m));
0411 
0412       completeSeed(*myPointer, newfts, thePropagatorAlongMomentum_, ilayer + 1);
0413       completeSeed(*myPointer, newfts, thePropagatorAlongMomentum_, ilayer + 2);
0414 
0415     } else {
0416       //Loop over compatible hits
0417       int mea = 0;
0418       for (std::vector<TrajectoryMeasurement>::iterator tmItr = theFirstMeasurements_.begin();
0419            tmItr != theFirstMeasurements_.end();
0420            ++tmItr) {
0421         mea++;
0422 
0423         if (tmItr->recHit()->isValid()) {
0424           // Make a new helix as in fillClusterSeeds() but using the hit position
0425           //std::cout << "InOutConversionSeedFinder::findSeeds hit  R " << tmItr->recHit()->globalPosition().perp() << " Z " <<  tmItr->recHit()->globalPosition().z() << " " <<  tmItr->recHit()->globalPosition() << "\n";
0426           GlobalPoint bcPos(
0427               (theSecondBC_.position()).x(), (theSecondBC_.position()).y(), (theSecondBC_.position()).z());
0428           GlobalVector dir = startingState.globalDirection();
0429           GlobalPoint back1mm = tmItr->recHit()->globalPosition();
0430 
0431           back1mm -= dir.unit() * 0.1;
0432           //std::cout << " InOutConversionSeedFinder:::findSeeds going to make the helix using back1mm " << back1mm << "\n";
0433           ConversionFastHelix helix(bcPos, tmItr->recHit()->globalPosition(), back1mm, &(*theMF_));
0434 
0435           helix.stateAtVertex();
0436           //std::cout << " InOutConversionSeedFinder:::findSeeds helix status " <<helix.isValid() << std::endl;
0437           if (!helix.isValid())
0438             continue;
0439 
0440           track2InitialMomentum_ = helix.stateAtVertex().momentum();
0441 
0442           //std::cout << "InOutConversionSeedFinder::findSeeds Updated estimatedPt = " << helix.stateAtVertex().momentum().perp()  << " curvature "  << helix.stateAtVertex().transverseCurvature() << "\n";
0443           //     << ", bcet = " << theBc->Et()
0444           //     << ", estimatedPt/bcet = " << estimatedPt/theBc->Et() << endl;
0445 
0446           // Make a new FTS
0447           FreeTrajectoryState newfts(GlobalTrajectoryParameters(tmItr->recHit()->globalPosition(),
0448                                                                 startingState.globalDirection(),
0449                                                                 helix.stateAtVertex().transverseCurvature(),
0450                                                                 0,
0451                                                                 &(*theMF_)),
0452                                      CurvilinearTrajectoryError(m));
0453 
0454           //std::cout <<  "InOutConversionSeedFinder::findSeeds  new FTS charge " << newfts.charge() << "\n";
0455 
0456           /*
0457       // Soome diagnostic output
0458       // may be useful - comparission of the basic cluster position 
0459       // with the ecal impact position of the track
0460       TrajectoryStateOnSurface stateAtECAL
0461       = forwardPropagator.propagate(newfts, ECALSurfaces::barrel());
0462       if (!stateAtECAL.isValid() || abs(stateAtECAL.globalPosition().eta())>1.479) {
0463       if (startingState.globalDirection().eta() > 0.) {
0464       stateAtECAL = forwardPropagator.propagate(newfts, 
0465       ECALSurfaces::positiveEtaEndcap());
0466       } else {
0467       stateAtECAL = forwardPropagator.propagate(newfts, 
0468       ECALSurfaces::negativeEtaEndcap());
0469       }
0470       }
0471       GlobalPoint ecalImpactPosition = stateAtECAL.isValid() ? stateAtECAL.globalPosition() : GlobalPoint(0.,0.,0.);
0472       cout << "Projected fts positon at ECAL surface: " << 
0473       ecalImpactPosition << " bc position: " << theBc->Position() << endl;
0474       */
0475 
0476           completeSeed(*tmItr, newfts, thePropagatorAlongMomentum_, ilayer + 1);
0477           completeSeed(*tmItr, newfts, thePropagatorAlongMomentum_, ilayer + 2);
0478         }
0479       }
0480     }
0481   }
0482 }
0483 
0484 void InOutConversionSeedFinder::completeSeed(const TrajectoryMeasurement& m1,
0485                                              const FreeTrajectoryState& fts,
0486                                              const Propagator* propagator,
0487                                              int ilayer) {
0488   //std::cout<<  "InOutConversionSeedFinder::completeSeed ilayer " << ilayer <<  "\n";
0489   // A seed is made from 2 Trajectory Measuremennts.  The 1st is the input
0490   // argument m1.  This routine looks for the 2nd measurement in layer ilayer
0491   // Begin by making a new much stricter MeasurementEstimator based on the
0492   // position errors of the 1st hit.
0493   printLayer(ilayer);
0494 
0495   MeasurementEstimator* newEstimator;
0496   std::vector<const DetLayer*> allLayers = layerList();
0497   const DetLayer* layer = allLayers[ilayer];
0498 
0499   ///// debug
0500   //  if ( layer->location() == GeomDetEnumerators::barrel ) {const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
0501   // //std::cout << "InOutConversionSeedFinder::completeSeed  ****  Barrel on layer " << ilayer  << " R= " << barrelLayer->specificSurface().radius() <<  "\n";
0502   // } else {
0503   // const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
0504   //  //std::cout << "InOutConversionSeedFinder::completeSeed ****  Forw on layer " << ilayer  << " Z= " << forwardLayer->specificSurface().position().z() <<  "\n";
0505   ///  }
0506   //// end debug
0507 
0508   if (layer->location() == GeomDetEnumerators::barrel) {
0509     float dz = sqrt(the2ndHitdznSigma_ * the2ndHitdznSigma_ * m1.recHit()->globalPositionError().czz() +
0510                     the2ndHitdzConst_ * the2ndHitdzConst_);
0511     newEstimator = new ConversionBarrelEstimator(-the2ndHitdphi_, the2ndHitdphi_, -dz, dz);
0512 
0513   } else {
0514     float m1dr = sqrt(m1.recHit()->localPositionError().yy());
0515     float dr = sqrt(the2ndHitdznSigma_ * the2ndHitdznSigma_ * m1dr * m1dr + the2ndHitdzConst_ * the2ndHitdznSigma_);
0516 
0517     newEstimator = new ConversionForwardEstimator(-the2ndHitdphi_, the2ndHitdphi_, dr);
0518   }
0519 
0520   //std::cout << "InOutConversionSeedFinder::completeSeed fts For the TSOS " << fts << "\n";
0521 
0522   TSOS tsos(fts, layer->surface());
0523 
0524   if (!tsos.isValid()) {
0525     //std::cout  << "InOutConversionSeedFinder::completeSeed TSOS is not valid " <<  "\n";
0526   }
0527 
0528   //std::cout << "InOutConversionSeedFinder::completeSeed TSOS " << tsos << "\n";
0529   //std::cout << "InOutConversionSeedFinder::completeSeed propagationDirection  " << int(propagator->propagationDirection() ) << "\n";
0530   //std::cout << "InOutConversionSeedFinder::completeSeed pointer to estimator " << newEstimator << "\n";
0531 
0532   LayerMeasurements theLayerMeasurements_(*this->getMeasurementTracker(), *theTrackerData_);
0533   std::vector<TrajectoryMeasurement> measurements =
0534       theLayerMeasurements_.measurements(*layer, tsos, *propagator, *newEstimator);
0535   //std::cout << "InOutConversionSeedFinder::completeSeed Found " << measurements.size() << " second hits " <<  "\n";
0536   delete newEstimator;
0537 
0538   for (unsigned int i = 0; i < measurements.size(); ++i) {
0539     if (measurements[i].recHit()->isValid()) {
0540       createSeed(m1, measurements[i]);
0541     }
0542   }
0543 }
0544 
0545 void InOutConversionSeedFinder::createSeed(const TrajectoryMeasurement& m1, const TrajectoryMeasurement& m2) {
0546   //std::cout << "InOutConversionSeedFinder::createSeed " << "\n";
0547 
0548   if (m1.predictedState().isValid()) {
0549     GlobalTrajectoryParameters newgtp(m1.recHit()->globalPosition(), track2InitialMomentum_, track2Charge_, &(*theMF_));
0550     CurvilinearTrajectoryError errors = m1.predictedState().curvilinearError();
0551     FreeTrajectoryState fts(newgtp, errors);
0552 
0553     TrajectoryStateOnSurface state1 = thePropagatorAlongMomentum_->propagate(fts, m1.recHit()->det()->surface());
0554 
0555     /*
0556    //std::cout << "hit surface " <<  m1.recHit()->det()->surface().position() << "\n";
0557    //std::cout << "prop to " << typeid( m1.recHit()->det()->surface() ).name() <<"\n";
0558    //std::cout << "prop to first hit " << state1 << "\n"; 
0559    //std::cout << "update to " <<  m1.recHit()->globalPosition() << "\n";
0560   */
0561 
0562     if (state1.isValid()) {
0563       TrajectoryStateOnSurface updatedState1 = theUpdator_.update(state1, *m1.recHit());
0564 
0565       if (updatedState1.isValid()) {
0566         TrajectoryStateOnSurface state2 =
0567             thePropagatorAlongMomentum_->propagate(*updatedState1.freeTrajectoryState(), m2.recHit()->det()->surface());
0568 
0569         if (state2.isValid()) {
0570           TrajectoryStateOnSurface updatedState2 = theUpdator_.update(state2, *m2.recHit());
0571           TrajectoryMeasurement meas1(state1, updatedState1, m1.recHit(), m1.estimate(), m1.layer());
0572           TrajectoryMeasurement meas2(state2, updatedState2, m2.recHit(), m2.estimate(), m2.layer());
0573 
0574           edm::OwnVector<TrackingRecHit> myHits;
0575           myHits.push_back(meas1.recHit()->hit()->clone());
0576           myHits.push_back(meas2.recHit()->hit()->clone());
0577 
0578           //std::cout << "InOutConversionSeedFinder::createSeed new seed " << "\n";
0579           if (nSeedsPerInputTrack_ >= maxNumberOfInOutSeedsPerInputTrack_)
0580             return;
0581 
0582           PTrajectoryStateOnDet const& ptsod =
0583               trajectoryStateTransform::persistentState(state2, meas2.recHit()->hit()->geographicalId().rawId());
0584           //std::cout << "  InOutConversionSeedFinder::createSeed New seed parameters " << state2 << "\n";
0585 
0586           theSeeds_.push_back(TrajectorySeed(ptsod, myHits, alongMomentum));
0587           nSeedsPerInputTrack_++;
0588 
0589           //std::cout << "InOutConversionSeedFinder::createSeed New seed hit 1 R " << m1.recHit()->globalPosition().perp() << "\n";
0590           //std::cout << "InOutConversionSeedFinder::createSeed New seed hit 2 R " << m2.recHit()->globalPosition().perp() << "\n";
0591         }
0592       }
0593     }
0594   }
0595 }