![]() |
|
|||
File indexing completed on 2023-05-08 23:19:12
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 for (std::vector<TrajectoryMeasurement>::iterator tmItr = theFirstMeasurements_.begin(); 0418 tmItr != theFirstMeasurements_.end(); 0419 ++tmItr) { 0420 if (tmItr->recHit()->isValid()) { 0421 // Make a new helix as in fillClusterSeeds() but using the hit position 0422 //std::cout << "InOutConversionSeedFinder::findSeeds hit R " << tmItr->recHit()->globalPosition().perp() << " Z " << tmItr->recHit()->globalPosition().z() << " " << tmItr->recHit()->globalPosition() << "\n"; 0423 GlobalPoint bcPos( 0424 (theSecondBC_.position()).x(), (theSecondBC_.position()).y(), (theSecondBC_.position()).z()); 0425 GlobalVector dir = startingState.globalDirection(); 0426 GlobalPoint back1mm = tmItr->recHit()->globalPosition(); 0427 0428 back1mm -= dir.unit() * 0.1; 0429 //std::cout << " InOutConversionSeedFinder:::findSeeds going to make the helix using back1mm " << back1mm << "\n"; 0430 ConversionFastHelix helix(bcPos, tmItr->recHit()->globalPosition(), back1mm, &(*theMF_)); 0431 0432 helix.stateAtVertex(); 0433 //std::cout << " InOutConversionSeedFinder:::findSeeds helix status " <<helix.isValid() << std::endl; 0434 if (!helix.isValid()) 0435 continue; 0436 0437 track2InitialMomentum_ = helix.stateAtVertex().momentum(); 0438 0439 //std::cout << "InOutConversionSeedFinder::findSeeds Updated estimatedPt = " << helix.stateAtVertex().momentum().perp() << " curvature " << helix.stateAtVertex().transverseCurvature() << "\n"; 0440 // << ", bcet = " << theBc->Et() 0441 // << ", estimatedPt/bcet = " << estimatedPt/theBc->Et() << endl; 0442 0443 // Make a new FTS 0444 FreeTrajectoryState newfts(GlobalTrajectoryParameters(tmItr->recHit()->globalPosition(), 0445 startingState.globalDirection(), 0446 helix.stateAtVertex().transverseCurvature(), 0447 0, 0448 &(*theMF_)), 0449 CurvilinearTrajectoryError(m)); 0450 0451 //std::cout << "InOutConversionSeedFinder::findSeeds new FTS charge " << newfts.charge() << "\n"; 0452 0453 /* 0454 // Soome diagnostic output 0455 // may be useful - comparission of the basic cluster position 0456 // with the ecal impact position of the track 0457 TrajectoryStateOnSurface stateAtECAL 0458 = forwardPropagator.propagate(newfts, ECALSurfaces::barrel()); 0459 if (!stateAtECAL.isValid() || abs(stateAtECAL.globalPosition().eta())>1.479) { 0460 if (startingState.globalDirection().eta() > 0.) { 0461 stateAtECAL = forwardPropagator.propagate(newfts, 0462 ECALSurfaces::positiveEtaEndcap()); 0463 } else { 0464 stateAtECAL = forwardPropagator.propagate(newfts, 0465 ECALSurfaces::negativeEtaEndcap()); 0466 } 0467 } 0468 GlobalPoint ecalImpactPosition = stateAtECAL.isValid() ? stateAtECAL.globalPosition() : GlobalPoint(0.,0.,0.); 0469 cout << "Projected fts positon at ECAL surface: " << 0470 ecalImpactPosition << " bc position: " << theBc->Position() << endl; 0471 */ 0472 0473 completeSeed(*tmItr, newfts, thePropagatorAlongMomentum_, ilayer + 1); 0474 completeSeed(*tmItr, newfts, thePropagatorAlongMomentum_, ilayer + 2); 0475 } 0476 } 0477 } 0478 } 0479 } 0480 0481 void InOutConversionSeedFinder::completeSeed(const TrajectoryMeasurement& m1, 0482 const FreeTrajectoryState& fts, 0483 const Propagator* propagator, 0484 int ilayer) { 0485 //std::cout<< "InOutConversionSeedFinder::completeSeed ilayer " << ilayer << "\n"; 0486 // A seed is made from 2 Trajectory Measuremennts. The 1st is the input 0487 // argument m1. This routine looks for the 2nd measurement in layer ilayer 0488 // Begin by making a new much stricter MeasurementEstimator based on the 0489 // position errors of the 1st hit. 0490 printLayer(ilayer); 0491 0492 MeasurementEstimator* newEstimator; 0493 std::vector<const DetLayer*> allLayers = layerList(); 0494 const DetLayer* layer = allLayers[ilayer]; 0495 0496 ///// debug 0497 // if ( layer->location() == GeomDetEnumerators::barrel ) {const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer); 0498 // //std::cout << "InOutConversionSeedFinder::completeSeed **** Barrel on layer " << ilayer << " R= " << barrelLayer->specificSurface().radius() << "\n"; 0499 // } else { 0500 // const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer); 0501 // //std::cout << "InOutConversionSeedFinder::completeSeed **** Forw on layer " << ilayer << " Z= " << forwardLayer->specificSurface().position().z() << "\n"; 0502 /// } 0503 //// end debug 0504 0505 if (layer->location() == GeomDetEnumerators::barrel) { 0506 float dz = sqrt(the2ndHitdznSigma_ * the2ndHitdznSigma_ * m1.recHit()->globalPositionError().czz() + 0507 the2ndHitdzConst_ * the2ndHitdzConst_); 0508 newEstimator = new ConversionBarrelEstimator(-the2ndHitdphi_, the2ndHitdphi_, -dz, dz); 0509 0510 } else { 0511 float m1dr = sqrt(m1.recHit()->localPositionError().yy()); 0512 float dr = sqrt(the2ndHitdznSigma_ * the2ndHitdznSigma_ * m1dr * m1dr + the2ndHitdzConst_ * the2ndHitdznSigma_); 0513 0514 newEstimator = new ConversionForwardEstimator(-the2ndHitdphi_, the2ndHitdphi_, dr); 0515 } 0516 0517 //std::cout << "InOutConversionSeedFinder::completeSeed fts For the TSOS " << fts << "\n"; 0518 0519 TSOS tsos(fts, layer->surface()); 0520 0521 if (!tsos.isValid()) { 0522 //std::cout << "InOutConversionSeedFinder::completeSeed TSOS is not valid " << "\n"; 0523 } 0524 0525 //std::cout << "InOutConversionSeedFinder::completeSeed TSOS " << tsos << "\n"; 0526 //std::cout << "InOutConversionSeedFinder::completeSeed propagationDirection " << int(propagator->propagationDirection() ) << "\n"; 0527 //std::cout << "InOutConversionSeedFinder::completeSeed pointer to estimator " << newEstimator << "\n"; 0528 0529 LayerMeasurements theLayerMeasurements_(*this->getMeasurementTracker(), *theTrackerData_); 0530 std::vector<TrajectoryMeasurement> measurements = 0531 theLayerMeasurements_.measurements(*layer, tsos, *propagator, *newEstimator); 0532 //std::cout << "InOutConversionSeedFinder::completeSeed Found " << measurements.size() << " second hits " << "\n"; 0533 delete newEstimator; 0534 0535 for (unsigned int i = 0; i < measurements.size(); ++i) { 0536 if (measurements[i].recHit()->isValid()) { 0537 createSeed(m1, measurements[i]); 0538 } 0539 } 0540 } 0541 0542 void InOutConversionSeedFinder::createSeed(const TrajectoryMeasurement& m1, const TrajectoryMeasurement& m2) { 0543 //std::cout << "InOutConversionSeedFinder::createSeed " << "\n"; 0544 0545 if (m1.predictedState().isValid()) { 0546 GlobalTrajectoryParameters newgtp(m1.recHit()->globalPosition(), track2InitialMomentum_, track2Charge_, &(*theMF_)); 0547 CurvilinearTrajectoryError errors = m1.predictedState().curvilinearError(); 0548 FreeTrajectoryState fts(newgtp, errors); 0549 0550 TrajectoryStateOnSurface state1 = thePropagatorAlongMomentum_->propagate(fts, m1.recHit()->det()->surface()); 0551 0552 /* 0553 //std::cout << "hit surface " << m1.recHit()->det()->surface().position() << "\n"; 0554 //std::cout << "prop to " << typeid( m1.recHit()->det()->surface() ).name() <<"\n"; 0555 //std::cout << "prop to first hit " << state1 << "\n"; 0556 //std::cout << "update to " << m1.recHit()->globalPosition() << "\n"; 0557 */ 0558 0559 if (state1.isValid()) { 0560 TrajectoryStateOnSurface updatedState1 = theUpdator_.update(state1, *m1.recHit()); 0561 0562 if (updatedState1.isValid()) { 0563 TrajectoryStateOnSurface state2 = 0564 thePropagatorAlongMomentum_->propagate(*updatedState1.freeTrajectoryState(), m2.recHit()->det()->surface()); 0565 0566 if (state2.isValid()) { 0567 TrajectoryStateOnSurface updatedState2 = theUpdator_.update(state2, *m2.recHit()); 0568 TrajectoryMeasurement meas1(state1, updatedState1, m1.recHit(), m1.estimate(), m1.layer()); 0569 TrajectoryMeasurement meas2(state2, updatedState2, m2.recHit(), m2.estimate(), m2.layer()); 0570 0571 edm::OwnVector<TrackingRecHit> myHits; 0572 myHits.push_back(meas1.recHit()->hit()->clone()); 0573 myHits.push_back(meas2.recHit()->hit()->clone()); 0574 0575 //std::cout << "InOutConversionSeedFinder::createSeed new seed " << "\n"; 0576 if (nSeedsPerInputTrack_ >= maxNumberOfInOutSeedsPerInputTrack_) 0577 return; 0578 0579 PTrajectoryStateOnDet const& ptsod = 0580 trajectoryStateTransform::persistentState(state2, meas2.recHit()->hit()->geographicalId().rawId()); 0581 //std::cout << " InOutConversionSeedFinder::createSeed New seed parameters " << state2 << "\n"; 0582 0583 theSeeds_.push_back(TrajectorySeed(ptsod, myHits, alongMomentum)); 0584 nSeedsPerInputTrack_++; 0585 0586 //std::cout << "InOutConversionSeedFinder::createSeed New seed hit 1 R " << m1.recHit()->globalPosition().perp() << "\n"; 0587 //std::cout << "InOutConversionSeedFinder::createSeed New seed hit 2 R " << m2.recHit()->globalPosition().perp() << "\n"; 0588 } 0589 } 0590 } 0591 } 0592 }
[ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
This page was automatically generated by the 2.2.1 LXR engine. The LXR team |
![]() ![]() |