CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
KalmanAlignmentAlgorithm.cc
Go to the documentation of this file.
1 
4 
5 // includes for alignment
10 
13 
19 
21 
24 
32 
34 
37 
38 // miscellaneous includes
40 #include "CLHEP/Random/RandGauss.h"
41 #include <fstream>
42 
43 using namespace std;
44 
45 
47  AlignmentAlgorithmBase( config ),
48  theConfiguration( config )
49 {}
50 
51 
53 
54 
58  AlignableExtras* extras,
60 {
62 
63  theMergerFlag = theConfiguration.getParameter<bool>( "MergeResults" );
64 
65  if ( theMergerFlag )
66  {
67  mergeResults();
68  } else {
69  theParameterStore = store;
70  theNavigator = new AlignableNavigator( tracker->components() );
71  theSelector = new AlignmentParameterSelector( tracker );
72 
74 
77 
79  }
80 }
81 
82 
84 {
85  if ( theMergerFlag )
86  {
87  return; // only merging mode. nothing to do here.
88  }
89 
90  cout << "[KalmanAlignmentAlgorithm::terminate] start ..." << endl;
91 
92  set< Alignable* > allAlignables;
93  vector< Alignable* > alignablesToWrite;
94 
95  AlignmentSetupCollection::const_iterator itSetup;
96  for ( itSetup = theAlignmentSetups.begin(); itSetup != theAlignmentSetups.end(); ++itSetup )
97  {
98  delete (*itSetup)->alignmentUpdator();
99 
100  const vector< Alignable* >& alignablesFromMetrics = (*itSetup)->metricsUpdator()->alignables();
101  cout << "[KalmanAlignmentAlgorithm::terminate] The metrics updator for setup \'" << (*itSetup)->id()
102  << "\' holds " << alignablesFromMetrics.size() << " alignables" << endl;
103  allAlignables.insert( alignablesFromMetrics.begin(), alignablesFromMetrics.end() );
104  }
105 
106  for ( set< Alignable* >::iterator it = allAlignables.begin(); it != allAlignables.end(); ++it )
107  {
108  AlignmentParameters* alignmentParameters = ( *it )->alignmentParameters();
109 
110  if ( alignmentParameters != 0 )
111  {
112  KalmanAlignmentUserVariables* userVariables =
113  dynamic_cast< KalmanAlignmentUserVariables* >( alignmentParameters->userVariables() );
114 
115  if ( userVariables != 0 && userVariables->numberOfUpdates() > 0 )
116  {
117  userVariables->update( true );
118  userVariables->histogramParameters( "KalmanAlignmentAlgorithm" );
119  alignablesToWrite.push_back( *it );
120  }
121  }
122  }
123 
124  if ( theConfiguration.getUntrackedParameter< bool >( "WriteAlignmentParameters", false ) )
125  {
126  AlignmentIORoot alignmentIO;
127  int ierr = 0;
128  string output = theConfiguration.getParameter< string >( "OutputFile" );
129 
130  cout << "Write data for " << alignablesToWrite.size() << " alignables ..." << endl;
131 
132  // Write output to "iteration 1", ...
133  alignmentIO.writeAlignmentParameters( alignablesToWrite, output.c_str(), 1, false, ierr );
134  // ... or, if "iteration 1" already exists, write it to "highest iteration + 1"
135  if ( ierr == -1 ) alignmentIO.writeAlignmentParameters( alignablesToWrite, output.c_str(), -1, false, ierr );
136  }
137 
139 
140  string timingLogFile = theConfiguration.getUntrackedParameter< string >( "TimingLogFile", "timing.log" );
141  cout << "The usage of TimingReport from Utilities/Timing! Timing not written to log file" << endl;
142 
143  delete theNavigator;
144 
145  cout << "[KalmanAlignmentAlgorithm::terminate] ... done." << endl;
146 }
147 
148 
150 {
151  if ( theMergerFlag ) return; // only merging mode. nothing to do here.
152 
153  static int iEvent = 1;
154  if ( iEvent % 100 == 0 ) cout << "[KalmanAlignmentAlgorithm::run] Event Nr. " << iEvent << endl;
155  iEvent++;
156 
157  edm::ESHandle< MagneticField > aMagneticField;
158  setup.get< IdealMagneticFieldRecord >().get( aMagneticField );
159 
160  try
161  {
162  // Run the refitter algorithm
164  const reco::BeamSpot &beamSpot = eventInfo.beamSpot();
165  TrackletCollection refittedTracklets = theRefitter->refitTracks( setup, theAlignmentSetups, tracks, &beamSpot );
166 
167  // Associate tracklets to alignment setups
168  map< AlignmentSetup*, TrackletCollection > setupToTrackletMap;
169  TrackletCollection::iterator itTracklet;
170  for ( itTracklet = refittedTracklets.begin(); itTracklet != refittedTracklets.end(); ++itTracklet )
171  setupToTrackletMap[(*itTracklet)->alignmentSetup()].push_back( *itTracklet );
172 
173  // Iterate on alignment setups
174  map< AlignmentSetup*, TrackletCollection >::iterator itMap;
175  for ( itMap = setupToTrackletMap.begin(); itMap != setupToTrackletMap.end(); ++itMap )
176  {
179 
180  TrackletCollection::iterator itTracklet;
181  for ( itTracklet = itMap->second.begin(); itTracklet != itMap->second.end(); ++itTracklet )
182  {
183  tracklets.push_back( (*itTracklet)->trajTrackPair() );
184  external.push_back( (*itTracklet)->externalPrediction() );
185  }
186 
187  // Construct reference trajectories
188  ReferenceTrajectoryCollection trajectories =
189  itMap->first->trajectoryFactory()->trajectories( setup, tracklets, external, eventInfo.beamSpot() );
190 
191  ReferenceTrajectoryCollection::iterator itTrajectories;
192 
193  // Run the alignment algorithm.
194  for ( itTrajectories = trajectories.begin(); itTrajectories != trajectories.end(); ++itTrajectories )
195  {
196  itMap->first->alignmentUpdator()->process( *itTrajectories, theParameterStore, theNavigator,
197  itMap->first->metricsUpdator(), aMagneticField.product() );
198 
199  KalmanAlignmentDataCollector::fillHistogram( "Trajectory_RecHits", (*itTrajectories)->recHits().size() );
200  }
201  }
202  }
203  catch( cms::Exception& exception )
204  {
205  cout << exception.what() << endl;
206  terminate(setup);
207  throw exception;
208  }
209 }
210 
211 
213 {
214  //Retrieve tracker topology from geometry
215  edm::ESHandle<TrackerTopology> tTopoHandle;
216  setup.get<IdealGeometryRecord>().get(tTopoHandle);
217  const TrackerTopology* const tTopo = tTopoHandle.product();
218 
219  // Just to be sure, set all APEs to zero ...
220  setAPEToZero();
221 
222  const edm::ParameterSet initConfig = theConfiguration.getParameter< edm::ParameterSet >( "ParameterConfig" );
223 
224  int updateGraph = initConfig.getUntrackedParameter< int >( "UpdateGraphs", 100 );
225 
226  bool addPositionError = false;// = initConfig.getUntrackedParameter< bool >( "AddPositionError", true );
227 
228  int seed = initConfig.getUntrackedParameter< int >( "RandomSeed", 1726354 );
229  CLHEP::HepRandom::createInstance();
230  CLHEP::HepRandom::setTheSeed( seed );
231 
232  bool applyXShifts = initConfig.getUntrackedParameter< bool >( "ApplyXShifts", false );
233  bool applyYShifts = initConfig.getUntrackedParameter< bool >( "ApplyYShifts", false );
234  bool applyZShifts = initConfig.getUntrackedParameter< bool >( "ApplyZShifts", false );
235  bool applyXRots = initConfig.getUntrackedParameter< bool >( "ApplyXRotations", false );
236  bool applyYRots = initConfig.getUntrackedParameter< bool >( "ApplyYRotations", false );
237  bool applyZRots = initConfig.getUntrackedParameter< bool >( "ApplyZRotations", false );
238 
239  bool applyRandomStartValues = initConfig.getUntrackedParameter< bool >( "ApplyRandomStartValues", false );
240  if ( applyRandomStartValues )
241  cout << "[KalmanAlignmentAlgorithm::initializeAlignmentParameters] ADDING RANDOM START VALUES!!!" << endl;
242 
243  bool applyCurl = initConfig.getUntrackedParameter< bool >( "ApplyCurl", false );
244  double curlConst = initConfig.getUntrackedParameter< double >( "CurlConstant", 1e-6 );
245 
246  bool applyShifts = applyXShifts || applyYShifts || applyZShifts;
247  bool applyRots = applyXRots || applyYRots || applyZRots;
248  //bool applyMisalignment = applyShifts || applyRots || applyCurl;
249 
250  AlgebraicVector displacement( 3 );
251  AlgebraicVector eulerAngles( 3 );
252 
253  AlgebraicVector startParameters( 6, 0 );
254  AlgebraicSymMatrix startError( 6, 0 );
255 
256  AlgebraicVector randSig( 6, 0 );
257 
258  vector< string > initSelection = initConfig.getParameter< vector<string> >( "InitializationSelector" );
259 
260  vector< string >::iterator itInitSel;
261  for ( itInitSel = initSelection.begin(); itInitSel != initSelection.end(); ++itInitSel )
262  {
263  const edm::ParameterSet config = initConfig.getParameter< edm::ParameterSet >( *itInitSel );
264 
265  addPositionError = initConfig.getUntrackedParameter< bool >( "AddPositionError", false );
266 
267  double xAPEError = initConfig.getUntrackedParameter< double >( "SigmaXPositionError", 2e-2 );
268  double yAPEError = initConfig.getUntrackedParameter< double >( "SigmaYPositionError", 2e-2 );
269  double zAPEError = initConfig.getUntrackedParameter< double >( "SigmaZPositionError", 2e-2 );
270 
271  double sigmaXShift = config.getUntrackedParameter< double >( "SigmaXShifts", 4e-2 );
272  double sigmaYShift = config.getUntrackedParameter< double >( "SigmaYShifts", 4e-2 );
273  double sigmaZShift = config.getUntrackedParameter< double >( "SigmaZShifts", 4e-2 );
274  double sigmaXRot = config.getUntrackedParameter< double >( "SigmaXRotations", 5e-4 );
275  double sigmaYRot = config.getUntrackedParameter< double >( "SigmaYRotations", 5e-4 );
276  double sigmaZRot = config.getUntrackedParameter< double >( "SigmaZRotations", 5e-4 );
277 
278  randSig[0] = sigmaXShift; randSig[1] = sigmaYShift; randSig[2] = sigmaZShift;
279  randSig[3] = sigmaXRot; randSig[4] = sigmaYRot; randSig[5] = sigmaZRot;
280 
281  startError[0][0] = config.getUntrackedParameter< double >( "XShiftsStartError", 4e-4 );
282  startError[1][1] = config.getUntrackedParameter< double >( "YShiftsStartError", 4e-4 );
283  startError[2][2] = config.getUntrackedParameter< double >( "ZShiftsStartError", 4e-4 );
284  startError[3][3] = config.getUntrackedParameter< double >( "XRotationsStartError", 3e-5 );
285  startError[4][4] = config.getUntrackedParameter< double >( "YRotationsStartError", 3e-5 );
286  startError[5][5] = config.getUntrackedParameter< double >( "ZRotationsStartError", 3e-5 );
287 
288  const vector< char > dummyParamSelector( 6, '0' );
289  const vector< string > alignableSelector = config.getParameter< vector<string> >( "AlignableSelection" );
290 
291  vector< string >::const_iterator itAliSel;
292  for ( itAliSel = alignableSelector.begin(); itAliSel != alignableSelector.end(); ++itAliSel )
293  {
294  theSelector->addSelection( *itAliSel, dummyParamSelector );
295  cout << "[" << *itInitSel << "] add selection: " << *itAliSel << endl;
296  }
297 
298  vector< Alignable* >::iterator itAlignable;
299  vector< Alignable* > alignables;
300  vector< Alignable* > alignablesFromSelector = theSelector->selectedAlignables();
301  for ( itAlignable = alignablesFromSelector.begin(); itAlignable != alignablesFromSelector.end(); ++itAlignable )
302  //if ( (*itAlignable)->alignmentParameters() )
303  alignables.push_back( *itAlignable );
304 
305  cout << "[" << *itInitSel << "] total number of selected alignables = " << alignables.size() << endl;
306 
307  sort( alignables.begin(), alignables.end(), *this );
308 
309  // Apply existing alignment parameters.
310  map< Alignable*, vector< AlignmentParameters* > > alignmentParametersMap;
311 
312  int iApply = 0;
313  bool readParam = config.getUntrackedParameter< bool >( "ReadParametersFromFile", false );
314  bool readCovar = config.getUntrackedParameter< bool >( "ReadCovarianceFromFile", false );
315  bool applyParam = config.getUntrackedParameter< bool >( "ApplyParametersFromFile", false );
316  bool applyCovar = config.getUntrackedParameter< bool >( "ApplyErrorFromFile", false );
317  if ( readParam || readCovar || applyParam || applyCovar )
318  {
319  string file = config.getUntrackedParameter< string >( "FileName", "Input.root" );
320  int ierr = 0;
321  int iter = 1;
322 
323  AlignmentIORoot alignmentIO;
324  while ( !ierr )
325  {
326  cout << "[" << *itInitSel << "] read alignment parameters. file / iteration = " << file << " / " << iter << endl;
327  vector< AlignmentParameters* > alignmentParameters = alignmentIO.readAlignmentParameters( alignables, file.c_str(), iter, ierr );
328  cout << "[" << *itInitSel << "] #param / ierr = " << alignmentParameters.size() << " / " << ierr << endl;
329 
330  vector< AlignmentParameters* >::iterator itParam;
331  for ( itParam = alignmentParameters.begin(); itParam != alignmentParameters.end(); ++itParam )
332  alignmentParametersMap[(*itParam)->alignable()].push_back( *itParam );
333 
334  ++iter;
335  }
336  }
337 
338  int iAlign = 0;
339 
340  for ( itAlignable = alignables.begin(); itAlignable != alignables.end(); itAlignable++ )
341  {
342  if ( (*itAlignable)->alignmentParameters() == 0 )
343  {
344  cout << "[KalmanAlignmentAlgorithm::initializeAlignmentParameters] "
345  << "alignable is not associated with alignment parameters --> skip" << endl;
346  continue;
347  }
348 
349  if ( (*itAlignable)->alignmentParameters()->type() != AlignmentParametersFactory::kRigidBody )
350  {
351  cout << "[KalmanAlignmentAlgorithm::initializeAlignmentParameters] "
352  << "alignable is not associated with rigid body alignment parameters --> skip" << endl;
353  continue;
354  }
355 
356  displacement[0] = applyXShifts ? sigmaXShift*CLHEP::RandGauss::shoot() : 0.;
357  displacement[1] = applyYShifts ? sigmaZShift*CLHEP::RandGauss::shoot() : 0.;
358  displacement[2] = applyZShifts ? sigmaYShift*CLHEP::RandGauss::shoot() : 0.;
359 
360  if ( applyShifts )
361  {
362  align::LocalVector localShift( displacement[0], displacement[1], displacement[2] );
363  align::GlobalVector globalShift = ( *itAlignable )->surface().toGlobal( localShift );
364  ( *itAlignable )->move( globalShift );
365  }
366 
367  align::EulerAngles eulerAngles( 3 );
368 
369  eulerAngles[0] = applyXRots ? sigmaXRot*CLHEP::RandGauss::shoot() : 0.;
370  eulerAngles[1] = applyYRots ? sigmaYRot*CLHEP::RandGauss::shoot() : 0.;
371  eulerAngles[2] = applyZRots ? sigmaZRot*CLHEP::RandGauss::shoot() : 0.;
372 
373  if ( applyRots )
374  {
375  align::RotationType localRotation = align::toMatrix( eulerAngles );
376  ( *itAlignable )->rotateInLocalFrame( localRotation );
377  }
378 
379  if ( applyCurl )
380  {
381  double radius = ( *itAlignable )->globalPosition().perp();
382  ( *itAlignable )->rotateAroundGlobalZ( curlConst*radius );
383  }
384 
385  if ( addPositionError )
386  {
387  LocalVector localError( sqrt(xAPEError), sqrt(yAPEError), sqrt(zAPEError) );
388  GlobalVector globalError = (*itAlignable)->surface().toGlobal( localError );
389  AlignmentPositionError ape( globalError.x(), globalError.y(), globalError.z() );
390  // FIXME (GF): The above does not seem to be a correct error transformation!
391  ( *itAlignable )->setAlignmentPositionError( ape, true ); // true: propagate down
392  }
393 
394  //AlgebraicVector trueParameters( 6 );
395  //trueParameters[0] = -displacement[0];
396  //trueParameters[1] = -displacement[1];
397  //trueParameters[2] = -displacement[2];
398  //trueParameters[3] = -eulerAngles[0];
399  //trueParameters[4] = -eulerAngles[1];
400  //trueParameters[5] = -eulerAngles[2];
401 
402  if ( (*itAlignable)->alignmentParameters() != 0 )
403  {
404  AlignmentParameters* alignmentParameters;
405  if ( readParam && readCovar )
406  {
407  if ( alignmentParametersMap.find( *itAlignable ) == alignmentParametersMap.end() )
408  {
409  //cout << "apply param and cov from FILE -> none stored, apply DEFAULT " << endl;
410  alignmentParameters = (*itAlignable)->alignmentParameters()->clone( startParameters, startError );
411  alignmentParameters->setUserVariables( new KalmanAlignmentUserVariables( *itAlignable, tTopo, updateGraph ) );
412  }
413  else
414  {
415  //cout << "apply param and cov from FILE" << endl;
416  alignmentParameters = alignmentParametersMap[*itAlignable].back();
417  KalmanAlignmentUserVariables* userVariables = new KalmanAlignmentUserVariables( *itAlignable, tTopo, updateGraph );
418  userVariables->update( alignmentParameters );
419  alignmentParameters->setUserVariables( userVariables );
420  }
421  }
422  else if ( readParam )
423  {
424  if ( alignmentParametersMap.find( *itAlignable ) == alignmentParametersMap.end() )
425  {
426  alignmentParameters = (*itAlignable)->alignmentParameters()->clone( startParameters, startError );
427  alignmentParameters->setUserVariables( new KalmanAlignmentUserVariables( *itAlignable, tTopo, updateGraph ) );
428  }
429  else
430  {
431  AlgebraicVector parameters = alignmentParametersMap[*itAlignable].back()->parameters();
432  alignmentParameters = (*itAlignable)->alignmentParameters()->clone( parameters, startError );
433  KalmanAlignmentUserVariables* userVariables = new KalmanAlignmentUserVariables( *itAlignable, tTopo, updateGraph );
434  userVariables->update( alignmentParameters );
435  alignmentParameters->setUserVariables( userVariables );
436  }
437  }
438  else
439  {
440  //cout << "apply DEFAULT param and cov" << endl;
441  alignmentParameters = (*itAlignable)->alignmentParameters()->clone( startParameters, startError );
442  //alignmentParameters = (*itAlignable)->alignmentParameters()->clone( trueParameters, startError );
443  alignmentParameters->setUserVariables( new KalmanAlignmentUserVariables( *itAlignable, tTopo, updateGraph ) );
444  }
445 
446  (*itAlignable)->setAlignmentParameters( alignmentParameters );
447  //if ( applyParam ) theParameterStore->applyParameters( *itAlignable );
448 
449  if ( applyRandomStartValues )
450  {
451  cout << "applying random start values" << endl;
452 
453  AlgebraicVector randomStartParameters = alignmentParameters->parameters();
454  AlgebraicSymMatrix randomStartErrors = alignmentParameters->covariance();
455 
456  for ( int iParam = 0; iParam < randomStartParameters.num_row(); ++iParam )
457  {
458  randomStartParameters[iParam] += sqrt(randSig[iParam])*CLHEP::RandGauss::shoot();
459  //randomStartErrors[iParam][iParam] += randSig[iParam]*randSig[iParam];
460  }
461 
462  //cout << randomStartParameters << endl;
463 
464  alignmentParameters = (*itAlignable)->alignmentParameters()->clone( randomStartParameters, randomStartErrors );
465  (*itAlignable)->setAlignmentParameters( alignmentParameters );
466  }
467 
468  (*itAlignable)->alignmentParameters()->setValid( true );
469 
470  }
471 
472  if ( ( applyParam || applyCovar ) && alignmentParametersMap.find( *itAlignable ) != alignmentParametersMap.end() )
473  {
474  ++iApply;
475 
476  vector< AlignmentParameters* > allAlignmentParameters = alignmentParametersMap[*itAlignable];
477  vector< AlignmentParameters* >::iterator itParam;
478 
479  for ( itParam = allAlignmentParameters.begin(); itParam != allAlignmentParameters.end(); ++itParam )
480  applyAlignmentParameters( *itAlignable, *itParam, applyParam, applyCovar );
481 
482  if ( ( *itAlignable )->alignmentParameters() )
483  {
484  KalmanAlignmentUserVariables* userVariables =
485  dynamic_cast< KalmanAlignmentUserVariables* >( ( *itAlignable )->alignmentParameters()->userVariables() );
486  if ( userVariables ) { ++iAlign; userVariables->setAlignmentFlag( true ); }
487  }
488  }
489  }
490 
491  cout << "[" << *itInitSel << "] Set the alignment flag for " << iAlign << " alignables." << endl;
492  cout << "[" << *itInitSel << "] number of applied parameters: " << iApply << endl;
493  theSelector->clear();
494  }
495 
496 }
497 
498 
500 {
501  // Read the setups from the config-file. They define which rec-hits are refitted to tracklets and whether they are used
502  // as an external measurement or not.
503 
504  const edm::ParameterSet initConfig = theConfiguration.getParameter< edm::ParameterSet >( "AlgorithmConfig" );
505  const vector<string> selSetup = initConfig.getParameter< vector<string> >( "Setups" );
506 
507  for ( vector<string>::const_iterator itSel = selSetup.begin(); itSel != selSetup.end(); ++itSel )
508  {
509  cout << "[KalmanAlignmentAlgorithm::initializeAlignmentSetups] Add AlignmentSetup: " << *itSel << endl;
510 
511  const edm::ParameterSet confSetup = initConfig.getParameter< edm::ParameterSet >( *itSel );
512 
513  string strPropDir = confSetup.getUntrackedParameter< string >( "PropagationDirection", "alongMomentum" );
514  string strSortingDir = confSetup.getUntrackedParameter< string >( "SortingDirection", "SortInsideOut" );
515  vector<int> trackingIDs = confSetup.getParameter< vector<int> >( "Tracking" );
516  unsigned int minTrackingHits = confSetup.getUntrackedParameter< unsigned int >( "MinTrackingHits", 0 );
517 
518  string strExternalPropDir = confSetup.getUntrackedParameter< string >( "ExternalPropagationDirection", "alongMomentum" );
519  string strExternalSortingDir = confSetup.getUntrackedParameter< string >( "ExternalSortingDirection", "SortInsideOut" );
520  vector<int> externalIDs = confSetup.getParameter< vector<int> >( "External" );
521  unsigned int minExternalHits = confSetup.getUntrackedParameter< unsigned int >( "MinExternalHits", 0 );
522 
523  edm::ESHandle< TrajectoryFitter > aTrajectoryFitter;
524  string fitterName = confSetup.getUntrackedParameter< string >( "Fitter", "KFFittingSmoother" );
525  setup.get<TrajectoryFitter::Record>().get( fitterName, aTrajectoryFitter );
526 
527  double outlierEstimateCut = 5.;
528 
529  const KFFittingSmoother* aFittingSmoother = dynamic_cast< const KFFittingSmoother* >( aTrajectoryFitter.product() );
530  if ( aFittingSmoother )
531  {
532  KFTrajectoryFitter* fitter = 0;
533  KFTrajectorySmoother* smoother = 0;
535 
536  KFTrajectoryFitter* externalFitter = 0;
537  KFTrajectorySmoother* externalSmoother = 0;
538 
539  PropagationDirection fitterDir = getDirection( strPropDir );
540  PropagationDirection externalFitterDir = getDirection( strExternalPropDir );
541 
542  PropagationDirection smootherDir = oppositeDirection( fitterDir );
543  PropagationDirection externalSmootherDir = oppositeDirection( externalFitterDir );
544 
545  const KFTrajectoryFitter* aKFFitter = dynamic_cast< const KFTrajectoryFitter* >( aFittingSmoother->fitter() );
546  if ( aKFFitter )
547  {
548  PropagatorWithMaterial propagator( fitterDir, 0.106, aKFFitter->propagator()->magneticField() );
549  Chi2MeasurementEstimator estimator( 30. );
550  fitter = new KFTrajectoryFitter( &propagator, updator, &estimator );
551 // fitter = new KFTrajectoryFitter( &propagator, aKFFitter->updator(), &estimator );
552 
553  AnalyticalPropagator externalPropagator( aKFFitter->propagator()->magneticField(), externalFitterDir );
554  Chi2MeasurementEstimator externalEstimator( 1000. );
555  externalFitter = new KFTrajectoryFitter( &externalPropagator, aKFFitter->updator(), &externalEstimator );
556  }
557 
558  const KFTrajectorySmoother* aKFSmoother = dynamic_cast< const KFTrajectorySmoother* >( aFittingSmoother->smoother() );
559  if ( aKFSmoother )
560  {
561 
562  PropagatorWithMaterial propagator( smootherDir, 0.106, aKFSmoother->alongPropagator()->magneticField() );
563  Chi2MeasurementEstimator estimator( 30. );
564  smoother = new KFTrajectorySmoother( &propagator, updator, &estimator );
565 // smoother = new KFTrajectorySmoother( &propagator, aKFFitter->updator(), &estimator );
566 
567  AnalyticalPropagator externalPropagator( aKFSmoother->alongPropagator()->magneticField(), externalSmootherDir );
568  Chi2MeasurementEstimator externalEstimator( 1000. );
569  externalSmoother = new KFTrajectorySmoother( &externalPropagator, aKFSmoother->updator(), &externalEstimator );
570  }
571 
572  if ( fitter && smoother )
573  {
574  KFFittingSmoother* fittingSmoother = new KFFittingSmoother( *fitter, *smoother, outlierEstimateCut );
575  KFFittingSmoother* externalFittingSmoother = new KFFittingSmoother( *externalFitter, *externalSmoother );
576 // KFFittingSmoother* fittingSmoother = aFittingSmoother->clone();
577 // KFFittingSmoother* externalFittingSmoother = aFittingSmoother->clone();
578 
579  string identifier;
581 
582  config = confSetup.getParameter< edm::ParameterSet >( "TrajectoryFactory" );
583  identifier = config.getParameter< string >( "TrajectoryFactoryName" );
584  cout << "TrajectoryFactoryPlugin::get() ... " << identifier << endl;
585  TrajectoryFactoryBase* trajectoryFactory = TrajectoryFactoryPlugin::get()->create( identifier, config );
586 
587  config = confSetup.getParameter< edm::ParameterSet >( "AlignmentUpdator" );
588  identifier = config.getParameter< string >( "AlignmentUpdatorName" );
589  KalmanAlignmentUpdator* alignmentUpdator = KalmanAlignmentUpdatorPlugin::get()->create( identifier, config );
590 
591  config = confSetup.getParameter< edm::ParameterSet >( "MetricsUpdator" );
592  identifier = config.getParameter< string >( "MetricsUpdatorName" );
593  KalmanAlignmentMetricsUpdator* metricsUpdator = KalmanAlignmentMetricsUpdatorPlugin::get()->create( identifier, config );
594 
595  KalmanAlignmentSetup::SortingDirection sortingDir = getSortingDirection( strSortingDir );
596  KalmanAlignmentSetup::SortingDirection externalSortingDir = getSortingDirection( strExternalSortingDir );
597 
598  AlignmentSetup* anAlignmentSetup
599  = new AlignmentSetup( *itSel,
600  fittingSmoother, fitter->propagator(), trackingIDs, minTrackingHits, sortingDir,
601  externalFittingSmoother, externalFitter->propagator(), externalIDs, minExternalHits, externalSortingDir,
602  trajectoryFactory, alignmentUpdator, metricsUpdator );
603 
604  theAlignmentSetups.push_back( anAlignmentSetup );
605 
606  delete fittingSmoother;
607  delete fitter;
608  delete smoother;
609 
610  delete externalFittingSmoother;
611  delete externalFitter;
612  delete externalSmoother;
613  }
614  else
615  {
616  throw cms::Exception( "BadConfig" ) << "[KalmanAlignmentAlgorithm::initializeAlignmentSetups] "
617  << "Instance of class KFFittingSmoother has no KFTrajectoryFitter/KFTrajectorySmoother.";
618  }
619 
620  delete updator;
621  }
622  else
623  {
624  throw cms::Exception( "BadConfig" ) << "[KalmanAlignmentAlgorithm::initializeAlignmentSetups] "
625  << "No instance of class KFFittingSmoother registered to the TrackingComponentsRecord.";
626  }
627  }
628 
629  cout << "[KalmanAlignmentAlgorithm::initializeAlignmentSetups] I'm using " << theAlignmentSetups.size() << " AlignmentSetup(s)." << endl;
630 
631 }
632 
633 
636 {
637  if ( sortDir == "SortInsideOut" ) { return KalmanAlignmentSetup::sortInsideOut; }
638  if ( sortDir == "SortOutsideIn" ) { return KalmanAlignmentSetup::sortOutsideIn; }
639  if ( sortDir == "SortUpsideDown" ) { return KalmanAlignmentSetup::sortUpsideDown; }
640  if ( sortDir == "SortDownsideUp" ) { return KalmanAlignmentSetup::sortDownsideUp; }
641 
642  throw cms::Exception( "BadConfig" ) << "[KalmanAlignmentAlgorithm::getSortingDirection] "
643  << "Unknown sorting direction: " << sortDir << std::endl;
644 }
645 
646 
647 void
649  bool applyPar, bool applyCov ) const
650 {
651  RigidBodyAlignmentParameters* rbap = dynamic_cast<RigidBodyAlignmentParameters*>( par );
652 
653  if ( !rbap )
654  throw cms::Exception( "BadConfig" ) << "applyParameters: provided alignable does not have rigid body alignment parameters";
655 
656  if ( applyPar )
657  {
659  const AlignableSurface& alignableSurface = ali->surface();
660  ali->move( alignableSurface.toGlobal( align::LocalVector( shift[0], shift[1], shift[2] ) ) );
661 
662  align::EulerAngles angles = rbap->rotation();
663  if ( angles.normsq() > 1e-10 ) ali->rotateInLocalFrame( align::toMatrix( angles ) );
664  }
665 
666  if ( applyCov )
667  {
668  const AlgebraicSymMatrix& aliCov = rbap->covariance();
669  LocalVector localError( sqrt(aliCov[0][0]), sqrt(aliCov[1][1]), sqrt(aliCov[2][2]) );
670  GlobalVector globalError = ali->surface().toGlobal( localError );
671  AlignmentPositionError ape( globalError.x(), globalError.y(), globalError.z() );
672  // FIXME (GF): The above does not seem to be a correct error transformation!
673  ali->setAlignmentPositionError( ape, true ); // true: propagate down
674  }
675 }
676 
677 
679  vector<Alignable*>& comps ) const
680 {
681  comps.push_back( ali );
682  vector<Alignable*> nextComps = ali->components();
683  vector<Alignable*>::iterator it;
684  for ( it = nextComps.begin(); it != nextComps.end(); ++it ) getComponents( *it, comps );
685 }
686 
687 
689 {
690  std::cout << "[KalmanAlignmentAlgorithm::mergeResults] START MERGING RESULTS" << std::endl;
691 
693 
694  vector<string> inFileNames = mergeConf.getParameter< vector<string> >( "InputMergeFileNames" );
695  string outFileName = mergeConf.getParameter<string>( "OutputMergeFileName" );
696 
697  bool applyPar = mergeConf.getParameter<bool>( "ApplyParameters" );
698  bool applyCov = mergeConf.getParameter<bool>( "ApplyErrors" );
699 
700  map< Alignable*, vector< AlignmentParameters* > > alignmentParametersMap;
701 
702  vector< Alignable* > allAlignables;
703  getComponents( theTracker, allAlignables );
704 
705  cout << "allAlignables.size() = " << allAlignables.size() << endl;
706 
707  AlignmentIORoot alignmentIO;
708 
709  for ( vector<string>::iterator itFile = inFileNames.begin(); itFile != inFileNames.end(); ++itFile )
710  {
711  int iter = 1;
712  int ierr = 0;
713 
714  while ( !ierr )
715  {
716  cout << "Read alignment parameters. file / iteration = " << *itFile << " / " << iter << endl;
717 
718  vector< AlignmentParameters* > alignmentParameters =
719  alignmentIO.readAlignmentParameters( allAlignables, (*itFile).c_str(), iter, ierr );
720 
721  cout << "#param / ierr = " << alignmentParameters.size() << " / " << ierr << endl;
722 
723  vector< AlignmentParameters* >::iterator itParam;
724  for ( itParam = alignmentParameters.begin(); itParam != alignmentParameters.end(); ++itParam )
725  alignmentParametersMap[(*itParam)->alignable()].push_back( *itParam );
726 
727  ++iter;
728  }
729  }
730 
731  vector< Alignable* > alignablesToWrite;
732  alignablesToWrite.reserve( alignmentParametersMap.size() );
733 
734  map< Alignable*, vector< AlignmentParameters* > >::iterator itMap;
735  for ( itMap = alignmentParametersMap.begin(); itMap != alignmentParametersMap.end(); ++itMap )
736  {
737  //cout << "merge param for alignable" << itMap->first << endl;
738 
739  AlgebraicVector mergedParam( 6, 0 );
740  AlgebraicSymMatrix mergedCov( 6, 0 );
741  int nMerge = 0;
742 
743  vector< AlignmentParameters* >& vecParam = itMap->second;
744  vector< AlignmentParameters* >::iterator itParam;
745  for ( itParam = vecParam.begin(); itParam != vecParam.end(); ++itParam, ++nMerge )
746  mergedParam += (*itParam)->parameters();
747 
748  mergedParam /= nMerge;
749 
750  // no merging of errors up to now
751  AlignmentParameters* mergedAliParam = vecParam.front()->clone( mergedParam, mergedCov );
752  itMap->first->setAlignmentParameters( mergedAliParam );
753 
754  alignablesToWrite.push_back( itMap->first );
755 
756  if ( applyPar || applyCov ) applyAlignmentParameters( itMap->first, mergedAliParam, applyPar, applyCov );
757  }
758 
759  cout << "alignablesToWrite.size() = " << alignablesToWrite.size() << endl;
760 
761  int ierr = 0;
762  // Write output to "iteration 1", ...
763  alignmentIO.writeAlignmentParameters( alignablesToWrite, outFileName.c_str(), 1, false, ierr );
764  // ... or, if "iteration 1" already exists, write it to "highest iteration + 1"
765  if ( ierr == -1 ) alignmentIO.writeAlignmentParameters( alignablesToWrite, outFileName.c_str(), -1, false, ierr );
766 
767  std::cout << "[KalmanAlignmentAlgorithm::mergeResults] DONE" << std::endl;
768 }
769 
770 
772 {
773  AlignmentPositionError zeroAPE( 0., 0., 0. );
774  theTracker->setAlignmentPositionError( zeroAPE, true );
775 }
776 
void initializeAlignmentSetups(const edm::EventSetup &setup)
virtual char const * what() const
Definition: Exception.cc:141
T getParameter(std::string const &) const
virtual void initialize(const edm::EventSetup &setup, AlignableTracker *tracker, AlignableMuon *muon, AlignableExtras *extras, AlignmentParameterStore *store)
Call at beginning of job (must be implemented in derived class)
T getUntrackedParameter(std::string const &, T const &) const
dictionary parameters
Definition: Parameters.py:2
TrackletCollection refitTracks(const edm::EventSetup &eventSetup, const AlignmentSetupCollection &algoSetups, const ConstTrajTrackPairCollection &tracks, const reco::BeamSpot *beamSpot)
void histogramParameters(std::string histoNamePrefix)
Histogram current estimate of the alignment parameters wrt. the true values.
virtual void rotateInLocalFrame(const RotationType &rotation)
Rotation intepreted in the local reference frame.
Definition: Alignable.cc:91
void getComponents(Alignable *ali, std::vector< Alignable * > &comps) const
static void fillHistogram(std::string histo_name, float data)
const TrajectoryStateUpdator * updator() const
KalmanAlignmentAlgorithm(const edm::ParameterSet &config)
T y() const
Definition: PV3DBase.h:63
virtual void move(const GlobalVector &displacement)=0
Movement with respect to the global reference frame.
PropagationDirection
virtual Alignables components() const =0
Return vector of all direct components.
virtual void setAlignmentPositionError(const AlignmentPositionError &ape, bool propagateDown)
void initializeAlignmentParameters(const edm::EventSetup &setup)
void update(bool enforceUpdate=false)
Call this function in case the associated Alignable was updated by the alignment algorithm.
const ConstTrajTrackPairCollection & trajTrackPairs() const
const TrajectorySmoother * smoother() const
define event information passed to algorithms
AlignmentParameterStore * theParameterStore
const AlgebraicVector & parameters(void) const
Get alignment parameters.
AlignmentUserVariables * userVariables(void) const
Get pointer to user variables.
void clear()
remove all selected Alignables and geometrical restrictions
int iEvent
Definition: GenABIO.cc:230
const PropagationDirection getDirection(const std::string &dir) const
void setValid(bool v)
Set validity flag.
void applyAlignmentParameters(Alignable *ali, AlignmentParameters *par, bool applyPar, bool applyCov) const
T sqrt(T t)
Definition: SSEVec.h:48
T z() const
Definition: PV3DBase.h:64
std::vector< TrackletPtr > TrackletCollection
virtual void run(const edm::EventSetup &setup, const EventInfo &eventInfo)
Run the algorithm (must be implemented in derived class)
virtual void setAlignmentPositionError(const AlignmentPositionError &ape, bool propagateDown)=0
Set the alignment position error - if (!propagateDown) do not affect daughters.
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Definition: Alignable.h:126
const align::Alignables & selectedAlignables() const
vector of alignables selected so far
CLHEP::HepVector AlgebraicVector
AlignmentParameterSelector * theSelector
AlgebraicVector EulerAngles
Definition: Definitions.h:36
void setUserVariables(AlignmentUserVariables *auv)
Set pointer to user variables.
KalmanAlignmentSetup AlignmentSetup
tuple tracks
Definition: testEve_cfg.py:39
const T & get() const
Definition: EventSetup.h:55
void writeAlignmentParameters(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
write AlignmentParameters
const KalmanAlignmentSetup::SortingDirection getSortingDirection(const std::string &sortDir) const
virtual AlignmentParameters * clone(const AlgebraicVector &par, const AlgebraicSymMatrix &cov) const =0
Enforce clone methods in derived classes.
const TrajectoryStateUpdator * updator() const
T const * product() const
Definition: ESHandle.h:86
virtual void terminate(const edm::EventSetup &setup)
Call at end of each loop (must be implemented in derived class)
AlignmentSetupCollection theAlignmentSetups
KalmanAlignmentTrackRefitter * theRefitter
AlgebraicVector translation(void) const
Get translation parameters.
unsigned int addSelection(const std::string &name, const std::vector< char > &paramSel)
virtual const MagneticField * magneticField() const =0
const Propagator * propagator() const
int numberOfUpdates(void) const
Return the number of updates.
const reco::BeamSpot & beamSpot() const
AlignableNavigator * theNavigator
align::GlobalPoints toGlobal(const align::LocalPoints &) const
Return in global coord given a set of local points.
RotationType toMatrix(const EulerAngles &)
Convert rotation angles about x-, y-, z-axes to matrix.
Definition: Utilities.cc:40
CLHEP::HepSymMatrix AlgebraicSymMatrix
static unsigned int const shift
const TrajectoryFitter * fitter() const
static void configure(const edm::ParameterSet &config)
tuple cout
Definition: gather_cfg.py:121
#define DEFINE_EDM_PLUGIN(factory, type, name)
virtual Alignables components() const
Return vector of direct components.
const Propagator * alongPropagator() const
const PropagationDirection oppositeDirection(const PropagationDirection dir) const
volatile std::atomic< bool > shutdown_flag false
align::Parameters readAlignmentParameters(const align::Alignables &alivec, const char *filename, int iter, int &ierr)
read AlignmentParameters
const AlgebraicSymMatrix & covariance(void) const
Get parameter covariance matrix.
AlgebraicVector rotation(void) const
Get rotation parameters.
Constructor of the full muon geometry.
Definition: AlignableMuon.h:36
T x() const
Definition: PV3DBase.h:62
TrajectoryFactoryBase::ReferenceTrajectoryCollection ReferenceTrajectoryCollection
TrajectoryFactoryBase::ExternalPredictionCollection ExternalPredictionCollection
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
T get(const Candidate &c)
Definition: component.h:55
std::vector< ConstTrajTrackPair > ConstTrajTrackPairCollection