CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HIPAlignmentAlgorithm.cc
Go to the documentation of this file.
1 #include <fstream>
2 
3 #include "TFile.h"
4 #include "TTree.h"
5 
10 
13 
30 
35 
36 
38 
39 // using namespace std;
40 
41 // Constructor ----------------------------------------------------------------
42 
45 {
46 
47  // parse parameters
48 
49  verbose = cfg.getParameter<bool>("verbosity");
50 
51  outpath = cfg.getParameter<std::string>("outpath");
52  outfile = cfg.getParameter<std::string>("outfile");
53  outfile2 = cfg.getParameter<std::string>("outfile2");
54  struefile = cfg.getParameter<std::string>("trueFile");
55  smisalignedfile = cfg.getParameter<std::string>("misalignedFile");
56  salignedfile = cfg.getParameter<std::string>("alignedFile");
57  siterationfile = cfg.getParameter<std::string>("iterationFile");
58  suvarfile = cfg.getParameter<std::string>("uvarFile");
59  sparameterfile = cfg.getParameter<std::string>("parameterFile");
60  ssurveyfile = cfg.getParameter<std::string>("surveyFile");
61 
62  outfile =outpath+outfile;//Eventwise tree
63  outfile2 =outpath+outfile2;//Alignablewise tree
64  struefile =outpath+struefile;
65  smisalignedfile=outpath+smisalignedfile;
66  salignedfile =outpath+salignedfile;
67  siterationfile =outpath+siterationfile;
68  suvarfile =outpath+suvarfile;
69  sparameterfile =outpath+sparameterfile;
70  ssurveyfile =outpath+ssurveyfile;
71 
72  // parameters for APE
73  theApplyAPE = cfg.getParameter<bool>("applyAPE");
74  theAPEParameterSet = cfg.getParameter<std::vector<edm::ParameterSet> >("apeParam");
75 
76  theMaxAllowedHitPull = cfg.getParameter<double>("maxAllowedHitPull");
77  theMinimumNumberOfHits = cfg.getParameter<int>("minimumNumberOfHits");
78  theMaxRelParameterError = cfg.getParameter<double>("maxRelParameterError");
79 
80  // for collector mode (parallel processing)
81  isCollector=cfg.getParameter<bool>("collectorActive");
82  theCollectorNJobs=cfg.getParameter<int>("collectorNJobs");
83  theCollectorPath=cfg.getParameter<std::string>("collectorPath");
84  theFillTrackMonitoring=cfg.getUntrackedParameter<bool>("fillTrackMonitoring");
85 
86  if (isCollector) edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Collector mode";
87 
88  theEventPrescale = cfg.getParameter<int>("eventPrescale");
90 
91  for (std::string &s : cfg.getUntrackedParameter<std::vector<std::string> >("surveyResiduals")) {
93  }
94 
95  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] constructed.";
96 
97 }
98 
99 // Call at beginning of job ---------------------------------------------------
100 
101 void
104  AlignmentParameterStore* store )
105 {
106  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Initializing...";
107 
108  edm::ESHandle<Alignments> globalPositionRcd;
109  // FIXME! temporary solution to get highest possible run number
110  const unsigned int MAX_VAL(std::numeric_limits<unsigned int>::max());
112  if (iov.first().eventID().run()!=1 || iov.last().eventID().run()!=MAX_VAL) {
113  throw cms::Exception("DatabaseError")
114  << "@SUB=AlignmentProducer::applyDB"
115  << "\nTrying to apply "<< setup.get<GlobalPositionRcd>().key().name()
116  << " with multiple IOVs in tag.\n"
117  << "Validity range is "
118  << iov.first().eventID().run() << " - " << iov.last().eventID().run();
119  }
120 
121  // accessor Det->AlignableDet
122  if ( !muon )
124  else if ( !tracker )
126  else
127  theAlignableDetAccessor = new AlignableNavigator(tracker, muon);
128 
129  // set alignmentParameterStore
131 
132  // get alignables
134 
135  // clear theAPEParameters, if necessary
136  theAPEParameters.clear();
137 
138  // get APE parameters
139  if(theApplyAPE){
140  AlignmentParameterSelector selector(tracker, muon);
141  for (std::vector<edm::ParameterSet>::const_iterator setiter = theAPEParameterSet.begin();
142  setiter != theAPEParameterSet.end();
143  ++setiter) {
144  std::vector<Alignable*> alignables;
145 
146  selector.clear();
147  edm::ParameterSet selectorPSet = setiter->getParameter<edm::ParameterSet>("Selector");
148  std::vector<std::string> alignParams = selectorPSet.getParameter<std::vector<std::string> >("alignParams");
149  if (alignParams.size() == 1 && alignParams[0] == std::string("selected")) {
150  alignables = theAlignables;
151  }
152  else {
153  selector.addSelections(selectorPSet);
154  alignables = selector.selectedAlignables();
155  }
156 
157  std::vector<double> apeSPar = setiter->getParameter<std::vector<double> >("apeSPar");
158  std::vector<double> apeRPar = setiter->getParameter<std::vector<double> >("apeRPar");
159  std::string function = setiter->getParameter<std::string>("function");
160 
161  if (apeSPar.size() != 3 || apeRPar.size() != 3)
162  throw cms::Exception("BadConfig") << "apeSPar and apeRPar must have 3 values each" << std::endl;
163 
164  for (std::vector<double>::const_iterator i = apeRPar.begin(); i != apeRPar.end(); ++i) {
165  apeSPar.push_back(*i);
166  }
167 
168  if (function == std::string("linear")) {
169  apeSPar.push_back(0.); // c.f. note in calcAPE
170  }
171  else if (function == std::string("exponential")) {
172  apeSPar.push_back(1.); // c.f. note in calcAPE
173  }
174  else if (function == std::string("step")) {
175  apeSPar.push_back(2.); // c.f. note in calcAPE
176  }
177  else {
178  throw cms::Exception("BadConfig") << "APE function must be \"linear\" or \"exponential\"." << std::endl;
179  }
180 
181  theAPEParameters.push_back(std::pair<std::vector<Alignable*>, std::vector<double> >(alignables, apeSPar));
182  }
183  }
184 }
185 
186 // Call at new loop -------------------------------------------------------------
188 {
189 
190  // iterate over all alignables and attach user variables
191  for( std::vector<Alignable*>::const_iterator it=theAlignables.begin();
192  it!=theAlignables.end(); it++ )
193  {
194  AlignmentParameters* ap = (*it)->alignmentParameters();
195  int npar=ap->numSelected();
196  HIPUserVariables* userpar = new HIPUserVariables(npar);
197  ap->setUserVariables(userpar);
198  }
199 
200  // try to read in alignment parameters from a previous iteration
201  AlignablePositions theAlignablePositionsFromFile =
203  salignedfile.c_str(),-1,ioerr);
204 
205  int numAlignablesFromFile = theAlignablePositionsFromFile.size();
206 
207  if (numAlignablesFromFile==0) { // file not there: first iteration
208 
209  // set iteration number to 1
210  if (isCollector) theIteration=0;
211  else theIteration=1;
212  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] File not found => iteration "<<theIteration;
213 
214  // get true (de-misaligned positions) and write to root file
215  // hardcoded iteration=1
217  struefile.c_str(),1,false,ioerr);
218 
219  // get misaligned positions and write to root file
220  // hardcoded iteration=1
222  smisalignedfile.c_str(),1,false,ioerr);
223 
224  }
225 
226  else { // there have been previous iterations
227 
228  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Alignables Read "
229  << numAlignablesFromFile;
230 
231  // get iteration number from file
233 
234  // increase iteration
235  theIteration++;
236  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Iteration increased by one!";
237 
238  // now apply psotions of file from prev iteration
239  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Apply positions from file ...";
241  theAlignablePositionsFromFile,ioerr);
242 
243  }
244 
245  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm] Current Iteration number: "
246  << theIteration;
247 
248 
249  // book root trees
250  bookRoot();
251 
252 
253  /*---------------------moved to terminate------------------------------
254  if (theLevels.size() > 0)
255  {
256  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Using survey constraint";
257 
258  unsigned int nAlignable = theAlignables.size();
259 
260  for (unsigned int i = 0; i < nAlignable; ++i)
261  {
262  const Alignable* ali = theAlignables[i];
263 
264  AlignmentParameters* ap = ali->alignmentParameters();
265 
266  HIPUserVariables* uservar =
267  dynamic_cast<HIPUserVariables*>(ap->userVariables());
268 
269  for (unsigned int l = 0; l < theLevels.size(); ++l)
270  {
271  SurveyResidual res(*ali, theLevels[l], true);
272 
273  if ( res.valid() )
274  {
275  AlgebraicSymMatrix invCov = res.inverseCovariance();
276 
277  // variable for tree
278  AlgebraicVector sensResid = res.sensorResidual();
279  m3_Id = ali->id();
280  m3_ObjId = theLevels[l];
281  m3_par[0] = sensResid[0]; m3_par[1] = sensResid[1]; m3_par[2] = sensResid[2];
282  m3_par[3] = sensResid[3]; m3_par[4] = sensResid[4]; m3_par[5] = sensResid[5];
283 
284  uservar->jtvj += invCov;
285  uservar->jtve += invCov * sensResid;
286 
287  theTree3->Fill();
288  }
289  }
290 
291  // align::LocalVectors residuals = res1.pointsResidual();
292 
293  // unsigned int nPoints = residuals.size();
294 
295  // for (unsigned int k = 0; k < nPoints; ++k)
296  // {
297  // AlgebraicMatrix J = term->survey()->derivatives(k);
298  // AlgebraicVector e(3); // local residual
299 
300  // const align::LocalVector& lr = residuals[k];
301 
302  // e(1) = lr.x(); e(2) = lr.y(); e(3) = lr.z();
303 
304  // uservar->jtvj += invCov1.similarity(J);
305  // uservar->jtve += J * (invCov1 * e);
306  // }
307 
308  }
309  }
310  //------------------------------------------------------*/
311 
312  // set alignment position error
314 
315  // run collector job if we are in parallel mode
316  if (isCollector) collector();
317 
318 }
319 
320 // Call at end of job ---------------------------------------------------------
321 
323 {
324 
325  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Terminating";
326 
327  // calculating survey residuals
328  if (theLevels.size() > 0)
329  {
330  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Using survey constraint";
331 
332  unsigned int nAlignable = theAlignables.size();
333 
334  for (unsigned int i = 0; i < nAlignable; ++i)
335  {
336  const Alignable* ali = theAlignables[i];
337 
339 
340  HIPUserVariables* uservar =
341  dynamic_cast<HIPUserVariables*>(ap->userVariables());
342 
343  for (unsigned int l = 0; l < theLevels.size(); ++l)
344  {
345  SurveyResidual res(*ali, theLevels[l], true);
346 
347  if ( res.valid() )
348  {
349  AlgebraicSymMatrix invCov = res.inverseCovariance();
350 
351  // variable for tree
352  AlgebraicVector sensResid = res.sensorResidual();
353  m3_Id = ali->id();
354  m3_ObjId = theLevels[l];
355  m3_par[0] = sensResid[0]; m3_par[1] = sensResid[1]; m3_par[2] = sensResid[2];
356  m3_par[3] = sensResid[3]; m3_par[4] = sensResid[4]; m3_par[5] = sensResid[5];
357 
358  uservar->jtvj += invCov;
359  uservar->jtve += invCov * sensResid;
360 
361  theTree3->Fill();
362  }
363  }
364 
365  // align::LocalVectors residuals = res1.pointsResidual();
366 
367  // unsigned int nPoints = residuals.size();
368 
369  // for (unsigned int k = 0; k < nPoints; ++k)
370  // {
371  // AlgebraicMatrix J = term->survey()->derivatives(k);
372  // AlgebraicVector e(3); // local residual
373 
374  // const align::LocalVector& lr = residuals[k];
375 
376  // e(1) = lr.x(); e(2) = lr.y(); e(3) = lr.z();
377 
378  // uservar->jtvj += invCov1.similarity(J);
379  // uservar->jtve += J * (invCov1 * e);
380  // }
381 
382  }
383  }
384 
385  // write user variables
388  theIteration,false,ioerr);
389 
390  // now calculate alignment corrections ...
391  int ialigned=0;
392  // iterate over alignment parameters
393  for(std::vector<Alignable*>::const_iterator
394  it=theAlignables.begin(); it!=theAlignables.end(); it++) {
395  Alignable* ali=(*it);
396  // Alignment parameters
398 
399  // try to calculate parameters
400  bool test = calcParameters(ali);
401 
402  // if successful, apply parameters
403  if (test) {
404  edm::LogInfo("Alignment") << "now apply params";
406  // set these parameters 'valid'
407  ali->alignmentParameters()->setValid(true);
408  // increase counter
409  ialigned++;
410  }
411  else par->setValid(false);
412  }
413  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::terminate] Aligned units: " << ialigned;
414 
415  // fill alignable wise root tree
416  fillRoot(iSetup);
417 
418  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm] Writing aligned parameters to file: " << theAlignables.size();
419 
420  // write new absolute positions to disk
422  salignedfile.c_str(),theIteration,false,ioerr);
423 
424  // write alignment parameters to disk
426  sparameterfile.c_str(),theIteration,false,ioerr);
427 
428  // write iteration number to file
430 
431  // write out trees and close root file
432 
433  // eventwise tree
434  theFile->cd();
435  theTree->Write();
436  delete theFile;
437 
438  if (theLevels.size() > 0){
439  theFile3->cd();
440  theTree3->Write();
441  delete theFile3;
442  }
443 
444  // alignable-wise tree is only filled once
445  if (theIteration==1) { // only for 1st iteration
446  theFile2->cd();
447  theTree2->Write();
448  delete theFile2;
449  }
450 
451 }
452 
454  const Alignable* ali,
455  const TrajectoryStateOnSurface & tsos,
457 {
458  static const unsigned int hitDim = 1;
459 
460  // get trajectory impact point
461  LocalPoint alvec = tsos.localPosition();
462  AlgebraicVector pos(hitDim);
463  pos[0] = alvec.x();
464 
465  // get impact point covariance
466  AlgebraicSymMatrix ipcovmat(hitDim);
467  ipcovmat[0][0] = tsos.localError().positionError().xx();
468 
469  // get hit local position and covariance
470  AlgebraicVector coor(hitDim);
471  coor[0] = hit->localPosition().x();
472 
473  AlgebraicSymMatrix covmat(hitDim);
474  covmat[0][0] = hit->localPositionError().xx();
475 
476  // add hit and impact point covariance matrices
477  covmat = covmat + ipcovmat;
478 
479  // calculate the x pull of this hit
480  double xpull = 0.;
481  if (covmat[0][0] != 0.) xpull = (pos[0] - coor[0])/sqrt(fabs(covmat[0][0]));
482 
483  // get Alignment Parameters
484  AlignmentParameters* params = ali->alignmentParameters();
485  // get derivatives
486  AlgebraicMatrix derivs2D = params->selectedDerivatives(tsos, alidet);
487  // calculate user parameters
488  int npar = derivs2D.num_row();
489  // std::cout << "Dumping derivsSEL: \n" << derivs2D <<std::endl;
490  AlgebraicMatrix derivs(npar, hitDim, 0);
491 
492  for (int ipar=0;ipar<npar;ipar++) {
493  for (unsigned int idim=0;idim<hitDim;idim++) {
494  derivs[ipar][idim] = derivs2D[ipar][idim];
495  }
496  }
497  // std::cout << "Dumping derivs: \n" << derivs << std::endl;
498  // invert covariance matrix
499  int ierr;
500  // if (covmat[1][1]>4.0) nhitDim = hitDim;
501 
502  covmat.invert(ierr);
503  if (ierr != 0) {
504  edm::LogError("Alignment") << "Matrix inversion failed!";
505  return false;
506  }
507 
508  bool useThisHit = (theMaxAllowedHitPull <= 0.);
509 
510  useThisHit |= (fabs(xpull) < theMaxAllowedHitPull);
511 
512  // bailing out
513  if (!useThisHit) return false;
514 
515  // std::cout << "We use this hit in " << subDet << std::endl;
516  // std::cout << "Npars= " << npar << " NhitDim=1 Size of derivs: "
517  // << derivs.num_row() << " x " << derivs.num_col() << std::endl;
518 
519  // AlgebraicMatrix jtvjtmp();
520  AlgebraicMatrix covtmp(covmat);
521  // std::cout << "Preapring JTVJTMP -> " << std::flush;
522  AlgebraicMatrix jtvjtmp(derivs * covtmp *derivs.T());
523  // std::cout << "TMP JTVJ= \n" << jtvjtmp << std::endl;
524  AlgebraicSymMatrix thisjtvj(npar);
525  AlgebraicVector thisjtve(npar);
526  // std::cout << "Preparing ThisJtVJ" << std::flush;
527  // thisjtvj = covmat.similarity(derivs);
528  thisjtvj.assign(jtvjtmp);
529  // std::cout<<" ThisJtVE"<<std::endl;
530  thisjtve=derivs * covmat * (pos-coor);
531 
532  AlgebraicVector hitresidual(hitDim);
533  hitresidual[0] = (pos[0] - coor[0]);
534 
535  AlgebraicMatrix hitresidualT;
536  hitresidualT = hitresidual.T();
537  // std::cout << "HitResidualT = \n" << hitresidualT << std::endl;
538 
539  // access user variables (via AlignmentParameters)
540  HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(params->userVariables());
541  uservar->jtvj += thisjtvj;
542  uservar->jtve += thisjtve;
543  uservar->nhit++;
544  // The following variable is needed for the extended 1D/2D hit fix using
545  // matrix shrinkage and expansion
546  // uservar->hitdim = hitDim;
547 
548  //for alignable chi squared
549  float thischi2;
550  thischi2 = (hitresidualT *covmat *hitresidual)[0];
551 
552  if ( verbose && ((thischi2/ static_cast<float>(uservar->nhit)) >10.0) ) {
553  edm::LogWarning("Alignment") << "Added to Chi2 the number " << thischi2 <<" having "
554  << uservar->nhit << " dof " << std::endl << "X-resid "
555  << hitresidual[0] << " Y-resid "
556  << hitresidual[1] << std::endl << " Cov^-1 matr (covmat): [0][0]= "
557  << covmat[0][0] << " [0][1]= "
558  << covmat[0][1] << " [1][0]= "
559  << covmat[1][0] << " [1][1]= "
560  << covmat[1][1] << std::endl;
561  }
562 
563  uservar->alichi2 += thischi2; // a bit weird, but vector.transposed * matrix * vector doesn't give a double in CMSSW's opinion
564  uservar->alindof += hitDim; // 2D hits contribute twice to the ndofs
565 
566  return true;
567 }
568 
570  const Alignable* ali,
571  const TrajectoryStateOnSurface & tsos,
573 {
574  static const unsigned int hitDim = 2;
575 
576  // get trajectory impact point
577  LocalPoint alvec = tsos.localPosition();
578  AlgebraicVector pos(hitDim);
579  pos[0] = alvec.x();
580  pos[1] = alvec.y();
581 
582  // get impact point covariance
583  AlgebraicSymMatrix ipcovmat(hitDim);
584  ipcovmat[0][0] = tsos.localError().positionError().xx();
585  ipcovmat[1][1] = tsos.localError().positionError().yy();
586  ipcovmat[0][1] = tsos.localError().positionError().xy();
587 
588  // get hit local position and covariance
589  AlgebraicVector coor(hitDim);
590  coor[0] = hit->localPosition().x();
591  coor[1] = hit->localPosition().y();
592 
593  AlgebraicSymMatrix covmat(hitDim);
594  covmat[0][0] = hit->localPositionError().xx();
595  covmat[1][1] = hit->localPositionError().yy();
596  covmat[0][1] = hit->localPositionError().xy();
597 
598  // add hit and impact point covariance matrices
599  covmat = covmat + ipcovmat;
600 
601  // calculate the x pull and y pull of this hit
602  double xpull = 0.;
603  double ypull = 0.;
604  if (covmat[0][0] != 0.) xpull = (pos[0] - coor[0])/sqrt(fabs(covmat[0][0]));
605  if (covmat[1][1] != 0.) ypull = (pos[1] - coor[1])/sqrt(fabs(covmat[1][1]));
606 
607  // get Alignment Parameters
608  AlignmentParameters* params = ali->alignmentParameters();
609  // get derivatives
610  AlgebraicMatrix derivs2D = params->selectedDerivatives(tsos, alidet);
611  // calculate user parameters
612  int npar = derivs2D.num_row();
613  // std::cout << "Dumping derivsSEL: \n"<< derivs2D <<std::endl;
614  AlgebraicMatrix derivs(npar, hitDim, 0);
615 
616  for (int ipar=0;ipar<npar;ipar++) {
617  for (unsigned int idim=0;idim<hitDim;idim++) {
618  derivs[ipar][idim] = derivs2D[ipar][idim];
619  }
620  }
621  // std::cout << "Dumping derivs: \n" << derivs << std::endl;
622  // invert covariance matrix
623  int ierr;
624  // if (covmat[1][1]>4.0) nhitDim = hitDim;
625 
626  covmat.invert(ierr);
627  if (ierr != 0) {
628  edm::LogError("Alignment") << "Matrix inversion failed!";
629  return false;
630  }
631 
632  bool useThisHit = (theMaxAllowedHitPull <= 0.);
633 
634  useThisHit |= (fabs(xpull) < theMaxAllowedHitPull && fabs(ypull) < theMaxAllowedHitPull);
635 
636  // bailing out
637  if (!useThisHit) return false;
638 
639  // std::cout << "We use this hit in " << subDet << std::endl;
640  // std::cout << "Npars= " << npar << " NhitDim=2 Size of derivs: "
641  // << derivs.num_row() << " x " << derivs.num_col() << std::endl;
642 
643  // AlgebraicMatrix jtvjtmp();
644  AlgebraicMatrix covtmp(covmat);
645  // std::cout << "Preapring JTVJTMP -> " << std::flush;
646  AlgebraicMatrix jtvjtmp(derivs * covtmp *derivs.T());
647  // std::cout << "TMP JTVJ= \n" << jtvjtmp << std::endl;
648  AlgebraicSymMatrix thisjtvj(npar);
649  AlgebraicVector thisjtve(npar);
650  // std::cout << "Preparing ThisJtVJ" << std::flush;
651  // thisjtvj = covmat.similarity(derivs);
652  thisjtvj.assign(jtvjtmp);
653  // std::cout<<" ThisJtVE"<<std::endl;
654  thisjtve=derivs * covmat * (pos-coor);
655 
656  AlgebraicVector hitresidual(hitDim);
657  hitresidual[0] = (pos[0] - coor[0]);
658  hitresidual[1] = (pos[1] - coor[1]);
659  // if(nhitDim>1) {
660  // hitresidual[1] =0.0;
661  // nhitDim=1;
662  // }
663 
664  AlgebraicMatrix hitresidualT;
665  hitresidualT = hitresidual.T();
666  // std::cout << "HitResidualT = \n" << hitresidualT << std::endl;
667  // access user variables (via AlignmentParameters)
668  HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(params->userVariables());
669  uservar->jtvj += thisjtvj;
670  uservar->jtve += thisjtve;
671  uservar->nhit++;
672  // The following variable is needed for the extended 1D/2D hit fix using
673  // matrix shrinkage and expansion
674  // uservar->hitdim = hitDim;
675 
676  //for alignable chi squared
677  float thischi2;
678  thischi2 = (hitresidualT *covmat *hitresidual)[0];
679 
680  if ( verbose && ((thischi2/ static_cast<float>(uservar->nhit)) >10.0) ) {
681  edm::LogWarning("Alignment") << "Added to Chi2 the number " << thischi2 <<" having "
682  << uservar->nhit << " dof " << std::endl << "X-resid "
683  << hitresidual[0] << " Y-resid "
684  << hitresidual[1] << std::endl << " Cov^-1 matr (covmat): [0][0]= "
685  << covmat[0][0] << " [0][1]= "
686  << covmat[0][1] << " [1][0]= "
687  << covmat[1][0] << " [1][1]= "
688  << covmat[1][1] << std::endl;
689  }
690 
691  uservar->alichi2 += thischi2; // a bit weird, but vector.transposed * matrix * vector doesn't give a double in CMSSW's opinion
692  uservar->alindof += hitDim; // 2D hits contribute twice to the ndofs
693 
694  return true;
695 }
696 
697 // Run the algorithm on trajectories and tracks -------------------------------
699 {
700  if (isCollector) return;
701 
702  TrajectoryStateCombiner tsoscomb;
703 
704  // AM: not really needed
705  // AM: m_Ntracks = 0 should be sufficient
706  int itr=0;
707  m_Ntracks=0;
708  for(itr=0;itr<MAXREC;++itr){
709  m_Nhits[itr]=0;
710  m_Pt[itr]=-5.0;
711  m_P[itr]=-5.0;
712  m_nhPXB[itr]=0;
713  m_nhPXF[itr]=0;
714  m_nhTIB[itr]=0;
715  m_nhTOB[itr]=0;
716  m_nhTIB[itr]=0;
717  m_nhTEC[itr]=0;
718  m_Eta[itr]=-99.0;
719  m_Phi[itr]=-4.0;
720  m_Chi2n[itr]=-11.0;
721  m_d0[itr]=-999;
722  m_dz[itr]=-999;
723  }
724  itr=0;
725 
726  // AM: what is this needed for?
727  theFile->cd();
728 
729  // loop over tracks
731  for (ConstTrajTrackPairCollection::const_iterator it=tracks.begin();
732  it!=tracks.end();
733  ++it) {
734 
735  const Trajectory* traj = (*it).first;
736  const reco::Track* track = (*it).second;
737 
738  float pt = track->pt();
739  float eta = track->eta();
740  float phi = track->phi();
741  float p = track->p();
742  float chi2n = track->normalizedChi2();
743  int nhit = track->numberOfValidHits();
744  float d0 = track->d0();
745  float dz = track->dz();
746 
747  int nhpxb = track->hitPattern().numberOfValidPixelBarrelHits();
748  int nhpxf = track->hitPattern().numberOfValidPixelEndcapHits();
749  int nhtib = track->hitPattern().numberOfValidStripTIBHits();
750  int nhtob = track->hitPattern().numberOfValidStripTOBHits();
751  int nhtid = track->hitPattern().numberOfValidStripTIDHits();
752  int nhtec = track->hitPattern().numberOfValidStripTECHits();
753 
754  if (verbose) edm::LogInfo("Alignment") << "New track pt,eta,phi,chi2n,hits: "
755  << pt << ","
756  << eta << ","
757  << phi << ","
758  << chi2n << ","
759  << nhit;
760  // edm::LogWarning("Alignment") << "New track pt,eta,phi,chi2n,hits: "
761  // << pt << ","
762  // << eta << ","
763  // << phi << ","
764  // << chi2n << ","
765  // << nhit;
766 
767  // fill track parameters in root tree
768  if (itr<MAXREC) {
769  m_Nhits[itr]=nhit;
770  m_Pt[itr]=pt;
771  m_P[itr]=p;
772  m_Eta[itr]=eta;
773  m_Phi[itr]=phi;
774  m_Chi2n[itr]=chi2n;
775  m_nhPXB[itr]=nhpxb;
776  m_nhPXF[itr]=nhpxf;
777  m_nhTIB[itr]=nhtib;
778  m_nhTOB[itr]=nhtob;
779  m_nhTID[itr]=nhtid;
780  m_nhTEC[itr]=nhtec;
781  m_d0[itr]=d0;
782  m_dz[itr]=dz;
783  itr++;
784  m_Ntracks=itr;
785  }
786  // AM: Can be simplified
787 
788  std::vector<const TransientTrackingRecHit*> hitvec;
789  std::vector<TrajectoryStateOnSurface> tsosvec;
790 
791  // loop over measurements
792  std::vector<TrajectoryMeasurement> measurements = traj->measurements();
793  for (std::vector<TrajectoryMeasurement>::iterator im=measurements.begin();
794  im!=measurements.end();
795  ++im) {
796 
797  TrajectoryMeasurement meas = *im;
798 
799  // const TransientTrackingRecHit* ttrhit = &(*meas.recHit());
800  // const TrackingRecHit *hit = ttrhit->hit();
801  const TransientTrackingRecHit* hit = &(*meas.recHit());
802 
804 
805  // this is the updated state (including the current hit)
806  // TrajectoryStateOnSurface tsos=meas.updatedState();
807  // combine fwd and bwd predicted state to get state
808  // which excludes current hit
809 
811  bool skiphit = false;
812  if (eventInfo.clusterValueMap()) {
813  // check from the PrescalingMap if the hit was taken.
814  // If not skip to the next TM
815  // bool hitTaken=false;
816  AlignmentClusterFlag myflag;
817 
818  int subDet = hit->geographicalId().subdetId();
819  //take the actual RecHit out of the Transient one
820  const TrackingRecHit *rechit=hit->hit();
821  if (subDet>2) { // AM: if possible use enum instead of hard-coded value
822  const std::type_info &type = typeid(*rechit);
823 
824  if (type == typeid(SiStripRecHit1D)) {
825 
826  const SiStripRecHit1D* stripHit1D = dynamic_cast<const SiStripRecHit1D*>(rechit);
827  if (stripHit1D) {
828  SiStripRecHit1D::ClusterRef stripclust(stripHit1D->cluster());
829  // myflag=PrescMap[stripclust];
830  myflag = (*eventInfo.clusterValueMap())[stripclust];
831  } else {
832  edm::LogError("HIPAlignmentAlgorithm")
833  << "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Strip RecHit failed! "
834  << "TypeId of the RecHit: " << className(*hit) <<std::endl;
835  }
836 
837  }//end if type = SiStripRecHit1D
838  else if(type == typeid(SiStripRecHit2D)){
839 
840  const SiStripRecHit2D* stripHit2D = dynamic_cast<const SiStripRecHit2D*>(rechit);
841  if (stripHit2D) {
842  SiStripRecHit2D::ClusterRef stripclust(stripHit2D->cluster());
843  // myflag=PrescMap[stripclust];
844  myflag = (*eventInfo.clusterValueMap())[stripclust];
845  } else {
846  edm::LogError("HIPAlignmentAlgorithm")
847  << "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Strip RecHit failed! "
848  // << "TypeId of the TTRH: " << className(*ttrhit) << std::endl;
849  << "TypeId of the TTRH: " << className(*hit) << std::endl;
850  }
851  } //end if type == SiStripRecHit2D
852  } //end if hit from strips
853  else {
854  const SiPixelRecHit* pixelhit= dynamic_cast<const SiPixelRecHit*>(rechit);
855  if (pixelhit) {
856  SiPixelClusterRefNew pixelclust(pixelhit->cluster());
857  // myflag=PrescMap[pixelclust];
858  myflag = (*eventInfo.clusterValueMap())[pixelclust];
859  }
860  else {
861  edm::LogError("HIPAlignmentAlgorithm")
862  << "ERROR in <HIPAlignmentAlgorithm::run>: Dynamic cast of Pixel RecHit failed! "
863  // << "TypeId of the TTRH: " << className(*ttrhit) << std::endl;
864  << "TypeId of the TTRH: " << className(*hit) << std::endl;
865  }
866  } //end 'else' it is a pixel hit
867  // bool hitTaken=myflag.isTaken();
868  if (!myflag.isTaken()) {
869  skiphit=true;
870  //std::cout<<"Hit from subdet "<<rechit->geographicalId().subdetId()<<" prescaled out."<<std::endl;
871  continue;
872  }
873  }//end if Prescaled Hits
875  if (skiphit) {
876  throw cms::Exception("LogicError")
877  << "ERROR in <HIPAlignmentAlgorithm::run>: this hit should have been skipped!"
878  << std::endl;
879  }
880 
882  meas.backwardPredictedState());
883 
884  if(tsos.isValid()){
885  // hitvec.push_back(ttrhit);
886  hitvec.push_back(hit);
887  tsosvec.push_back(tsos);
888  }
889 
890  } //hit valid
891  }
892 
893  // transform RecHit vector to AlignableDet vector
894  std::vector <AlignableDetOrUnitPtr> alidetvec = theAlignableDetAccessor->alignablesFromHits(hitvec);
895 
896  // get concatenated alignment parameters for list of alignables
898 
899  std::vector<TrajectoryStateOnSurface>::const_iterator itsos=tsosvec.begin();
900  std::vector<const TransientTrackingRecHit*>::const_iterator ihit=hitvec.begin();
901 
902  // loop over vectors(hit,tsos)
903  while (itsos != tsosvec.end()) {
904 
905  // get AlignableDet for this hit
906  const GeomDet* det = (*ihit)->det();
907  // int subDet= (*ihit)->geographicalId().subdetId();
908  uint32_t nhitDim = (*ihit)->dimension();
909 
911 
912  // get relevant Alignable
913  Alignable* ali = aap.alignableFromAlignableDet(alidet);
914 
915  if (ali!=0) {
916  if (nhitDim==1) {
917  processHit1D(alidet, ali, *itsos, *ihit);
918  } else if (nhitDim==2) {
919  processHit2D(alidet, ali, *itsos, *ihit);
920  }
921  }
922 
923  itsos++;
924  ihit++;
925  }
926  } // end of track loop
927 
928  // fill eventwise root tree (with prescale defined in pset)
931  if (theCurrentPrescale<=0) {
932  theTree->Fill();
934  }
935  }
936 }
937 
938 // ----------------------------------------------------------------------------
939 
941 {
942  int result;
943 
944  std::ifstream inIterFile(filename.c_str(), std::ios::in);
945  if (!inIterFile) {
946  edm::LogError("Alignment") << "[HIPAlignmentAlgorithm::readIterationFile] ERROR! "
947  << "Unable to open Iteration file";
948  result = -1;
949  }
950  else {
951  inIterFile >> result;
952  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::readIterationFile] "
953  << "Read last iteration number from file: " << result;
954  }
955  inIterFile.close();
956 
957  return result;
958 }
959 
960 // ----------------------------------------------------------------------------
961 
963 {
964  std::ofstream outIterFile((filename.c_str()), std::ios::out);
965  if (!outIterFile) {
966  edm::LogError("Alignment") << "[HIPAlignmentAlgorithm::writeIterationFile] ERROR: Unable to write Iteration file";
967  }
968  else {
969  outIterFile << iter;
970  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::writeIterationFile] writing iteration number to file: " << iter;
971  }
972  outIterFile.close();
973 }
974 
975 
976 // ----------------------------------------------------------------------------
977 // set alignment position error
978 
980 {
981 
982 
983  // Check if user wants to override APE
984  if ( !theApplyAPE )
985  {
986  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::setAlignmentPositionError] No APE applied";
987  return; // NO APE APPLIED
988  }
989 
990 
991  edm::LogWarning("Alignment") <<"[HIPAlignmentAlgorithm::setAlignmentPositionError] Apply APE!";
992 
993  double apeSPar[3], apeRPar[3];
994  for (std::vector<std::pair<std::vector<Alignable*>, std::vector<double> > >::const_iterator alipars = theAPEParameters.begin();
995  alipars != theAPEParameters.end();
996  ++alipars) {
997  const std::vector<Alignable*> &alignables = alipars->first;
998  const std::vector<double> &pars = alipars->second;
999 
1000  apeSPar[0] = pars[0];
1001  apeSPar[1] = pars[1];
1002  apeSPar[2] = pars[2];
1003  apeRPar[0] = pars[3];
1004  apeRPar[1] = pars[4];
1005  apeRPar[2] = pars[5];
1006 
1007  double function = pars[6];
1008 
1009  // Printout for debug
1010  printf("APE: %u alignables\n", (unsigned int)alignables.size());
1011  for ( int i=0; i<21; ++i ) {
1012  double apelinstest=calcAPE(apeSPar,i,0.);
1013  double apeexpstest=calcAPE(apeSPar,i,1.);
1014  double apelinrtest=calcAPE(apeRPar,i,0.);
1015  double apeexprtest=calcAPE(apeRPar,i,1.);
1016  printf("APE: iter slin sexp rlin rexp: %5d %12.5f %12.5f %12.5f %12.5f\n",
1017  i,apelinstest,apeexpstest,apelinrtest,apeexprtest);
1018  }
1019 
1020  // set APE
1021  double apeshift=calcAPE(apeSPar,theIteration,function);
1022  double aperot =calcAPE(apeRPar,theIteration,function);
1023  theAlignmentParameterStore->setAlignmentPositionError( alignables, apeshift, aperot );
1024  }
1025 
1026 }
1027 
1028 // ----------------------------------------------------------------------------
1029 // calculate APE
1030 
1031 double
1032 HIPAlignmentAlgorithm::calcAPE(double* par, int iter, double function)
1033 {
1034  double diter=(double)iter;
1035 
1036  // The following floating-point equality check is safe because this
1037  // "0." and this "1." are generated by the compiler, in the very
1038  // same file. Whatever approximization scheme it uses to turn "1."
1039  // into 0.9999999999998 in the HIPAlignmentAlgorithm::initialize is
1040  // also used here. If I'm wrong, you'll get an assertion.
1041  if (function == 0.) {
1042  return std::max(par[1],par[0]+((par[1]-par[0])/par[2])*diter);
1043  }
1044  else if (function == 1.) {
1045  return std::max(0.,par[0]*(exp(-pow(diter,par[1])/par[2])));
1046  }
1047  else if (function == 2.) {
1048  int ipar2 = (int) par[2];
1049  int step = iter/ipar2;
1050  double dstep = (double) step;
1051  return std::max(0.0, par[0] - par[1]*dstep);
1052  }
1053  else assert(false); // should have been caught in the constructor
1054 }
1055 
1056 
1057 // ----------------------------------------------------------------------------
1058 // book root trees
1059 
1061 {
1062  // create ROOT files
1063  theFile = new TFile(outfile.c_str(),"update");
1064  theFile->cd();
1065 
1066  // book event-wise ROOT Tree
1067 
1068  TString tname="T1";
1069  char iterString[5];
1070  snprintf(iterString, sizeof(iterString), "%i",theIteration);
1071  tname.Append("_");
1072  tname.Append(iterString);
1073 
1074  theTree = new TTree(tname,"Eventwise tree");
1075 
1076  //theTree->Branch("Run", &m_Run, "Run/I");
1077  //theTree->Branch("Event", &m_Event, "Event/I");
1078  theTree->Branch("Ntracks", &m_Ntracks, "Ntracks/I");
1079  theTree->Branch("Nhits", m_Nhits, "Nhits[Ntracks]/I");
1080  theTree->Branch("nhPXB", m_nhPXB, "nhPXB[Ntracks]/I");
1081  theTree->Branch("nhPXF", m_nhPXF, "nhPXF[Ntracks]/I");
1082  theTree->Branch("nhTIB", m_nhTIB, "nhTIB[Ntracks]/I");
1083  theTree->Branch("nhTOB", m_nhTOB, "nhTOB[Ntracks]/I");
1084  theTree->Branch("nhTID", m_nhTID, "nhTID[Ntracks]/I");
1085  theTree->Branch("nhTEC", m_nhTEC, "nhTEC[Ntracks]/I");
1086  theTree->Branch("Pt", m_Pt, "Pt[Ntracks]/F");
1087  theTree->Branch("P", m_P, "P[Ntracks]/F");
1088  theTree->Branch("Eta", m_Eta, "Eta[Ntracks]/F");
1089  theTree->Branch("Phi", m_Phi, "Phi[Ntracks]/F");
1090  theTree->Branch("Chi2n", m_Chi2n, "Chi2n[Ntracks]/F");
1091  theTree->Branch("d0", m_d0, "d0[Ntracks]/F");
1092  theTree->Branch("dz", m_dz, "dz[Ntracks]/F");
1093 
1094  // book Alignable-wise ROOT Tree
1095 
1096  theFile2 = new TFile(outfile2.c_str(),"update");
1097  theFile2->cd();
1098 
1099  theTree2 = new TTree("T2","Alignablewise tree");
1100 
1101  theTree2->Branch("Nhit", &m2_Nhit, "Nhit/I");
1102  theTree2->Branch("Type", &m2_Type, "Type/I");
1103  theTree2->Branch("Layer", &m2_Layer, "Layer/I");
1104  theTree2->Branch("Xpos", &m2_Xpos, "Xpos/F");
1105  theTree2->Branch("Ypos", &m2_Ypos, "Ypos/F");
1106  theTree2->Branch("Zpos", &m2_Zpos, "Zpos/F");
1107  theTree2->Branch("Eta", &m2_Eta, "Eta/F");
1108  theTree2->Branch("Phi", &m2_Phi, "Phi/F");
1109  theTree2->Branch("Id", &m2_Id, "Id/i");
1110  theTree2->Branch("ObjId", &m2_ObjId, "ObjId/I");
1111 
1112  // book survey-wise ROOT Tree only if survey is enabled
1113  if (theLevels.size() > 0){
1114 
1115  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::bookRoot] Survey trees booked.";
1116  theFile3 = new TFile(ssurveyfile.c_str(),"update");
1117  theFile3->cd();
1118  theTree3 = new TTree(tname, "Survey Tree");
1119  theTree3->Branch("Id", &m3_Id, "Id/i");
1120  theTree3->Branch("ObjId", &m3_ObjId, "ObjId/I");
1121  theTree3->Branch("Par", &m3_par, "Par[6]/F");
1122  }
1123 
1124  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::bookRoot] Root trees booked.";
1125 
1126 }
1127 
1128 // ----------------------------------------------------------------------------
1129 // fill alignable-wise root tree
1130 
1132 {
1133  using std::setw;
1134  theFile2->cd();
1135 
1136  int naligned=0;
1137 
1138  //Retrieve tracker topology from geometry
1139  edm::ESHandle<TrackerTopology> tTopoHandle;
1140  iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
1141  const TrackerTopology* const tTopo = tTopoHandle.product();
1142 
1143  for (std::vector<Alignable*>::const_iterator it=theAlignables.begin();
1144  it!=theAlignables.end();
1145  ++it) {
1146  Alignable* ali = (*it);
1148 
1149  // consider only those parameters classified as 'valid'
1150  if (dap->isValid()) {
1151 
1152  // get number of hits from user variable
1153  HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(dap->userVariables());
1154  m2_Nhit = uservar->nhit;
1155 
1156  // get type/layer
1157  std::pair<int,int> tl = theAlignmentParameterStore->typeAndLayer(ali, tTopo);
1158  m2_Type = tl.first;
1159  m2_Layer = tl.second;
1160 
1161  // get identifier (as for IO)
1162  m2_Id = ali->id();
1163  m2_ObjId = ali->alignableObjectId();
1164 
1165  // get position
1166  GlobalPoint pos = ali->surface().position();
1167  m2_Xpos = pos.x();
1168  m2_Ypos = pos.y();
1169  m2_Zpos = pos.z();
1170  m2_Eta = pos.eta();
1171  m2_Phi = pos.phi();
1172 
1173  AlgebraicVector pars = dap->parameters();
1174 
1175  if (verbose) {
1176  edm::LogVerbatim("Alignment")
1177  << "------------------------------------------------------------------------\n"
1178  << " ALIGNABLE: " << setw(6) << naligned
1179  << '\n'
1180  << "hits: " << setw(4) << m2_Nhit
1181  << " type: " << setw(4) << m2_Type
1182  << " layer: " << setw(4) << m2_Layer
1183  << " id: " << setw(4) << m2_Id
1184  << " objId: " << setw(4) << m2_ObjId
1185  << '\n'
1186  << std::fixed << std::setprecision(5)
1187  << "x,y,z: "
1188  << setw(12) << m2_Xpos
1189  << setw(12) << m2_Ypos
1190  << setw(12) << m2_Zpos
1191  << " eta,phi: "
1192  << setw(12) << m2_Eta
1193  << setw(12) << m2_Phi
1194  << '\n'
1195  << "params: "
1196  << setw(12) << pars[0]
1197  << setw(12) << pars[1]
1198  << setw(12) << pars[2]
1199  << setw(12) << pars[3]
1200  << setw(12) << pars[4]
1201  << setw(12) << pars[5];
1202  }
1203 
1204  naligned++;
1205  theTree2->Fill();
1206  }
1207  }
1208 }
1209 
1210 // ----------------------------------------------------------------------------
1211 
1213 {
1214  // Alignment parameters
1216  // access user variables
1217  HIPUserVariables* uservar = dynamic_cast<HIPUserVariables*>(par->userVariables());
1218  int nhit = uservar->nhit;
1219  // The following variable is needed for the extended 1D/2D hit fix using
1220  // matrix shrinkage and expansion
1221  // int hitdim = uservar->hitdim;
1222 
1223  if (nhit < theMinimumNumberOfHits) {
1224  par->setValid(false);
1225  return false;
1226  }
1227 
1228  AlgebraicSymMatrix jtvj = uservar->jtvj;
1229  AlgebraicVector jtve = uservar->jtve;
1230 
1231  // Shrink input in case of 1D hits and 'v' selected
1232  // in alignment parameters
1233  // if (hitdim==1 && selector[1]==true) {
1234  // int iremove = 1;
1235  // if (selector[0]==false) iremove--;
1236 
1237  // AlgebraicSymMatrix tempjtvj(jtvj.num_row()-1);
1238  // int nr = 0, nc = 0;
1239  // for (int r=0;r<jtvj.num_row();r++) {
1240  // if (r==iremove) continue;
1241  // nc = 0;
1242  // for (int c=0;c<jtvj.num_col();c++) {
1243  // if (c==iremove) continue;
1244  // tempjtvj[nr][nc] = jtvj[r][c];
1245  // nc++;
1246  // }
1247  // nr++;
1248  // }
1249  // jtvj = tempjtvj;
1250 
1251  // AlgebraicVector tempjtve(jtve.num_row()-1);
1252  // nr = 0;
1253  // for (int r=0;r<jtve.num_row();r++) {
1254  // if (r==iremove) continue;
1255  // tempjtve[nr] = jtve[r];
1256  // nr++;
1257  // }
1258  // jtve = tempjtve;
1259  // }
1260 
1261  int ierr;
1262  AlgebraicSymMatrix jtvjinv = jtvj.inverse(ierr);
1263 
1264  if (ierr !=0) {
1265  edm::LogError("Alignment") << "Matrix inversion failed!";
1266  return false;
1267  }
1268 
1269  // these are the alignment corrections+covariance (for selected params)
1270  AlgebraicVector params = - (jtvjinv * jtve);
1271  AlgebraicSymMatrix cov = jtvjinv;
1272 
1273  edm::LogInfo("Alignment") << "parameters " << params;
1274 
1275  // errors of parameters
1276  int npar = params.num_row();
1277  AlgebraicVector paramerr(npar);
1278  AlgebraicVector relerr(npar);
1279  for (int i=0;i<npar;i++) {
1280  if (abs(cov[i][i])>0) paramerr[i] = sqrt(abs(cov[i][i]));
1281  else paramerr[i] = params[i];
1282  relerr[i] = abs(paramerr[i]/params[i]);
1283  if (relerr[i] >= theMaxRelParameterError) {
1284  params[i] = 0;
1285  paramerr[i]=0;
1286  }
1287  }
1288 
1289  // expand output in case of 1D hits and 'v' selected
1290  // in alignment parameters
1291  // if (hitdim==1 && selector[1]==true) {
1292  // int iremove = 1;
1293  // if (selector[0]==false) iremove--;
1294 
1295  // AlgebraicSymMatrix tempcov(cov.num_row()+1,0);
1296  // int nr = 0, nc = 0;
1297  // for (int r=0;r<cov.num_row();r++) {
1298  // if (r==iremove) nr++;
1299  // nc = 0;
1300  // for (int c=0;c<cov.num_col();c++) {
1301  // if (c==iremove) nc++;
1302  // tempcov[nr][nc] = cov[r][c];
1303  // nc++;
1304  // }
1305  // nr++;
1306  // }
1307  // cov = tempcov;
1308 
1309  // AlgebraicVector tempparams(params.num_row()+1,0);
1310  // nr = 0;
1311  // for (int r=0;r<params.num_row();r++) {
1312  // if (r==iremove) nr++;
1313  // tempparams[nr] = params[r];
1314  // nr++;
1315  // }
1316  // params = tempparams;
1317  // }
1318 
1319  // store alignment parameters
1320  AlignmentParameters* parnew = par->cloneFromSelected(params,cov);
1321  ali->setAlignmentParameters(parnew);
1322  parnew->setValid(true);
1323 
1324  return true;
1325 }
1326 
1327 //-----------------------------------------------------------------------------
1328 
1330 {
1331  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::collector] called for iteration "
1332  << theIteration << std::endl;
1333 
1334  HIPUserVariablesIORoot HIPIO;
1335 
1336  for (int ijob=1;ijob<=theCollectorNJobs;ijob++) {
1337 
1338  edm::LogWarning("Alignment") << "reading uservar for job " << ijob;
1339 
1340  std::stringstream ss;
1341  std::string str;
1342  ss << ijob;
1343  ss >> str;
1344  std::string uvfile = theCollectorPath+"/job"+str+"/IOUserVariables.root";
1345 
1346  std::vector<AlignmentUserVariables*> uvarvec =
1347  HIPIO.readHIPUserVariables(theAlignables, uvfile.c_str(),
1348  theIteration, ioerr);
1349 
1350  if (ioerr!=0) {
1351  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::collector] could not read user variable files for job "
1352  << ijob;
1353  continue;
1354  }
1355 
1356  // add
1357  std::vector<AlignmentUserVariables*> uvarvecadd;
1358  std::vector<AlignmentUserVariables*>::const_iterator iuvarnew=uvarvec.begin();
1359  for (std::vector<Alignable*>::const_iterator it=theAlignables.begin();
1360  it!=theAlignables.end();
1361  ++it) {
1362  Alignable* ali = *it;
1364 
1365  HIPUserVariables* uvarold = dynamic_cast<HIPUserVariables*>(ap->userVariables());
1366  HIPUserVariables* uvarnew = dynamic_cast<HIPUserVariables*>(*iuvarnew);
1367 
1368  HIPUserVariables* uvar = uvarold->clone();
1369  if (uvarnew!=0) {
1370  uvar->nhit = (uvarold->nhit)+(uvarnew->nhit);
1371  uvar->jtvj = (uvarold->jtvj)+(uvarnew->jtvj);
1372  uvar->jtve = (uvarold->jtve)+(uvarnew->jtve);
1373  uvar->alichi2 = (uvarold->alichi2)+(uvarnew->alichi2);
1374  uvar->alindof = (uvarold->alindof)+(uvarnew->alindof);
1375  delete uvarnew;
1376  }
1377 
1378  uvarvecadd.push_back(uvar);
1379  iuvarnew++;
1380  }
1381 
1383 
1384  //fill Eventwise Tree
1385  if (theFillTrackMonitoring) {
1386  uvfile = theCollectorPath+"/job"+str+"/HIPAlignmentEvents.root";
1387  edm::LogWarning("Alignment") << "Added to the tree "
1388  << fillEventwiseTree(uvfile.c_str(), theIteration, ioerr)
1389  << "tracks";
1390  }
1391 
1392  }//end loop on jobs
1393 }
1394 
1395 //------------------------------------------------------------------------------------
1396 
1397 int HIPAlignmentAlgorithm::fillEventwiseTree(const char* filename, int iter, int ierr)
1398 {
1399  int totntrk = 0;
1400  char treeName[64];
1401  snprintf(treeName, sizeof(treeName), "T1_%d", iter);
1402  //open the file "HIPAlignmentEvents.root" in the job directory
1403  TFile *jobfile = new TFile(filename, "READ");
1404  //grab the tree corresponding to this iteration
1405  TTree *jobtree = (TTree*)jobfile->Get(treeName);
1406  //address and read the variables
1407  static const int nmaxtrackperevent = 1000;
1408  int jobNtracks, jobNhitspertrack[nmaxtrackperevent], jobnhPXB[nmaxtrackperevent], jobnhPXF[nmaxtrackperevent],jobnhTIB[nmaxtrackperevent], jobnhTOB[nmaxtrackperevent],jobnhTID[nmaxtrackperevent], jobnhTEC[nmaxtrackperevent];
1409  float jobP[nmaxtrackperevent], jobPt[nmaxtrackperevent], jobEta[nmaxtrackperevent] , jobPhi[nmaxtrackperevent];
1410  float jobd0[nmaxtrackperevent], jobdz[nmaxtrackperevent] , jobChi2n[nmaxtrackperevent];
1411 
1412  jobtree->SetBranchAddress("Ntracks", &jobNtracks);
1413  jobtree->SetBranchAddress("Nhits", jobNhitspertrack);
1414  jobtree->SetBranchAddress("nhPXB", jobnhPXB);
1415  jobtree->SetBranchAddress("nhPXF", jobnhPXF);
1416  jobtree->SetBranchAddress("nhTIB", jobnhTIB);
1417  jobtree->SetBranchAddress("nhTOB", jobnhTOB);
1418  jobtree->SetBranchAddress("nhTID", jobnhTID);
1419  jobtree->SetBranchAddress("nhTEC", jobnhTEC);
1420  jobtree->SetBranchAddress("Pt", jobPt);
1421  jobtree->SetBranchAddress("P", jobP);
1422  jobtree->SetBranchAddress("d0", jobd0);
1423  jobtree->SetBranchAddress("dz", jobdz);
1424  jobtree->SetBranchAddress("Eta", jobEta);
1425  jobtree->SetBranchAddress("Phi", jobPhi);
1426  jobtree->SetBranchAddress("Chi2n", jobChi2n);
1427  int ievent = 0;
1428  for (ievent=0;ievent<jobtree->GetEntries();++ievent) {
1429  jobtree->GetEntry(ievent);
1430 
1431  //fill the collector tree with them
1432 
1433  // TO BE IMPLEMENTED: a prescale factor like in run()
1434  m_Ntracks = jobNtracks;
1435  int ntrk = 0;
1436  while (ntrk<m_Ntracks) {
1437  if (ntrk<MAXREC) {
1438  totntrk = ntrk+1;
1439  m_Nhits[ntrk] = jobNhitspertrack[ntrk];
1440  m_Pt[ntrk] = jobPt[ntrk];
1441  m_P[ntrk] = jobP[ntrk];
1442  m_nhPXB[ntrk] = jobnhPXB[ntrk];
1443  m_nhPXF[ntrk] = jobnhPXF[ntrk];
1444  m_nhTIB[ntrk] = jobnhTIB[ntrk];
1445  m_nhTOB[ntrk] = jobnhTOB[ntrk];
1446  m_nhTID[ntrk] = jobnhTID[ntrk];
1447  m_nhTEC[ntrk] = jobnhTEC[ntrk];
1448  m_Eta[ntrk] = jobEta[ntrk];
1449  m_Phi[ntrk] = jobPhi[ntrk];
1450  m_Chi2n[ntrk] = jobChi2n[ntrk];
1451  m_d0[ntrk] = jobd0[ntrk];
1452  m_dz[ntrk] = jobdz[ntrk];
1453  }//end if j<MAXREC
1454  else{
1455  edm::LogWarning("Alignment") << "[HIPAlignmentAlgorithm::fillEventwiseTree] Number of tracks in Eventwise tree exceeds MAXREC: "
1456  << m_Ntracks << " Skipping exceeding tracks.";
1457  ntrk = m_Ntracks+1;
1458  }
1459  ++ntrk;
1460  }//end while loop
1461  theTree->Fill();
1462  }//end loop on i - entries in the job tree
1463 
1464  //clean up
1465  delete jobtree;
1466  delete jobfile;
1467 
1468  return totntrk;
1469 }//end fillEventwiseTree
ClusterRef cluster() const
void attachUserVariables(const align::Alignables &alivec, const std::vector< AlignmentUserVariables * > &uvarvec, int &ierr)
Attach User Variables to given alignables.
double p() const
momentum vector magnitude
Definition: TrackBase.h:602
type
Definition: HCALResponse.h:21
T getParameter(std::string const &) const
align::ID id() const
Return the ID of Alignable, i.e. DetId of &#39;first&#39; component GeomDet(Unit).
Definition: Alignable.h:180
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
float xx() const
Definition: LocalError.h:24
tuple cfg
Definition: looper.py:259
ConstRecHitPointer const & recHit() const
double d0() const
dxy parameter in perigee convention (d0 = -dxy)
Definition: TrackBase.h:584
AlignableDetOrUnitPtr alignableFromGeomDet(const GeomDet *geomDet)
Returns AlignableDetOrUnitPtr corresponding to given GeomDet.
AlignmentParameterStore * theAlignmentParameterStore
std::pair< int, int > typeAndLayer(const Alignable *ali, const TrackerTopology *tTopo) const
Obtain type and layer from Alignable.
bool calcParameters(Alignable *ali)
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:548
TSOS combine(const TSOS &pTsos1, const TSOS &pTsos2) const
int fillEventwiseTree(const char *filename, int iter, int ierr)
void setAlignmentPositionError(const align::Alignables &alivec, double valshift, double valrot)
Set Alignment position error.
static char const * tname
Definition: GTSchema.h:13
void writeAlignableOriginalPositions(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
write Alignable original (before misalignment) absolute positions
assert(m_qm.get())
align::StructureType m2_ObjId
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:632
void writeIterationFile(std::string filename, int iter)
void applyParameters(void)
Obsolete: Use AlignableNavigator::alignableDetFromDetId and alignableFromAlignableDet.
AlgebraicVector jtve
Alignable * alignableFromAlignableDet(const AlignableDetOrUnitPtr &adet) const
Get relevant Alignable from AlignableDet.
HIPUserVariables * clone(void) const
void applyAlignableAbsolutePositions(const align::Alignables &alis, const AlignablePositions &newpos, int &ierr)
apply absolute positions to alignables
LocalError positionError() const
int numberOfValidStripTOBHits() const
Definition: HitPattern.h:846
AlignmentParameters * alignmentParameters() const
Get the AlignmentParameters.
Definition: Alignable.h:57
const ConstTrajTrackPairCollection & trajTrackPairs() const
define event information passed to algorithms
void writeHIPUserVariables(const Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
static align::StructureType stringToId(const char *)
AlgebraicSymMatrix jtvj
std::vector< std::pair< std::vector< Alignable * >, std::vector< double > > > theAPEParameters
const AlgebraicVector & parameters(void) const
Get alignment parameters.
DataContainer const & measurements() const
Definition: Trajectory.h:203
AlignmentUserVariables * userVariables(void) const
Get pointer to user variables.
AlignablePositions readAlignableAbsolutePositions(const align::Alignables &alivec, const char *filename, int iter, int &ierr)
read Alignable current absolute positions
float xy() const
Definition: LocalError.h:25
void clear()
remove all selected Alignables and geometrical restrictions
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:638
double calcAPE(double *par, int iter, double function)
void startNewLoop(void)
Called at start of new loop.
int numberOfValidPixelBarrelHits() const
Definition: HitPattern.h:821
std::vector< AlignmentUserVariables * > readHIPUserVariables(const Alignables &alivec, const char *filename, int iter, int &ierr)
void setAlignmentParameters(AlignmentParameters *dap)
Set the AlignmentParameters.
Definition: Alignable.cc:81
CLHEP::HepMatrix AlgebraicMatrix
HIPAlignmentAlgorithm(const edm::ParameterSet &cfg)
Constructor.
tuple iov
Definition: o2o.py:307
float yy() const
Definition: LocalError.h:26
void setValid(bool v)
Set validity flag.
T sqrt(T t)
Definition: SSEVec.h:48
double pt() const
track transverse momentum
Definition: TrackBase.h:608
T z() const
Definition: PV3DBase.h:64
tuple result
Definition: query.py:137
void initialize(const edm::EventSetup &setup, AlignableTracker *tracker, AlignableMuon *muon, AlignableExtras *extras, AlignmentParameterStore *store)
Call at beginning of job.
virtual StructureType alignableObjectId() const =0
Return the alignable type identifier.
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool processHit2D(const AlignableDetOrUnitPtr &alidet, const Alignable *ali, const TrajectoryStateOnSurface &tsos, const TrackingRecHit *hit)
AlignableNavigator * theAlignableDetAccessor
int numberOfValidStripTIDHits() const
Definition: HitPattern.h:841
ClusterRef cluster() const
void fillRoot(const edm::EventSetup &setup)
int readIterationFile(std::string filename)
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:807
int numberOfValidStripTECHits() const
Definition: HitPattern.h:851
bool processHit1D(const AlignableDetOrUnitPtr &alidet, const Alignable *ali, const TrajectoryStateOnSurface &tsos, const TrackingRecHit *hit)
const LocalTrajectoryError & localError() const
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
string key
FastSim: produces sample of signal events, overlayed with premixed minbias events.
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Definition: Alignable.h:126
TrajectoryStateOnSurface const & forwardPredictedState() const
Access to forward predicted state (from fitter or builder)
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:596
tuple out
Definition: dbtoconf.py:99
const align::Alignables & selectedAlignables() const
vector of alignables selected so far
int numSelected(void) const
Get number of selected parameters.
virtual TrackingRecHit const * hit() const
CLHEP::HepVector AlgebraicVector
void setUserVariables(AlignmentUserVariables *auv)
Set pointer to user variables.
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:437
tuple tracks
Definition: testEve_cfg.py:39
virtual LocalError localPositionError() const =0
ClusterRef cluster() const
Definition: SiPixelRecHit.h:49
const T & get() const
Definition: EventSetup.h:56
void writeAlignmentParameters(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
write AlignmentParameters
CompositeAlignmentParameters selectParameters(const std::vector< AlignableDet * > &alignabledets) const
T const * product() const
Definition: ESHandle.h:86
AlgebraicSymMatrix inverseCovariance() const
Get inverse of survey covariance wrt given structure type in constructor.
bool isValid() const
int numberOfValidStripTIBHits() const
Definition: HitPattern.h:836
std::vector< edm::ParameterSet > theAPEParameterSet
align::StructureType m3_ObjId
int numberOfValidPixelEndcapHits() const
Definition: HitPattern.h:826
T eta() const
Definition: PV3DBase.h:76
virtual void terminate()
Called at end of job (must be implemented in derived class)
AlgebraicVector sensorResidual() const
std::vector< AlignableDetOrUnitPtr > alignablesFromHits(const std::vector< const TransientTrackingRecHit * > &hitvec)
Returns vector AlignableDetOrUnitPtr for given vector of Hits.
const AliClusterValueMap * clusterValueMap() const
bool isValid(void) const
Get validity flag.
tuple filename
Definition: lut2db_cfg.py:20
CLHEP::HepSymMatrix AlgebraicSymMatrix
unsigned int addSelections(const edm::ParameterSet &pSet)
std::vector< AlignableAbsData > AlignablePositions
Definition: AlignableData.h:51
void writeAlignableAbsolutePositions(const align::Alignables &alivec, const char *filename, int iter, bool validCheck, int &ierr)
write Alignable current absolute positions
DetId geographicalId() const
virtual AlignmentParameters * cloneFromSelected(const AlgebraicVector &par, const AlgebraicSymMatrix &cov) const =0
bool detAndSubdetInMap(const DetId &detid) const
Given a DetId, returns true if DetIds with this detector and subdetector id are in the map (not neces...
Constructor of the full muon geometry.
Definition: AlignableMuon.h:36
T x() const
Definition: PV3DBase.h:62
ValidityInterval const & validityInterval() const
virtual LocalPoint localPosition() const =0
const PositionType & position() const
std::vector< align::StructureType > theLevels
const Time_t MAX_VAL(std::numeric_limits< Time_t >::max())
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
std::vector< Alignable * > theAlignables
bool valid() const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
const align::Alignables & alignables(void) const
get all alignables
std::vector< ConstTrajTrackPair > ConstTrajTrackPairCollection
std::string className(const T &t)
Definition: ClassName.h:30
TrajectoryStateOnSurface const & backwardPredictedState() const
Access to backward predicted state (from smoother)
Our base class.
Definition: SiPixelRecHit.h:23
void run(const edm::EventSetup &setup, const EventInfo &eventInfo)
Run the algorithm.
virtual AlgebraicMatrix selectedDerivatives(const TrajectoryStateOnSurface &tsos, const AlignableDetOrUnitPtr &alidet) const