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