CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
InOutConversionSeedFinder.cc
Go to the documentation of this file.
1 #include "RecoEgamma/EgammaPhotonAlgos/interface/InOutConversionSeedFinder.h"
2 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionBarrelEstimator.h"
3 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionForwardEstimator.h"
4 
5 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionFastHelix.h"
6 
7 // Field
8 //
9 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
10 // Geometry
11 //
12 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
13 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
14 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
15 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
16 
17 //
18 //
19 //
20 
21 InOutConversionSeedFinder::InOutConversionSeedFinder( const edm::ParameterSet& conf,edm::ConsumesCollector && iC ):
22  ConversionSeedFinder( conf,iC ), conf_(conf)
23 {
24 
25 
26  // std::cout << " InOutConversionSeedFinder CTOR " << "\n";
27 
28  maxNumberOfInOutSeedsPerInputTrack_ = conf_.getParameter<int>("maxNumOfSeedsInOut");
29  //the2ndHitdphi_ = 0.008;
30  the2ndHitdphi_ = 0.01;
31  the2ndHitdzConst_ = 5.;
32  the2ndHitdznSigma_ = 2.;
33 
34 
35 }
36 
37 
38 
39 InOutConversionSeedFinder::~InOutConversionSeedFinder() {
40  //std::cout << " InOutConversionSeedFinder DTOR " << "\n";
41 }
42 
43 
44 
45 void InOutConversionSeedFinder::makeSeeds( const edm::Handle<edm::View<reco::CaloCluster> > & allBC ) const {
46 
47 
48  //std::cout << " InOutConversionSeedFinder::makeSeeds() " << "\n";
49  theSeeds_.clear();
50  //std::cout << " Check Calo cluster collection size " << allBC->size() << "\n";
51  bcCollection_= allBC;
52 
53 
54  findLayers();
55 
56 
57  fillClusterSeeds();
58  //std::cout << "Built vector of seeds of size " << theSeeds_.size() << "\n" ;
59 
60 
61  theOutInTracks_.clear();
62  inputTracks_.clear();
63  theFirstMeasurements_.clear();
64 
65 }
66 
67 
68 void InOutConversionSeedFinder::fillClusterSeeds() const {
69 
70  std::vector<Trajectory>::const_iterator outInTrackItr;
71 
72  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds outInTracks_.size " << theOutInTracks_.size() << "\n";
73  //Start looking for seeds for both of the 2 best tracks from the inward tracking
74 
76  /*
77  for(outInTrackItr = theOutInTracks_.begin(); outInTrackItr != theOutInTracks_.end(); ++outInTrackItr) {
78 
79 
80  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds out in input track hits " << (*outInTrackItr).foundHits() << "\n";
81  DetId tmpId = DetId( (*outInTrackItr).seed().startingState().detId());
82  const GeomDet* tmpDet = this->getMeasurementTracker()->geomTracker()->idToDet( tmpId );
83  GlobalVector gv = tmpDet->surface().toGlobal( (*outInTrackItr).seed().startingState().parameters().momentum() );
84 
85 
86  //std::cout << " InOutConversionSeedFinder::fillClusterSeed was built from seed position " <<gv << " charge " << (*outInTrackItr).seed().startingState().parameters().charge() << "\n";
87 
88  Trajectory::DataContainer m= outInTrackItr->measurements();
89  int nHit=0;
90  for (Trajectory::DataContainer::iterator itm = m.begin(); itm != m.end(); ++itm) {
91  if ( itm->recHit()->isValid() ) {
92  nHit++;
93  //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";
94  }
95 
96  }
97 
98  }
99 
100  */
101 
102 
103  //Start looking for seeds for both of the 2 best tracks from the inward tracking
104  for(outInTrackItr = theOutInTracks_.begin(); outInTrackItr != theOutInTracks_.end(); ++outInTrackItr) {
105  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds out in input track hits " << (*outInTrackItr).foundHits() << "\n";
106  nSeedsPerInputTrack_=0;
107 
108  //Find the first valid hit of the track
109  // Measurements are ordered according to the direction in which the trajectories were built
110  std::vector<TrajectoryMeasurement> measurements = (*outInTrackItr).measurements();
111 
112  std::vector<const DetLayer*> allLayers=layerList();
113 
114  //std::cout << " InOutConversionSeedFinder::fill clusterSeed allLayers.size " << allLayers.size() << "\n";
115  for(unsigned int i = 0; i < allLayers.size(); ++i) {
116  //std::cout << " allLayers " << allLayers[i] << "\n";
117  printLayer(i);
118  }
119 
120 
121 
122  std::vector<const DetLayer*> myLayers;
123  myLayers.clear();
124  std::vector<TrajectoryMeasurement>::reverse_iterator measurementItr;
125  std::vector<TrajectoryMeasurement*> myItr;
126  // TrajectoryMeasurement* myPointer=0;
127  myPointer=0;
128  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds measurements.size " << measurements.size() <<"\n";
129 
130  for(measurementItr = measurements.rbegin() ; measurementItr != measurements.rend(); ++measurementItr) {
131 
132 
133  if( (*measurementItr).recHit()->isValid()) {
134 
135  //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";
136 
137 
138  myLayers.push_back( measurementItr->layer() ) ;
139  myItr.push_back( &(*measurementItr) );
140 
141 
142  }
143  }
144 
145 
146 
147  //std::cout << " InOutConversionSeedFinder::fillClusterSeed myLayers.size " << myLayers.size() << "\n";
148  // for( unsigned int i = 0; i < myLayers.size(); ++i) {
150  // }
151 
152 
153  if ( myItr.size()==0 ) {
154  //std::cout << "HORRENDOUS ERROR! No meas on track!" << "\n";
155  }
156  unsigned int ilayer;
157  for(ilayer = 0; ilayer < allLayers.size(); ++ilayer) {
158  //std::cout << " allLayers in the search loop " << allLayers[ilayer] << " " << myLayers[0] << "\n";
159  if ( allLayers[ilayer] == myLayers[0]) {
160 
161  myPointer=myItr[0];
162 
163  //std::cout << " allLayers in the search loop allLayers[ilayer] == myLayers[0]) " << allLayers[ilayer] << " " << myLayers[0] << " myPointer " << myPointer << "\n";
164 
165  //std::cout << "Layer " << ilayer << " contains the first valid measurement " << "\n";
166  printLayer(ilayer);
167 
168  if ( (myLayers[0])->location() == GeomDetEnumerators::barrel ) {
169  // const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(myLayers[0]);
170  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds **** firstHit found in Barrel on layer " << ilayer << " R= " << barrelLayer->specificSurface().radius() << "\n";
171  } else {
172  //const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(myLayers[0]);
173  //std::cout << " InOutwardConversionSeedFinder::fillClusterSeeds **** firstHit found in Forw on layer " << ilayer << " Z= " << forwardLayer->specificSurface().position().z() << "\n";
174  }
175 
176 
177  break;
178 
179  } else if ( allLayers[ilayer] == myLayers[1] ) {
180  myPointer=myItr[1];
181 
182  //std::cout << " allLayers in the search loop allLayers[ilayer] == myLayers[1]) " << allLayers[ilayer] << " " << myLayers[1] << " myPointer " << myPointer << "\n";
183 
184  //std::cout << "Layer " << ilayer << " contains the first innermost valid measurement " << "\n";
185  if ( (myLayers[1])->location() == GeomDetEnumerators::barrel ) {
186  // const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(myLayers[1]);
187  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds **** 2ndHit found in Barrel on layer " << ilayer << " R= " << barrelLayer->specificSurface().radius() << "\n";
188  } else {
189  //const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(myLayers[1]);
190  //std::cout << " InOutwardConversionSeedFinder::fillClusterSeeds **** 2ndHitfound on forw layer " << ilayer << " Z= " << forwardLayer->specificSurface().position().z() << "\n";
191  }
192 
193 
194 
195  break;
196 
197  }
198  }
199 
200 
201 
202  if(ilayer == allLayers.size()) {
203  //std::cout << "InOutConversionSeedFinder::fillClusterSeeds ERROR could not find layer on list" << "\n";
204  return;
205  }
206 
207  //PropagatorWithMaterial reversePropagator(oppositeToMomentum, 0.000511, &(*theMF_) );
208  const FreeTrajectoryState * fts = myPointer->updatedState().freeTrajectoryState();
209 
210  //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";
211 
212 
213  while (ilayer > 0) {
214 
215  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds looking for 2nd seed from layer " << ilayer << "\n";
216 
217  // if ( (allLayers[ilayer])->location() == GeomDetEnumerators::barrel ) {const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(allLayers[ilayer]);
218  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds **** Barrel on layer " << ilayer << " R= " << barrelLayer->specificSurface().radius() << "\n";
219  // } else {
220  //const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(allLayers[ilayer]);
221  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds **** Forw on layer " << ilayer << " Z= " << forwardLayer->specificSurface().position().z() << "\n";
222  // }
223 
224 
225  const DetLayer * previousLayer = allLayers[ilayer];
226  TrajectoryStateOnSurface stateAtPreviousLayer;
227  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds previousLayer->surface() position before " <<allLayers[ilayer] << " " << previousLayer->surface().position() << " layer location " << previousLayer->location() << "\n";
228  // Propagate to the previous layer
229  // The present layer is actually included in the loop so that a partner can be searched for
230  // Applying the propagator to the same layer does not do any harm. It simply does nothing
231 
232  // const Propagator& newProp=thePropagatorOppositeToMomentum_;
233  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds reversepropagator direction " << thePropagatorOppositeToMomentum_->propagationDirection() << "\n";
234  if (ilayer-1>0) {
235 
236  if ( allLayers[ilayer] == myLayers[0] ) {
237  //std::cout << " innermost hit R " << myPointer->recHit()->globalPosition().perp() << " Z " << myPointer->recHit()->globalPosition().z() << " phi " <<myPointer->recHit()->globalPosition().phi() << "\n";
238  //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";
239 
240  stateAtPreviousLayer= thePropagatorOppositeToMomentum_->propagate(*fts, theTrackerGeom_->idToDet( myPointer->recHit() ->geographicalId())->surface() );
241 
242  } else {
243 
244  stateAtPreviousLayer= thePropagatorOppositeToMomentum_->propagate(*fts, previousLayer->surface() );
245  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds previousLayer->surface() position after " << previousLayer->surface().position() << " layer location " << previousLayer->location() << "\n";
246 
247  }
248 
249  } else if ( ilayer-1==0) {
250 
251 
254 
255  //stateAtPreviousLayer= thePropagatorOppositeToMomentum_->propagate(*fts, theTrackerGeom_->idToDet( myPointer->recHit() ->geographicalId())->surface() );
256  stateAtPreviousLayer= thePropagatorOppositeToMomentum_->propagate(*fts, previousLayer->surface() );
257 
258  }
259 
260  if(!stateAtPreviousLayer.isValid()) {
261  //std::cout << "InOutConversionSeedFinder::fillClusterSeeds ERROR:could not propagate back to layer " << ilayer << "\n";
263  } else {
264  //std::cout << "InOutConversionSeedFinder::fillClusterSeeds stateAtPreviousLayer is valid. Propagating back to layer " << ilayer << "\n";
265  //std::cout << "InOutConversionSeedFinder::fillClusterSeeds stateAtPreviousLayer R " << stateAtPreviousLayer.globalPosition().perp() << " Z " << stateAtPreviousLayer.globalPosition().z() << " phi " << stateAtPreviousLayer.globalPosition().phi() << "\n";
266 
267  startSeed(fts, stateAtPreviousLayer, -1, ilayer );
268 
269 
270  }
271 
272  --ilayer;
273 
274  }
275 
276  if ( ilayer == 0) {
277 
278 
279  // if ( (allLayers[ilayer])->location() == GeomDetEnumerators::barrel ) {const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(allLayers[ilayer]);
280  // //std::cout << " InOutConversionSeedFinder::fillClusterSeeds **** Barrel on layer " << ilayer << " R= " << barrelLayer->specificSurface().radius() << "\n";
281  // } else {
282  //const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(allLayers[ilayer]);
283  //std::cout << " InOutConversionSeedFinder::fillClusterSeeds **** Forw on layer " << ilayer << " Z= " << forwardLayer->specificSurface().position().z() << "\n";
284  // }
285  const DetLayer * previousLayer = allLayers[ilayer];
286  TrajectoryStateOnSurface stateAtPreviousLayer;
287  stateAtPreviousLayer= thePropagatorOppositeToMomentum_->propagate(*fts, previousLayer->surface() );
288 
289  if(!stateAtPreviousLayer.isValid()) {
290  //std::cout << "InOutConversionSeedFinder::fillClusterSeeds ERROR:could not propagate back to layer " << ilayer << "\n";
292  } else {
293  //std::cout << "InOutConversionSeedFinder::fillClusterSeeds stateAtPreviousLayer is valid. Propagating back to layer " << ilayer << "\n";
294  //std::cout << "InOutConversionSeedFinder::fillClusterSeeds stateAtPreviousLayer R " << stateAtPreviousLayer.globalPosition().perp() << " Z " << stateAtPreviousLayer.globalPosition().z() << " phi " << stateAtPreviousLayer.globalPosition().phi() << "\n";
295 
296  startSeed(fts, stateAtPreviousLayer, -1, ilayer );
297  }
298 
299 
300  }
301 
302 
303 
304 
305  } // End loop over Out In tracks
306 
307 
308 
309 }
310 
311 
312 
313 void InOutConversionSeedFinder::startSeed( const FreeTrajectoryState * fts, const TrajectoryStateOnSurface & stateAtPreviousLayer, int charge, int ilayer ) const {
314 
315  //std::cout << "InOutConversionSeedFinder::startSeed ilayer " << ilayer << "\n";
316  // Get a list of basic clusters that are consistent with a track
317  // starting at the assumed conversion point with opp. charge to the
318  // inward track. Loop over these basic clusters.
319  track2Charge_ = charge*fts->charge();
320  std::vector<const reco::CaloCluster*> bcVec;
321  //std::cout << "InOutConversionSeedFinder::startSeed charge assumed for the in-out track " << track2Charge_ << "\n";
322 
323  // Geom::Phi<float> theConvPhi( stateAtPreviousLayer.globalPosition().phi());
324  //std::cout << "InOutConversionSeedFinder::startSeed stateAtPreviousLayer phi " << stateAtPreviousLayer.globalPosition().phi() << " R " << stateAtPreviousLayer.globalPosition().perp() << " Z " << stateAtPreviousLayer.globalPosition().z() << "\n";
325 
326  bcVec = getSecondCaloClusters(stateAtPreviousLayer.globalPosition(),track2Charge_);
327 
328  std::vector<const reco::CaloCluster*>::iterator bcItr;
329  //std::cout << "InOutConversionSeedFinder::startSeed bcVec.size " << bcVec.size() << "\n";
330 
331  // debug
332  // for(bcItr = bcVec.begin(); bcItr != bcVec.end(); ++bcItr) {
333  // //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";
334  // }
335 
336 
337  for(bcItr = bcVec.begin(); bcItr != bcVec.end(); ++bcItr) {
338 
339  theSecondBC_ = **bcItr;
340  GlobalPoint bcPos((theSecondBC_.position()).x(),
341  (theSecondBC_.position()).y(),
342  (theSecondBC_.position()).z());
343 
344  //std::cout << "InOutConversionSeedFinder::startSeed for bc position x " << bcPos.x() << " y " << bcPos.y() << " z " << bcPos.z() << " eta " << bcPos.eta() << " phi " << bcPos.phi() << "\n";
345  GlobalVector dir = stateAtPreviousLayer.globalDirection();
346  GlobalPoint back1mm = stateAtPreviousLayer.globalPosition();
347  //std::cout << "InOutConversionSeedFinder::startSeed stateAtPreviousLayer.globalPosition() " << back1mm << "\n";
348 
349  back1mm -= dir.unit()*0.1;
350  //std::cout << " InOutConversionSeedFinder:::startSeed going to make the helix using back1mm " << back1mm <<"\n";
351  ConversionFastHelix helix(bcPos, stateAtPreviousLayer.globalPosition(), back1mm, &(*theMF_));
352  helix.stateAtVertex();
353 
354  //std::cout << " InOutConversionSeedFinder:::startSeed helix status " <<helix.isValid() << std::endl;
355  if ( !helix.isValid() ) continue;
356  findSeeds(stateAtPreviousLayer, helix.stateAtVertex().transverseCurvature(), ilayer);
357 
358 
359  }
360 
361 
362 
363 }
364 
365 
366 
367 std::vector<const reco::CaloCluster*> InOutConversionSeedFinder::getSecondCaloClusters(const GlobalPoint & conversionPosition, float charge) const {
368 
369 
370  std::vector<const reco::CaloCluster*> result;
371 
372  //std::cout << "InOutConversionSeedFinder::getSecondCaloClusters" << "\n";
373 
374  Geom::Phi<float> theConvPhi(conversionPosition.phi() );
375 
376  for (unsigned i = 0; i < bcCollection_->size(); ++i ) {
377 
378  Geom::Phi<float> theBcPhi( bcCollection_->ptrAt(i)->position().phi() );
379  //std::cout<< "InOutConversionSeedFinder::getSecondCaloClusters BC energy " << bcCollection_->ptrAt(i)->energy() << " Calo cluster phi " << theBcPhi << " " << bcCollection_->ptrAt(i)->position().phi()<< " theConvPhi " << theConvPhi << "\n";
380 
381  // Require phi of cluster to be consistent with the conversion
382  // position and the track charge
383 
384 
385  if (fabs(theBcPhi-theConvPhi ) < .5 &&
386  ((charge<0 && theBcPhi-theConvPhi >-.5) ||
387  (charge>0 && theBcPhi-theConvPhi <.5))){
389 
390  //result.push_back(&(*bcItr));
391 
392  result.push_back(&(*(bcCollection_->ptrAt(i)) ));
393 
394  }
395 
396 
397 
398 
399  }
400 
401 
402 
403  return result;
404 
405 
406 }
407 
408 
409 
410 void InOutConversionSeedFinder::findSeeds(const TrajectoryStateOnSurface & startingState,
411  float transverseCurvature,
412  unsigned int startingLayer) const {
413 
414 
415  std::vector<const DetLayer*> allLayers=layerList();
416  //std::cout << "InOutConversionSeedFinder::findSeeds starting forward propagation from startingLayer " << startingLayer << "\n";
417 
418 
419  // create error matrix
420  AlgebraicSymMatrix55 m = AlgebraicMatrixID();
421  m(0,0) = 0.1; m(1,1) = 0.0001 ; m(2,2) = 0.0001 ;
422  m(3,3) = 0.0001 ; m(4,4) = 0.001;
423 
424  // Make an FTS consistent with the start point, start direction and curvature
425  FreeTrajectoryState fts(GlobalTrajectoryParameters(startingState.globalPosition(),
426  startingState.globalDirection(),
427  double(transverseCurvature), 0, &(*theMF_) ),
428  CurvilinearTrajectoryError(m));
429  if (fts.momentum().mag2() == 0){
430  edm::LogWarning("FailedToInitiateSeeding")<< " initial FTS has a zero momentum, probably because of the zero field. ";
431  return;
432  }
433  //std::cout << " InOutConversionSeedFinder::findSeeds startingState R "<< startingState.globalPosition().perp() << " Z " << startingState.globalPosition().z() << " phi " << startingState.globalPosition().phi() << " position " << startingState.globalPosition() << "\n";
434  //std::cout << " InOutConversionSeedFinder::findSeeds Initial FTS charge " << fts.charge() << " curvature " << transverseCurvature << "\n";
435  //std::cout << " InOutConversionSeedFinder::findSeeds Initial FTS parameters " << fts << "\n";
436 
437 
438  //float dphi = 0.01;
439  float dphi = 0.03;
440  float zrange = 5.;
441  for( unsigned int ilayer = startingLayer; ilayer <= startingLayer+1 && (ilayer < allLayers.size()-2); ++ilayer) {
442  const DetLayer * layer = allLayers[ilayer];
443 
444 
445 
447  // if ( layer->location() == GeomDetEnumerators::barrel ) {const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
448  // //std::cout << "InOutConversionSeedFinder::findSeeds **** Barrel on layer " << ilayer << " R= " << barrelLayer->specificSurface().radius() << "\n";
449  // } else {
450  // const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
452  // }
453  // // end debug
454 
455 
456 
457  MeasurementEstimator * newEstimator=0;
458  if (layer->location() == GeomDetEnumerators::barrel ) {
459  //std::cout << "InOutConversionSeedFinder::findSeeds Barrel ilayer " << ilayer << "\n";
460  newEstimator = new ConversionBarrelEstimator(-dphi, dphi, -zrange, zrange);
461  }
462  else {
463  //std::cout << "InOutConversionSeedFinder::findSeeds Forward ilayer " << ilayer << "\n";
464  newEstimator = new ConversionForwardEstimator(-dphi, dphi, 15.);
465  }
466 
467 
468  theFirstMeasurements_.clear();
469  // Get measurements compatible with the FTS and Estimator
470  TSOS tsos(fts, layer->surface() );
471 
472  //std::cout << "InOutConversionSeedFinder::findSeed propagationDirection " << int(thePropagatorAlongMomentum_->propagationDirection() ) << "\n";
474  LayerMeasurements theLayerMeasurements_( *this->getMeasurementTracker(), *theTrackerData_ );
475 
476  theFirstMeasurements_ = theLayerMeasurements_.measurements( *layer, tsos, *thePropagatorAlongMomentum_, *newEstimator);
477 
478  delete newEstimator;
479  //std::cout << "InOutConversionSeedFinder::findSeeds Found " << theFirstMeasurements_.size() << " first hits" << "\n";
480 
481  if ( theFirstMeasurements_.size() == 1 ) { // only dummy hit found: start finding the seed from the innermost hit of the OutIn track
482 
483 
484  GlobalPoint bcPos((theSecondBC_.position()).x(),(theSecondBC_.position()).y(),(theSecondBC_.position()).z());
485  GlobalVector dir = startingState.globalDirection();
486  GlobalPoint back1mm = myPointer->recHit()->globalPosition();
487 
488  back1mm -= dir.unit()*0.1;
489  //std::cout << " InOutConversionSeedFinder:::findSeeds going to make the helix using back1mm " << back1mm << "\n";
490  ConversionFastHelix helix(bcPos, myPointer->recHit()->globalPosition(), back1mm, &(*theMF_));
491 
492  helix.stateAtVertex();
493  //std::cout << " InOutConversionSeedFinder:::findSeeds helix status " <<helix.isValid() << std::endl;
494  if ( !helix.isValid() ) continue;
495 
496  track2InitialMomentum_= helix.stateAtVertex().momentum();
497 
498  // Make a new FTS
499  FreeTrajectoryState newfts(GlobalTrajectoryParameters(
500  myPointer->recHit()->globalPosition(), startingState.globalDirection(),
501  helix.stateAtVertex().transverseCurvature(), 0, &(*theMF_)),
502  CurvilinearTrajectoryError(m));
503 
504 
505  completeSeed(*myPointer, newfts, thePropagatorAlongMomentum_, ilayer+1);
506  completeSeed(*myPointer, newfts, thePropagatorAlongMomentum_, ilayer+2);
507 
508 
509  } else {
510 
511 
512 
513  //Loop over compatible hits
514  int mea=0;
515  for(std::vector<TrajectoryMeasurement>::iterator tmItr = theFirstMeasurements_.begin(); tmItr !=theFirstMeasurements_.end(); ++tmItr) {
516 
517  mea++;
518 
519 
520  if (tmItr->recHit()->isValid() ) {
521  // Make a new helix as in fillClusterSeeds() but using the hit position
522  //std::cout << "InOutConversionSeedFinder::findSeeds hit R " << tmItr->recHit()->globalPosition().perp() << " Z " << tmItr->recHit()->globalPosition().z() << " " << tmItr->recHit()->globalPosition() << "\n";
523  GlobalPoint bcPos((theSecondBC_.position()).x(),(theSecondBC_.position()).y(),(theSecondBC_.position()).z());
524  GlobalVector dir = startingState.globalDirection();
525  GlobalPoint back1mm = tmItr->recHit()->globalPosition();
526 
527  back1mm -= dir.unit()*0.1;
528  //std::cout << " InOutConversionSeedFinder:::findSeeds going to make the helix using back1mm " << back1mm << "\n";
529  ConversionFastHelix helix(bcPos, tmItr->recHit()->globalPosition(), back1mm, &(*theMF_));
530 
531  helix.stateAtVertex();
532  //std::cout << " InOutConversionSeedFinder:::findSeeds helix status " <<helix.isValid() << std::endl;
533  if ( !helix.isValid() ) continue;
534 
535  track2InitialMomentum_= helix.stateAtVertex().momentum();
536 
537  //std::cout << "InOutConversionSeedFinder::findSeeds Updated estimatedPt = " << helix.stateAtVertex().momentum().perp() << " curvature " << helix.stateAtVertex().transverseCurvature() << "\n";
538  // << ", bcet = " << theBc->Et()
539  // << ", estimatedPt/bcet = " << estimatedPt/theBc->Et() << endl;
540 
541 
542  // Make a new FTS
543  FreeTrajectoryState newfts(GlobalTrajectoryParameters(
544  tmItr->recHit()->globalPosition(), startingState.globalDirection(),
545  helix.stateAtVertex().transverseCurvature(), 0, &(*theMF_)),
546  CurvilinearTrajectoryError(m));
547 
548  //std::cout << "InOutConversionSeedFinder::findSeeds new FTS charge " << newfts.charge() << "\n";
549 
550 
551  /*
552  // Soome diagnostic output
553  // may be useful - comparission of the basic cluster position
554  // with the ecal impact position of the track
555  TrajectoryStateOnSurface stateAtECAL
556  = forwardPropagator.propagate(newfts, ECALSurfaces::barrel());
557  if (!stateAtECAL.isValid() || abs(stateAtECAL.globalPosition().eta())>1.479) {
558  if (startingState.globalDirection().eta() > 0.) {
559  stateAtECAL = forwardPropagator.propagate(newfts,
560  ECALSurfaces::positiveEtaEndcap());
561  } else {
562  stateAtECAL = forwardPropagator.propagate(newfts,
563  ECALSurfaces::negativeEtaEndcap());
564  }
565  }
566  GlobalPoint ecalImpactPosition = stateAtECAL.isValid() ? stateAtECAL.globalPosition() : GlobalPoint(0.,0.,0.);
567  cout << "Projected fts positon at ECAL surface: " <<
568  ecalImpactPosition << " bc position: " << theBc->Position() << endl;
569  */
570 
571 
572  completeSeed(*tmItr, newfts, thePropagatorAlongMomentum_, ilayer+1);
573  completeSeed(*tmItr, newfts, thePropagatorAlongMomentum_, ilayer+2);
574 
575 
576  }
577 
578  }
579 
580  }
581 
582 
583  }
584 
585 
586 
587 }
588 
589 
590 
591 
592 
593 void InOutConversionSeedFinder::completeSeed(const TrajectoryMeasurement & m1,
594  FreeTrajectoryState & fts, const Propagator* propagator, int ilayer) const {
595 
596  //std::cout<< "InOutConversionSeedFinder::completeSeed ilayer " << ilayer << "\n";
597  // A seed is made from 2 Trajectory Measuremennts. The 1st is the input
598  // argument m1. This routine looks for the 2nd measurement in layer ilayer
599  // Begin by making a new much stricter MeasurementEstimator based on the
600  // position errors of the 1st hit.
601  printLayer(ilayer);
602 
603  MeasurementEstimator * newEstimator;
604  std::vector<const DetLayer*> allLayers=layerList();
605  const DetLayer * layer = allLayers[ilayer];
606 
608  // if ( layer->location() == GeomDetEnumerators::barrel ) {const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
609  // //std::cout << "InOutConversionSeedFinder::completeSeed **** Barrel on layer " << ilayer << " R= " << barrelLayer->specificSurface().radius() << "\n";
610  // } else {
611  // const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
612  // //std::cout << "InOutConversionSeedFinder::completeSeed **** Forw on layer " << ilayer << " Z= " << forwardLayer->specificSurface().position().z() << "\n";
615 
616 
617 
618 
619  if (layer->location() == GeomDetEnumerators::barrel ) {
620 
621  float dz = sqrt(the2ndHitdznSigma_*the2ndHitdznSigma_*m1.recHit()->globalPositionError().czz()
622  + the2ndHitdzConst_*the2ndHitdzConst_);
623  newEstimator = new ConversionBarrelEstimator(-the2ndHitdphi_, the2ndHitdphi_, -dz, dz);
624 
625  }
626  else {
627  float m1dr = sqrt(m1.recHit()->localPositionError().yy());
628  float dr = sqrt(the2ndHitdznSigma_*the2ndHitdznSigma_*m1dr*m1dr
629  + the2ndHitdzConst_*the2ndHitdznSigma_);
630 
631  newEstimator = new ConversionForwardEstimator(-the2ndHitdphi_, the2ndHitdphi_, dr);
632  }
633 
634 
635  //std::cout << "InOutConversionSeedFinder::completeSeed fts For the TSOS " << fts << "\n";
636 
637  TSOS tsos(fts, layer->surface() );
638 
639  if ( !tsos.isValid() ) {
640  //std::cout << "InOutConversionSeedFinder::completeSeed TSOS is not valid " << "\n";
641  }
642 
643  //std::cout << "InOutConversionSeedFinder::completeSeed TSOS " << tsos << "\n";
644  //std::cout << "InOutConversionSeedFinder::completeSeed propagationDirection " << int(propagator->propagationDirection() ) << "\n";
645  //std::cout << "InOutConversionSeedFinder::completeSeed pointer to estimator " << newEstimator << "\n";
646 
647  LayerMeasurements theLayerMeasurements_( *this->getMeasurementTracker(), *theTrackerData_ );
648  std::vector<TrajectoryMeasurement> measurements = theLayerMeasurements_.measurements( *layer, tsos, *propagator, *newEstimator);
649  //std::cout << "InOutConversionSeedFinder::completeSeed Found " << measurements.size() << " second hits " << "\n";
650  delete newEstimator;
651 
652  for(unsigned int i = 0; i < measurements.size(); ++i) {
653  if( measurements[i].recHit()->isValid() ) {
654  createSeed(m1, measurements[i]);
655  }
656  }
657 
658 
659 
660 
661 
662 
663 }
664 
665 
666 
667 void InOutConversionSeedFinder::createSeed(const TrajectoryMeasurement & m1, const TrajectoryMeasurement & m2) const {
668 
669  //std::cout << "InOutConversionSeedFinder::createSeed " << "\n";
670 
671  if ( m1.predictedState().isValid() ) {
672  GlobalTrajectoryParameters newgtp( m1.recHit()->globalPosition(), track2InitialMomentum_, track2Charge_, &(*theMF_) );
673  CurvilinearTrajectoryError errors = m1.predictedState().curvilinearError();
674  FreeTrajectoryState fts(newgtp, errors);
675 
676  TrajectoryStateOnSurface state1 = thePropagatorAlongMomentum_->propagate(fts, m1.recHit()->det()->surface());
677 
678  /*
679  //std::cout << "hit surface " << m1.recHit()->det()->surface().position() << "\n";
680  //std::cout << "prop to " << typeid( m1.recHit()->det()->surface() ).name() <<"\n";
681  //std::cout << "prop to first hit " << state1 << "\n";
682  //std::cout << "update to " << m1.recHit()->globalPosition() << "\n";
683  */
684 
685 
686 
687 
688  if ( state1.isValid() ) {
689 
690  TrajectoryStateOnSurface updatedState1 = theUpdator_.update(state1, *m1.recHit() );
691 
692  if ( updatedState1.isValid() ) {
693 
694  TrajectoryStateOnSurface state2 = thePropagatorAlongMomentum_->propagate(*updatedState1.freeTrajectoryState(), m2.recHit()->det()->surface());
695 
696  if ( state2.isValid() ) {
697 
698  TrajectoryStateOnSurface updatedState2 = theUpdator_.update(state2, *m2.recHit() );
699  TrajectoryMeasurement meas1(state1, updatedState1, m1.recHit() , m1.estimate(), m1.layer());
700  TrajectoryMeasurement meas2(state2, updatedState2, m2.recHit() , m2.estimate(), m2.layer());
701 
702  edm::OwnVector<TrackingRecHit> myHits;
703  myHits.push_back(meas1.recHit()->hit()->clone());
704  myHits.push_back(meas2.recHit()->hit()->clone());
705 
706  //std::cout << "InOutConversionSeedFinder::createSeed new seed " << "\n";
707  if ( nSeedsPerInputTrack_ >= maxNumberOfInOutSeedsPerInputTrack_ ) return;
708 
709 
710 
711  PTrajectoryStateOnDet const & ptsod= trajectoryStateTransform::persistentState(state2, meas2.recHit()->hit()->geographicalId().rawId() );
712  //std::cout << " InOutConversionSeedFinder::createSeed New seed parameters " << state2 << "\n";
713 
714 
715 
716  theSeeds_.push_back(TrajectorySeed( ptsod, myHits, alongMomentum ));
717  nSeedsPerInputTrack_++;
718 
719  //std::cout << "InOutConversionSeedFinder::createSeed New seed hit 1 R " << m1.recHit()->globalPosition().perp() << "\n";
720  //std::cout << "InOutConversionSeedFinder::createSeed New seed hit 2 R " << m2.recHit()->globalPosition().perp() << "\n";
721 
722  }
723  }
724  }
725  }
726 
727 
728 }
ConversionSeedFinder::thePropagatorAlongMomentum_
const Propagator * thePropagatorAlongMomentum_
Definition: ConversionSeedFinder.h:128
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
reco::CaloCluster::position
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:124
LayerMeasurements::measurements
std::vector< TrajectoryMeasurement > measurements(const DetLayer &layer, const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est) const
Definition: LayerMeasurements.cc:110
ConversionSeedFinder::layerList
std::vector< const DetLayer * > const & layerList() const
Definition: ConversionSeedFinder.h:69
GeometricSearchDet::surface
virtual const BoundSurface & surface() const =0
The surface of the GeometricSearchDet.
Propagator::propagate
virtual FreeTrajectoryState propagate(const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest) const final
Definition: Propagator.h:119
i
int i
Definition: DBlmapReader.cc:9
ConversionSeedFinder::getMeasurementTracker
const MeasurementTracker * getMeasurementTracker() const
Definition: ConversionSeedFinder.h:73
TrajectoryMeasurement::predictedState
TrajectoryStateOnSurface const & predictedState() const
Definition: TrajectoryMeasurement.h:174
TrajectoryMeasurement::recHit
ConstRecHitPointer const & recHit() const
Definition: TrajectoryMeasurement.h:209
benchmark_cfg.errors
tuple errors
Definition: benchmark_cfg.py:136
InOutConversionSeedFinder::findSeeds
virtual void findSeeds(const TrajectoryStateOnSurface &startingState, float signedpt, unsigned int startingLayer) const
Definition: InOutConversionSeedFinder.cc:410
DetLayer::location
virtual Location location() const =0
Which part of the detector (barrel, endcap)
AlgebraicMatrixID
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
Definition: AlgebraicROOTObjects.h:80
TrajectoryStateOnSurface::curvilinearError
const CurvilinearTrajectoryError & curvilinearError() const
Definition: TrajectoryStateOnSurface.h:110
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
trajectoryStateTransform::persistentState
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
Definition: TrajectoryStateTransform.cc:16
TrajectoryStateOnSurface::globalPosition
GlobalPoint globalPosition() const
Definition: TrajectoryStateOnSurface.h:89
AlgebraicSymMatrix55
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
Definition: AlgebraicROOTObjects.h:31
FreeTrajectoryState::charge
TrackCharge charge() const
Definition: FreeTrajectoryState.h:90
InOutConversionSeedFinder::theFirstMeasurements_
std::vector< TrajectoryMeasurement > theFirstMeasurements_
Definition: InOutConversionSeedFinder.h:90
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
ConversionSeedFinder::theSeeds_
TrajectorySeedCollection theSeeds_
Definition: ConversionSeedFinder.h:100
detailsBasic3DVector::z
float float float z
Definition: extBasic3DVector.h:15
ConversionSeedFinder::thePropagatorOppositeToMomentum_
const Propagator * thePropagatorOppositeToMomentum_
Definition: ConversionSeedFinder.h:129
InOutConversionSeedFinder::inputTracks_
std::vector< Trajectory > inputTracks_
Definition: InOutConversionSeedFinder.h:88
InOutConversionSeedFinder::myPointer
TrajectoryMeasurement * myPointer
Definition: InOutConversionSeedFinder.h:86
InOutConversionSeedFinder::getSecondCaloClusters
std::vector< const reco::CaloCluster * > getSecondCaloClusters(const GlobalPoint &conversionPosition, float charge) const
Definition: InOutConversionSeedFinder.cc:367
InOutConversionSeedFinder::startSeed
void startSeed(const FreeTrajectoryState *fts, const TrajectoryStateOnSurface &stateAtPreviousLayer, int charge, int layer) const
Definition: InOutConversionSeedFinder.cc:313
edm::OwnVector::push_back
void push_back(D *&d)
Definition: OwnVector.h:274
KFUpdator::update
TrajectoryStateOnSurface update(const TrajectoryStateOnSurface &, const TrackingRecHit &) const
Definition: KFUpdator.cc:10
TrajectoryStateOnSurface::freeTrajectoryState
FreeTrajectoryState const * freeTrajectoryState(bool withErrors=true) const
Definition: TrajectoryStateOnSurface.h:80
ConversionSeedFinder::printLayer
void printLayer(int i) const
Definition: ConversionSeedFinder.cc:151
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
query.result
tuple result
Definition: query.py:137
TrajectoryMeasurement::layer
const DetLayer * layer() const
Definition: TrajectoryMeasurement.h:215
ConversionSeedFinder::theTrackerGeom_
const TrackingGeometry * theTrackerGeom_
Definition: ConversionSeedFinder.h:106
dbtoconf.conf
tuple conf
Definition: dbtoconf.py:185
InOutConversionSeedFinder::bcCollection_
edm::Handle< edm::View< reco::CaloCluster > > bcCollection_
Definition: InOutConversionSeedFinder.h:93
InOutConversionSeedFinder::completeSeed
void completeSeed(const TrajectoryMeasurement &m1, FreeTrajectoryState &fts, const Propagator *propagator, int ilayer) const
Definition: InOutConversionSeedFinder.cc:593
TrackingGeometry::idToDet
virtual const GeomDet * idToDet(DetId) const =0
InOutConversionSeedFinder::InOutConversionSeedFinder
InOutConversionSeedFinder(const edm::ParameterSet &config, edm::ConsumesCollector &&iC)
Definition: InOutConversionSeedFinder.cc:21
TrajectoryMeasurement::updatedState
TrajectoryStateOnSurface const & updatedState() const
Definition: TrajectoryMeasurement.h:190
dir
dbl *** dir
Definition: mlp_gen.cc:35
x
Definition: DDAxes.h:10
ConversionSeedFinder::theTrackerData_
edm::Handle< MeasurementTrackerEvent > theTrackerData_
Definition: ConversionSeedFinder.h:132
InOutConversionSeedFinder::createSeed
void createSeed(const TrajectoryMeasurement &m1, const TrajectoryMeasurement &m2) const
Definition: InOutConversionSeedFinder.cc:667
InOutConversionSeedFinder::theOutInTracks_
std::vector< Trajectory > theOutInTracks_
Definition: InOutConversionSeedFinder.h:89
InOutConversionSeedFinder::makeSeeds
virtual void makeSeeds(const edm::Handle< edm::View< reco::CaloCluster > > &allBc) const
Definition: InOutConversionSeedFinder.cc:45
InOutConversionSeedFinder::fillClusterSeeds
virtual void fillClusterSeeds() const
Definition: InOutConversionSeedFinder.cc:68
TrajectoryStateOnSurface::globalDirection
GlobalVector globalDirection() const
Definition: TrajectoryStateOnSurface.h:95
Generated for CMSSW Reference Manual by doxygen 1.8.5