CMS 3D CMS Logo

Fit.cc
Go to the documentation of this file.
2 // COCOA class implementation file
3 //Id: Fit.cc
4 //CAT: Fit
5 //
6 // History: v1.0
7 // Pedro Arce
8 
9 #include <cstdlib>
10 #include <iomanip>
11 #include <cmath> // among others include also floating-point std::abs functions
12 #include <ctime>
13 #include <set>
14 
18 
28 #ifdef COCOA_VIS
29 #include "Alignment/CocoaVisMgr/interface/ALIVRMLMgr.h"
30 #include "Alignment/IgCocoaFileWriter/interface/IgCocoaFileMgr.h"
31 #endif
38 
39 Fit* Fit::instance = nullptr;
40 
45 //op ALIMatrix* Fit::VaMatrix;
47 //op ALIMatrix* Fit::PDMatrix;
48 //-ALIMatrix* Fit::VyMatrix;
50 //ALIMatrix* Fit::fMatrix;
51 
54 //op ALIMatrix* Fit::thePropagationMatrix;
55 
63 
65 
66 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
67 //@@ Gets the only instance of Model
68 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
70  if (!instance) {
71  instance = new Fit;
72  ALIdouble go;
74 
75  gomgr->getGlobalOptionValue("maxDeviDerivative", go);
77  if (ALIUtils::debug >= 3)
78  std::cout << " Fit::maximum_deviation_derivative " << ALIUtils::getMaximumDeviationDerivative() << std::endl;
79 
80  gomgr->getGlobalOptionValue("maxNoFitIterations", go);
81  MaxNoFitIterations = int(go);
82 
83  gomgr->getGlobalOptionValue("fitQualityCut", go);
84  theFitQualityCut = go;
85  if (ALIUtils::debug >= 3)
86  std::cout << " theFitQualityCut " << theFitQualityCut << std::endl;
87 
88  gomgr->getGlobalOptionValue("RelativeFitQualityCut", go);
90  if (ALIUtils::debug >= 3)
91  std::cout << " theRelativeFitQualityCut " << theRelativeFitQualityCut << std::endl;
92 
93  gomgr->getGlobalOptionValue("minDaFactor", go);
94  theMinDaFactor = go;
95  }
96 
97  return *instance;
98 }
99 
100 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
101 //@@ startFit: steering method to make the fit
102 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
104  // Model::setCocoaStatus( COCOA_InitFit );
106  if (GlobalOptionMgr::getInstance()->GlobalOptions()["rootResults"] > 0) {
107  NTmgr->BookNtuple();
108  }
109 
111 
113 
115  for (;;) {
116  bool bend = fitNextEvent(nEvent);
117  if (gomgr->GlobalOptions()["writeDBOptAlign"] > 0 || gomgr->GlobalOptions()["writeDBAlign"] > 0) {
119  }
120 
121  if (!bend) {
122  if (ALIUtils::debug >= 1)
123  std::cout << "@@@ Fit::startFit ended n events = " << nEvent << std::endl;
124  break;
125  }
126 
127  //- if ( ALIUtils::debug >= 0) std::cout << " FIT STATUS " << Model::printCocoaStatus( Model::getCocoaStatus() ) << std::endl;
128 
129  nEvent++;
130  }
131 
132  //---------- Program ended, fill histograms of fitted entries
133  if (gomgr->GlobalOptions()["histograms"] > 0) {
135  FEmgr->MakeHistos();
136  }
137 
138  if (GlobalOptionMgr::getInstance()->GlobalOptions()["rootResults"] > 0) {
139  NTmgr->WriteNtuple();
140  }
141 }
142 
143 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
145  if (Model::getFittedEntriesReader() != nullptr)
147 
148  //----- Reset coordinates to those read at the start
149  std::vector<OpticalObject*>::iterator voite;
150  for (voite = Model::OptOList().begin(); voite != Model::OptOList().end(); ++voite) {
151  (*voite)->resetOriginalOriginalCoordinates();
152  }
153 
154  //----- Reset entries displacements to 0.
155  std::vector<Entry*>::iterator veite;
156  for (veite = Model::EntryList().begin(); veite != Model::EntryList().end(); ++veite) {
157  (*veite)->resetValueDisplacementByFitting();
158  }
159 
160  ALIbool lastEvent = false;
161 
162  //- DeviationsFromFileSensor2D::setApply( 1 );
163 
164  //m ALIbool moreDataSets = Model::readMeasurementsFromFile( Measurement::only1Date, Measurement::only1Time );
165 
166  //----- Check if there are more data sets
167  ALIbool moreDataSets = true;
168  if (CocoaDaqReader::GetDaqReader() != nullptr)
169  moreDataSets = CocoaDaqReader::GetDaqReader()->ReadNextEvent();
170 
171  if (ALIUtils::debug >= 5)
172  std::cout << CocoaDaqReader::GetDaqReader() << "$$$$$$$$$$$$$$$ More Data Sets to be processed: " << moreDataSets
173  << std::endl;
174 
175  if (moreDataSets) {
176  if (ALIUtils::debug >= 2)
177  std::cout << std::endl << "@@@@@@@@@@@@@@@@@@ Starting data set fit : " << nEvent << std::endl;
178 
179  //----- Count entries to be fitted, and set their order in theFitPos
181 
182  //----- Dump dimensions of output in 'report.out' file
184  fileout << std::endl << "@@@@@@@ NEW MEASUREMENT SET " << nEvent << std::endl;
185  if (ALIUtils::report >= 1)
186  ALIUtils::dumpDimensions(fileout);
187 
188  //----- reset Number of iterations of non linear fit
189  theNoFitIterations = 0;
190 
192  ALIdouble dumpMat;
193  gomgr->getGlobalOptionValue("save_matrices", dumpMat);
194 
195  //----- Fit parameters
196  double daFactor = 1.;
198  for (;;) {
199  if (ALIUtils::debug >= 2) {
200  std::cout << std::endl << "Fit iteration " << theNoFitIterations << " ..." << std::endl;
201  }
202 
203  //---------- Calculate the original simulated values of each Measurement (when all entries have their read in values)
204  calculateSimulatedMeasurementsWithOriginalValues(); //?? original changed atfer each iteration
205 
206  FitQuality fq = fitParameters(daFactor);
207  if (dumpMat > 1)
208  dumpMatrices();
209 
210  //- evaluateFitQuality( fq, daFactor );
211 
212  if (ALIUtils::debug >= 2) {
213  std::cout << std::endl << "@@@@ Check fit quality for iteration " << theNoFitIterations << std::endl;
214  }
215 
216  //----- Check if new iteration must be done
217  if (fq == FQsmallDistanceToMinimum) {
218  if (ALIUtils::debug >= 2)
219  std::cout << std::endl << "@@@@ Fit quality: distance SMALLER than mininum " << std::endl;
221  if (ALIUtils::report >= 1)
223  //--- Print entries in all ancestor frames
224  ALIdouble go;
225  gomgr->getGlobalOptionValue("dumpInAllFrames", go);
226  if (go >= 1)
228 
229  break; // No more iterations
230  } else if (fq == FQbigDistanceToMinimum) {
231  if (ALIUtils::debug >= 2)
232  std::cout << std::endl << "@@@@ Fit quality: distance BIGGER than mininum " << std::endl;
234  if (ALIUtils::report >= 1)
236 
237  //----- Next iteration (if not too many already)
239  daFactor = 1.;
240 
241  //----- Too many iterations: end event here
243  if (ALIUtils::debug >= 1)
244  std::cerr << "!!!! WARNING: Too many iterations " << theNoFitIterations << " and fit DOES NOT CONVERGE "
245  << std::endl;
246 
247  if (ALIUtils::report >= 2) {
249  fileout << "!!!! WARNING: Too many iterations " << theNoFitIterations << " and fit DOES NOT CONVERGE "
250  << std::endl;
251  }
252  // Model::setCocoaStatus( COCOA_FitCannotImprove );
253  break; // No more iterations
254  }
255 
256  } else if (fq == FQchiSquareWorsened) {
257  if (ALIUtils::debug >= 1) {
258  //----- Recalculate fit quality with decreasing values of Da
259  std::cerr << "!! WARNING: fit quality has worsened, Recalculate fit quality with decreasing values of Da "
260  << std::endl;
261  std::cout << " quality daFactor= " << daFactor << " minimum= " << theMinDaFactor << std::endl;
262  }
263  daFactor *= 0.5;
264  if (daFactor > theMinDaFactor) {
266 
267  if (ALIUtils::report >= 2) {
269  fileout << " Redoing iteration with Da factor " << daFactor << std::endl;
270  }
271  } else {
272  daFactor *= 2.;
273  std::cerr << " !!!ERROR: not possible to get good fit quality even multiplying Da by " << daFactor
274  << std::endl;
275  if (ALIUtils::report >= 2) {
277  fileout << " !!!ERROR: not possible to get good fit quality even multiplying Da by " << daFactor
278  << std::endl;
279  }
280  // Model::setCocoaStatus( COCOA_FitCannotImprove );
281  //- std::cout << "fdsaf FIT STATUS " << Model::printCocoaStatus( Model::getCocoaStatus() ) << std::endl;
282  break; // No more iterations
283  }
284  }
286  }
287 
288  //----- Iteration is finished: dump fitted entries
289  if (ALIUtils::debug >= 1)
291  if (gomgr->GlobalOptions()["histograms"] > 0) {
293  }
294 
295  if (GlobalOptionMgr::getInstance()->GlobalOptions()["rootResults"] > 0) {
297  ntupleMgr->InitNtuple();
298  ntupleMgr->FillChi2();
299  ntupleMgr->FillOptObjects(AtWAMatrix);
300  ntupleMgr->FillMeasurements();
301  ntupleMgr->FillFitParameters(AtWAMatrix);
302  ntupleMgr->FillNtupleTree();
303  }
304 
305  //- only if not stopped in worsening quality state if(ALIUtils::report >= 0) dumpFittedValues( ALIFileOut::getInstance( Model::ReportFName() ));
306 
307  /*- std::vector< OpticalObject* >::iterator voite;
308  for( voite = Model::OptOList().begin(); voite != Model::OptOList().end(); voite++ ) {
309  //-?? (*voite)->resetOriginalOriginalCoordinates();
310  }*/
311 
312  //---- If no measurement file, break after looping once
313  //- std::cout << " Measurement::measurementsFileName() " << Measurement::measurementsFileName() << " Measurement::measurementsFileName()" <<std::endl;
314  if (CocoaDaqReader::GetDaqReader() == nullptr) {
315  //m if( Measurement::measurementsFileName() == "" ) {
316  if (ALIUtils::debug >= 1)
317  std::cout << std::endl << "@@@@@@@@@@@@@@@@@@ Fit has ended : only one measurement " << nEvent << std::endl;
318  lastEvent = true;
319  return !lastEvent;
320  }
321 
322  //- std::cout << " Measurement::only1" << Measurement::only1 << std::endl;
323  if (Measurement::only1) {
324  if (ALIUtils::debug >= 1)
325  std::cout << std::endl << "@@@@@@@@@@@@@@@@@@ Fit has ended : 'Measurement::only1' is set" << std::endl;
326 
327  lastEvent = true;
328  return !lastEvent;
329  }
330 
331  if (GlobalOptionMgr::getInstance()->GlobalOptions()["maxEvents"] <= nEvent) {
332  if (ALIUtils::debug >= 1)
333  std::cout << std::endl
334  << "@@@@@@@@@@@@@@@@@@ Fit has ended : 'Number of events exhausted " << nEvent << std::endl;
335 
336  lastEvent = true;
337  return !lastEvent;
338  }
339 
340  } else {
341  lastEvent = true;
342  if (ALIUtils::debug >= 1)
343  std::cout << std::endl << "@@@@@@@@@@@@@@@@@@ Fit has ended : ??no more data sets' " << nEvent << std::endl;
344  return !lastEvent;
345  }
346 
347  if (ALIUtils::debug >= 1)
348  std::cout << std::endl << "@@@@@@@@@@@@@@@@@@ Fit has ended : " << nEvent << std::endl;
349 
350  return !lastEvent;
351 }
352 
353 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
355 #ifdef COCOA_VIS
356  if (gomgr->GlobalOptions()["VisOnly"] == 1) {
357  calculateSimulatedMeasurementsWithOriginalValues(); //?? original changed atfer each iteration
358  }
359 
361  if (gomgr->GlobalOptions()["VisWriteVRML"] > 0) {
363  ALIVRMLMgr::getInstance().writeFile();
364  }
365  if (gomgr->GlobalOptions()["VisWriteIguana"] > 0) {
367  IgCocoaFileMgr::getInstance().writeFile();
368  }
369 
370  if (gomgr->GlobalOptions()["VisOnly"] == 1) {
371  if (ALIUtils::debug >= 1)
372  std::cout << " Visualiation file(s) succesfully written. Ending.... " << std::endl;
373  exit(1);
374  }
375 #endif
376 }
377 
378 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
379 //@@ Count how many entries are going to be fitted (have quality >= theMinimumEntryQuality)
380 //@@ Set for this entries the value of theFitPos
381 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
383  std::vector<Entry*>::const_iterator vecite;
384 
386  theMinimumEntryQuality = int(gomgr->GlobalOptions()[ALIstring("calcul_type")]) + 1;
387  if (ALIUtils::debug >= 3)
388  std::cout << "@@@ Fit::setFittableEntries: total Entry List size= " << Model::EntryList().size() << std::endl;
389 
390  int No_entry_to_fit = 0;
391  for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
392  // Number the parameters that are going to be fitted
393  if ((*vecite)->quality() >= theMinimumEntryQuality) {
394  (*vecite)->setFitPos(No_entry_to_fit);
395  if (ALIUtils::debug >= 4)
396  std::cout << " Entry To Fit= " << No_entry_to_fit << " " << (*vecite)->OptOCurrent()->name() << " "
397  << (*vecite)->name() << " with quality= " << (*vecite)->quality() << std::endl;
398  No_entry_to_fit++;
399  }
400  }
401 }
402 
403 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
404 //@@ Main method in class Fit
405 //@@ fitParameters: get the parameters through the chi square fit
406 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
407 FitQuality Fit::fitParameters(const double daFactor) {
408  if (ALIUtils::debug >= 3)
409  std::cout << "@@@ Fit::fitParameters: Fit quality daFactor " << daFactor << std::endl;
410 
411  redoMatrices();
412 
413  //---- Get chi2 of first iteration
416 
418  if (gomgr->GlobalOptions()[ALIstring("stopAfter1stIteration")] == 1) {
419  std::cout << "@!! STOPPED by user after 1st iteration " << std::endl;
420  exit(1);
421  }
422  }
423 
424  /* //---------- Open output file
425  if( ALIUtils::report >= 1 ) {
426  ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
427  //t fileout << " REPORT.OUT " << std::endl;
428  //t ALIUtils::dumpDimensions( fileout );
429  fileout << std::endl << "Fit iteration " << theNoFitIterations << " ..." << std::endl;
430  }*/
431 
432  //- std::cout << "2 FIT STATUS " << Model::printCocoaStatus( Model::getCocoaStatus() ) << std::endl;
433 
434  if (ALIUtils::debug >= 10) {
435  std::cout << std::endl << " End fitParameters " << theNoFitIterations << " ..." << std::endl;
436  }
437 
438  return getFitQuality();
439 }
440 
441 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
443  if (ALIUtils::debug >= 3)
444  std::cout << "@@@ Fit::redoMatrices" << std::endl;
445 
446  deleteMatrices();
447 
449 
450  PropagateErrors();
451 }
452 
453 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
454 //@@ Propagate the Errors from the entries to the measurements
455 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
456 //cocoaStatus
458  if (ALIUtils::debug >= 3)
459  std::cout << "@@@ Fit::PropagateErrors" << std::endl;
460 
461  //----- Create empty matrices of appropiate size
462  CreateMatrices();
463 
464  //---- count running time
465  time_t now;
466  now = clock();
467  if (ALIUtils::debug >= 2)
468  std::cout << "TIME:CREATE_MAT : " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
470 
471  //----- Fill the A, W & y matrices with the measurements
473 
474  //---- count running time
475  now = clock();
476  if (ALIUtils::debug >= 2)
477  std::cout << "TIME:MAT_MEAS_FILLED: " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
479 
480  //----- Fill the A, W & y matrices with the calibrated parameters
482  if (gomgr->GlobalOptions()[ALIstring("calcul_type")] == 0) {
484 
485  //---- count running time
486  now = clock();
487  if (ALIUtils::debug >= 0)
488  std::cout << "TIME:MAT_CAL_FILLED : " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
490  }
491 
492  //----- Put by hand some correlations if known previously
494 
495  if (ALIUtils::debug >= 3)
496  WMatrix->Dump("WMatrix before inverse");
497 
498  //----- Check first that matrix can be inverted
499  if (m_norm1(WMatrix->MatNonConst()) == 0) {
500  // Model::setCocoaStatus( COCOA_FitMatrixNonInversable );
501  return; // Model::getCocoaStatus();
502  } else {
503  WMatrix->inverse();
504  }
505 
506  if (ALIUtils::debug >= 3)
507  AMatrix->Dump("AMatrix");
508  if (ALIUtils::debug >= 3)
509  WMatrix->Dump("WMatrix");
510  if (ALIUtils::debug >= 3)
511  yfMatrix->Dump("yfMatrix");
512 
513  if (gomgr->GlobalOptions()["onlyDeriv"] >= 1) {
514  std::cout << "ENDING after derivatives are calculated ('onlyDeriv' option set)" << std::endl;
515  exit(1);
516  }
517 
519 
520  now = clock();
521  if (ALIUtils::debug >= 0)
522  std::cout << "TIME:MAT_MULTIPLIED : " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
524 
525  if (ALIUtils::getFirstTime() == 1)
526  ALIUtils::setFirstTime(false);
527 }
528 
529 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
530 //@@ Calculate the simulated value of each Measurement propagating the LightRay when all the entries have their original values
531 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
533  // if( ALIUtils::debug >= 4) OpticalObjectMgr::getInstance()->dumpOptOs();
534 
535  //---------- Set DeviationsFromFileSensor2D::apply true
537 
538  if (ALIUtils::debug >= 3)
539  std::cout << "@@@ Fit::calculateSimulatedMeasurementsWithOriginalValues" << std::endl;
540  //---------- Loop Measurements
541  std::vector<Measurement*>::const_iterator vmcite;
542  for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
543  //----- Calculate Simulated Value Original
544  (*vmcite)->calculateOriginalSimulatedValue();
545  }
546 
547  //---------- Set DeviationsFromFileSensor2D::apply false
548  // It cannot be applied when calculating derivatives, because after a displacement the laser could hit another square in matrix and then cause a big step in the derivative
550 }
551 
552 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
554  //delete matrices created in previous iteration
555  delete DaMatrix;
556  delete AMatrix;
557  delete WMatrix;
558  delete yfMatrix;
559  //op delete fMatrix;
560  delete AtMatrix;
561  delete AtWAMatrix;
562  //op delete VaMatrix;
563  //- delete VyMatrix;
564  //op delete PDMatrix;
565 }
566 
567 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
568 //@@ Calculate the NoLines & NoColumns and create matrices
569 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
571  //---------- Count number of measurements
572  ALIint NoMeas = 0;
573  std::vector<Measurement*>::const_iterator vmcite;
574  for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
575  NoMeas += (*vmcite)->dim();
576  }
577  if (ALIUtils::debug >= 9)
578  std::cout << "NOMEAS" << NoMeas << std::endl;
579 
580  //-------- Count number of 'cal'ibrated parameters
581  ALIint nEnt_cal = 0;
582  ALIint noent = 0;
583  //- std::cout << Model::EntryList().size() << std::endl;
585  if (gomgr->GlobalOptions()["calcul_type"] == 0) { // fit also 'cal' parameters
586  //- if( ALIUtils::debug >= 9) std::cout << "NOENTCALLL " << nEnt_cal << std::endl;
587  if (ALIUtils::debug >= 5)
588  std::cout << " Count number of 'cal'ibrated parameters " << std::endl;
589  std::vector<Entry*>::iterator veite;
590  for (veite = Model::EntryList().begin(); veite != Model::EntryList().end(); ++veite) {
591  if ((*veite)->quality() == 1)
592  nEnt_cal++;
593  noent++;
594  if (ALIUtils::debug >= 6) {
595  std::cout << (*veite)->quality() << " " << (*veite)->OptOCurrent()->name() << " " << (*veite)->name()
596  << " # ENT CAL " << nEnt_cal << " # ENT " << noent << std::endl;
597  }
598  }
599  }
600 
601  //---------- Count number parameters to be fitted ('cal' + 'unk')
602  ALIint NoParamFit = 0;
603  std::vector<Entry*>::const_iterator vecite;
604  for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
605  if ((*vecite)->quality() >= theMinimumEntryQuality) {
606  NoParamFit++;
607  if (ALIUtils::debug >= 99)
608  std::cout << (*vecite)->quality() << (*vecite)->OptOCurrent()->name() << (*vecite)->name() << "NoParamFit"
609  << NoParamFit << std::endl;
610  // break;
611  }
612  }
613 
614  //---------- Create Matrices
615  ALIint NoLinesA = NoMeas + nEnt_cal;
616  ALIint NoColumnsA = NoParamFit;
617  AMatrix = new ALIMatrix(NoLinesA, NoColumnsA);
618 
619  ALIint NoLinesW = NoLinesA;
620  ALIint NoColumnsW = NoLinesA;
621  WMatrix = new ALIMatrix(NoLinesW, NoColumnsW);
622 
623  ALIint NoLinesY = NoLinesA;
624  //op yMatrix = new ALIMatrix( NoLinesY, 1 );
625  yfMatrix = new ALIMatrix(NoLinesY, 1);
626 
627  //op fMatrix = new ALIMatrix( NoLinesY, 1 );
628 
629  if (ALIUtils::debug >= 4)
630  std::cout << "CreateMatrices: NoLinesA = " << NoLinesA << " NoColumnsA = " << NoColumnsA << std::endl;
631 }
632 
633 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
634 //@@ Loop Measurements:
635 //@@ Fill Matrix A with derivatives respect to affecting entries
636 //@@ Fill Matrix W, y & f with values and sigmas of measurement coordinate
637 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
639  if (ALIUtils::debug >= 3)
640  std::cout << "@@@ Fit::FillMatricesWithMeasurements" << std::endl;
641 
642  int Aline = 0;
643 
644  //---------- Loop Measurements
645  std::vector<Measurement*>::const_iterator vmcite;
646  std::vector<Entry*>::const_iterator vecite;
647  for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
648  if (ALIUtils::debug >= 5)
649  std::cout << "FillMatricesWithMeasurements: measurement " << (*vmcite)->name() << " # entries affecting "
650  << (*vmcite)->affectingEntryList().size() << std::endl;
651 
652  //-------- Array of derivatives with respect to entries
653  ALIint measdim = (*vmcite)->dim();
654  std::vector<ALIdouble> derivRE;
655  // derivRE = new ALIdouble[measdim];
656 
657  //-------- Fill matrix A:
658  //------ Loop only Entries affecting this Measurement
659  //-std::cout << "number affecting entries: " << (*vmcite)->affectingEntryList().size() << std::endl;
660  for (vecite = (*vmcite)->affectingEntryList().begin(); vecite != (*vmcite)->affectingEntryList().end(); ++vecite) {
661  //-------- If good quality, get derivative of measurement with respect to this Entry
662  if ((*vecite)->quality() >= theMinimumEntryQuality) {
663  if (ALIUtils::debug >= 4) {
664  // std::cout << "FillMatricesWithMeasurements: filling element ( " << Aline << " - " << Aline+measdim-1 << " , " << (*vecite)->fitPos() << std::endl;
665  std::cout << "entry affecting: " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << std::endl;
666  }
667  derivRE = (*vmcite)->DerivativeRespectEntry(*vecite);
668  //---------- Fill matrix A with derivatives
669  for (ALIuint jj = 0; jj < ALIuint(measdim); jj++) {
670  AMatrix->AddData(Aline + jj, (*vecite)->fitPos(), derivRE[jj]);
671  if (ALIUtils::debug >= 5)
672  std::cout << "FillMatricesWithMeasurements: AMATRIX (" << Aline + jj << "," << (*vecite)->fitPos() << " = "
673  << derivRE[jj] << std::endl;
674  //---------- Reset Measurement simulated_value
675  (*vmcite)->setValueSimulated(jj, (*vmcite)->valueSimulated_orig(jj));
676  }
677  }
678  }
679  // delete[] derivRE;
680 
681  //---------- Fill matrices W, y and f:
682  //------ Loop Measurement coordinates
683  for (ALIuint jj = 0; jj < ALIuint((*vmcite)->dim()); jj++) {
684  ALIdouble sigma = (*vmcite)->sigma()[jj];
685  if (sigma == 0.) {
686  std::cerr << "EXITING: Measurement number " << vmcite - Model::MeasurementList().begin() << "has 0 error!!"
687  << std::endl;
688  } else {
689  //----- Fill W Matrix with inverse of sigma squared
690  // multiply error by cameraScaleFactor
691  ALIdouble sigmanew = sigma * Measurement::cameraScaleFactor;
692  // std::cout << Aline+jj << " WMATRIX FILLING " << sigmanew << " * " << Measurement::cameraScaleFactor << std::endl;
693  WMatrix->AddData(Aline + jj, Aline + jj, (sigmanew * sigmanew));
694  }
695  //op //----- Fill Matrices y with measurement value
696  //op yMatrix->AddData( Aline+jj, 0, (*vmcite)->value()[jj] );
697  //op //----- Fill f Matrix with simulated_value
698  //op fMatrix->AddData( Aline+jj, 0, (*vmcite)->valueSimulated_orig(jj) );
699  //----- Fill Matrix y - f with measurement value - simulated value
700  yfMatrix->AddData(Aline + jj, 0, (*vmcite)->value()[jj] - (*vmcite)->valueSimulated_orig(jj));
701  // std::cout << " yfMatrix FILLING " << Aline+jj << " + " << (*vmcite)->value()[jj] - (*vmcite)->valueSimulated_orig(jj) << " meas " << (*vmcite)->name() << " val " << (*vmcite)->value()[jj] << " simu val " << (*vmcite)->valueSimulated_orig(jj) << std::endl;
702  }
703  if (ALIUtils::debug >= 99)
704  std::cout << "change line" << Aline << std::endl;
705  Aline += measdim;
706  if (ALIUtils::debug >= 99)
707  std::cout << "change line" << Aline << std::endl;
708  }
709 
710  if (ALIUtils::debug >= 4)
711  AMatrix->Dump("Matrix A with meas");
712  if (ALIUtils::debug >= 4)
713  WMatrix->Dump("Matrix W with meas");
714  if (ALIUtils::debug >= 4)
715  yfMatrix->Dump("Matrix y with meas");
716 }
717 
718 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
719 //@@ Loop Measurements:
720 //@@ Fill Matrix A with derivatives respect to affecting entries
721 //@@ Fill Matrix W, y & f with values and sigmas of measurement coordinate
722 //@@
723 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
725  if (ALIUtils::debug >= 3)
726  std::cout << "@@@ Fit::FillMatricesWithCalibratedParameters" << std::endl;
727 
728  //---------- Count how many measurements
729  ALIint NolinMes = 0;
730  std::vector<Measurement*>::const_iterator vmcite;
731  for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
732  NolinMes += (*vmcite)->dim();
733  }
734  if (ALIUtils::debug >= 4)
735  std::cout << "@@FillMatricesWithCalibratedParameters" << std::endl;
736 
737  std::vector<Entry*>::const_iterator vecite;
738  ALIint nEntcal = 0;
739  //---------- Loop entries
740  for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
741  // (= No parameters to be fitted - No parameters 'unk' )
742  //- std::cout << "entry" << (*veite) << std::endl;
743  //----- Take entries of quality = 'cal'
744  if ((*vecite)->quality() == 1) {
745  //--- Matrix A: fill diagonals with 1. (derivatives of entry w.r.t itself)
746  ALIint lineNo = NolinMes + nEntcal;
747  ALIint columnNo = (*vecite)->fitPos(); //=? nEntcal
748  AMatrix->AddData(lineNo, columnNo, 1.);
749  if (ALIUtils::debug >= 4)
750  std::cout << "Fit::FillMatricesWithCalibratedParameters: AMatrix ( " << lineNo << " , " << columnNo
751  << ") = " << 1. << std::endl;
752 
753  //--- Matrix W: sigma*sigma
754  ALIdouble entsig = (*vecite)->sigma();
755  if (ALIUtils::debug >= 4)
756  std::cout << "Fit::FillMatricesWithCalibratedParameters: WMatrix ( " << lineNo << " , " << columnNo
757  << ") = " << entsig * entsig << std::endl;
758  WMatrix->AddData(lineNo, lineNo, entsig * entsig);
759 
760  //--- Matrix y & f: fill it with 0.
761  //op yMatrix->AddData( lineNo, 0, (*vecite)->value());
762  //op yfMatrix->AddData( lineNo, 0, (*vecite)->value());
763  ALIdouble calFit;
765  gomgr->getGlobalOptionValue("calParamInyfMatrix", calFit);
766  if (calFit) {
767  yfMatrix->AddData(lineNo, 0, -(*vecite)->valueDisplacementByFitting());
768  //- yfMatrix->AddData( lineNo, 0, (*vecite)->value() );
769  //- yfMatrix->AddData( lineNo, 0, (*vecite)->lastAdditionToValueDisplacementByFitting() );
770  //- ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
771  // fileout << "cal to yf " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " " << (*vecite)->valueDisplacementByFitting() << endl;
772  // std::cout << "call to yf " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " " << (*vecite)->valueDisplacementByFitting() << std::endl;
773 
774  } else {
775  yfMatrix->AddData(lineNo, 0, 0.);
776  }
777  //t if(ALIUtils::debug >= 5) std::cout << "Fit::FillMatricesWithCalibratedParameters: yfMatrix ( " << lineNo << " , " << columnNo << ") = " << (*yfMatrix)(lineNo)(0) << std::endl;
778  nEntcal++;
779  }
780  }
781 }
782 
783 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
784 //@@ Gets the only instance of Model
785 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
787  if (ALIUtils::debug >= 3)
788  std::cout << "@@@ Fit::setCorrelationsInWMatrix" << std::endl;
789 
790  //----- Check if there are correlations to input
792  ALIint siz = corrMgr->getNumberOfCorrelations();
793  if (siz == 0)
794  return;
795 
796  //----- Set correlations
797  ALIuint ii;
798  for (ii = 0; ii < ALIuint(siz); ii++) {
799  //t if(ALIUtils::debug >= 5) std::cout << "globaloption cmslink fit" << Model::GlobalOptions()["cms_link"] << std::endl;
800  ErrorCorrelation* corr = corrMgr->getCorrelation(ii);
801  setCorrelationFromParamFitted(corr->getEntry1(), corr->getEntry2(), corr->getCorrelation());
802  }
803 }
804 
805 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
806 //@@ set correlation between two entries of two OptOs
807 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
808 void Fit::setCorrelationFromParamFitted(const pss& entry1, const pss& entry2, ALIdouble correl) {
809  ALIint pmsize = WMatrix->NoLines();
810  ALIint fit_pos1 = Model::getEntryByName(entry1.first, entry1.second)->fitPos();
811  ALIint fit_pos2 = Model::getEntryByName(entry2.first, entry2.second)->fitPos();
812  std::cout << "CHECKsetCorrelatiFPF " << fit_pos1 << " " << fit_pos2 << std::endl;
813 
814  if (fit_pos1 >= 0 && fit_pos1 < pmsize && fit_pos2 >= 0 && fit_pos2 < pmsize) {
815  setCorrelationFromParamFitted(fit_pos1, fit_pos2, correl);
816  }
817 }
818 
819 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
820 void Fit::setCorrelationFromParamFitted(const ALIint fit_pos1, const ALIint fit_pos2, ALIdouble correl) {
821  // ALIdouble error1 = sqrt( (*WMatrix)(fit_pos1, fit_pos1) );
822  // ALIdouble error2 = sqrt( (*WMatrix)(fit_pos2, fit_pos2) );
823  WMatrix->SetCorrelation(fit_pos1, fit_pos2, correl);
824  std::cout << "setCorrelatiFPF " << fit_pos1 << " " << fit_pos2 << " " << correl << std::endl;
825 }
826 
827 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
828 //@@ multiply matrices needed for fit
829 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
831  if (ALIUtils::debug >= 3)
832  std::cout << "@@@ Fit::multiplyMatrices " << std::endl;
833  //---------- Calculate transpose of A
834  AtMatrix = new ALIMatrix(*AMatrix);
835  if (ALIUtils::debug >= 5)
836  AtMatrix->Dump("AtMatrix=A");
837  //- std::cout << "call transpose";
838  AtMatrix->transpose();
839  if (ALIUtils::debug >= 4)
840  AtMatrix->Dump("AtMatrix");
841 
842  //---------- Calculate At * W * A
843  AtWAMatrix = new ALIMatrix(0, 0);
844  // if(ALIUtils::debug >= 5) AtWAMatrix->Dump("AtWAMatrix=0");
846  if (ALIUtils::debug >= 5)
847  AtWAMatrix->Dump("AtWAMatrix");
848 
850 
851  //t AtWAMatrix->EliminateLines(0,48);
852  //t AtWAMatrix->EliminateColumns(0,48);
853  time_t now;
854  now = clock();
855  if (ALIUtils::debug >= 0)
856  std::cout << "TIME:BEFORE_INVERSE : " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
858 
859  /* std::cout << " DETERMINANT W " << m_norm1( AtWAMatrix->MatNonConst() ) << std::endl;
860  if( m_norm1( AtWAMatrix->MatNonConst() ) == 0 ) {
861  std::cout << " DETERMINANT W " << m_norm1( AtWAMatrix->MatNonConst() ) << std::endl;
862  std::exception();
863  } */
864 
865  AtWAMatrix->inverse();
866  if (ALIUtils::debug >= 4)
867  AtWAMatrix->Dump("inverse AtWAmatrix");
868  now = clock();
869  if (ALIUtils::debug >= 0)
870  std::cout << "TIME:AFTER_INVERSE : " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
872 
873  //op thePropagationMatrix = AtWAMatrix;
874 
875  //op VaMatrix = new ALIMatrix( *AtWAMatrix );
876 
877  //----- Print out propagated errors of parameters (=AtWA diagonal elements)
878  std::vector<Entry*>::const_iterator vecite;
879 
880  if (ALIUtils::debug >= 4) {
881  std::cout << "PARAM"
882  << " Optical Object "
883  << " entry name "
884  << " Param.Value "
885  << " Prog.Error"
886  << " Orig.Error" << std::endl;
887  }
888 
889  ALIint nEnt = 0;
890  for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
891  //------------------ Number of parameters 'cal'
892  // (= No parameters to be fitted - No parameters 'unk' )
893  if ((*vecite)->quality() >= theMinimumEntryQuality) {
894  if (ALIUtils::debug >= 4) {
895  std::cout << nEnt << "PARAM" << std::setw(26) << (*vecite)->OptOCurrent()->name().c_str() << std::setw(8) << " "
896  << (*vecite)->name().c_str() << " " << std::setw(8) << " " << (*vecite)->value() << " "
897  << std::setw(8) << sqrt(AtWAMatrix->Mat()->me[nEnt][nEnt]) / (*vecite)->OutputSigmaDimensionFactor()
898  << " " << (*vecite)->sigma() / (*vecite)->OutputSigmaDimensionFactor() << " Q" << (*vecite)->quality()
899  << std::endl;
900  }
901  nEnt++;
902  }
903  }
904 
905  if (ALIUtils::debug >= 5)
906  yfMatrix->Dump("PD(y-f)Matrix final");
907 }
908 
909 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
910 //@@ check also that the fit_quality = SMat(0,0) is smaller for each new iteration
911 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
913  if (ALIUtils::debug >= 3)
914  std::cout << "@@@ Fit::getFitQuality" << std::endl;
915 
916  double fit_quality = GetSChi2(true);
917 
918  //---------- Calculate DS = Variable to recognize convergence (distance to minimum)
919  ALIMatrix* DatMatrix = new ALIMatrix(*DaMatrix);
920  // delete DaMatrix; //op
921  DatMatrix->transpose();
922  if (ALIUtils::debug >= 5)
923  DatMatrix->Dump("DatMatrix");
924  //op ALIMatrix* DSMat = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix * *PDMatrix);
925  ALIMatrix* DSMat = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix * *yfMatrix);
926  if (ALIUtils::debug >= 5) {
927  ALIMatrix* DSMattemp = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix);
928  DSMattemp->Dump("DSMattempMatrix=Dat*At*W");
929  ALIMatrix* DSMattemp2 = new ALIMatrix(*AtMatrix * *WMatrix * *yfMatrix);
930  DSMattemp2->Dump("DSMattempMatrix2=At*W*yf");
931  ALIMatrix* DSMattemp3 = new ALIMatrix(*AtMatrix * *WMatrix);
932  DSMattemp3->Dump("DSMattempMatrix3=At*W");
933  AtMatrix->Dump("AtMatrix");
934  }
935 
936  /* for( int ii = 0; ii < DatMatrix->NoColumns(); ii++ ){
937  std::cout << ii << " DS term " << (*DatMatrix)(0,ii) * (*DSMattemp2)(ii,0) << std::endl;
938  }*/
939  // delete AtMatrix; //op
940  // delete WMatrix; //op
941 
942  //op if(ALIUtils::debug >= 5) (*PDMatrix).Dump("PDMatrix");
943  if (ALIUtils::debug >= 5)
944  (*yfMatrix).Dump("yfMatrix");
945  if (ALIUtils::debug >= 5)
946  DSMat->Dump("DSMatrix final");
947  // delete yfMatrix; //op
948 
949  ALIdouble fit_quality_cut = (*DSMat)(0, 0);
950  //- ALIdouble fit_quality_cut =std::abs( (*DSMat)(0,0) );
951  delete DSMat;
952  if (ALIUtils::debug >= 0)
954  << " Fit quality predicted improvement in distance to minimum is = " << fit_quality_cut << std::endl;
955  if (ALIUtils::report >= 2) {
957  fileout << theNoFitIterations
958  << " Fit quality predicted improvement in distance to minimum is = " << fit_quality_cut << std::endl;
959  }
960 
961  //- double fit_quality_cut = thePreviousIterationFitQuality - fit_quality;
962  //- double fit_quality_cut = fit_quality;
963  //- std::cout << " fit_quality_cut " << fit_quality_cut << " fit_quality " << fit_quality << std::endl;
964 
965  //----- Check quality
966  time_t now;
967  now = clock();
968  if (ALIUtils::debug >= 0)
969  std::cout << "TIME:QUALITY_CHECKED: " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
971 
972  FitQuality fitQuality;
973 
974  //----- Chi2 is bigger, bad
975  // if( theNoFitIterations != 0 && fit_quality_cut > 0. ) {
976  if (fit_quality_cut < 0.) {
977  fitQuality = FQchiSquareWorsened;
978  if (ALIUtils::debug >= 1)
979  std::cerr << "!!WARNING: Fit quality has worsened: Fit Quality now = " << fit_quality << " before "
980  << thePreviousIterationFitQuality << " diff " << fit_quality - thePreviousIterationFitQuality
981  << std::endl;
982 
983  //----- Chi2 is smaller, check if we make another iteration
984  } else {
985  ALIdouble rel_fit_quality = std::abs(thePreviousIterationFitQuality - fit_quality) / fit_quality;
986  //----- Small chi2 change: end
987  if ((fit_quality_cut < theFitQualityCut || rel_fit_quality < theRelativeFitQualityCut) && canBeGood) {
988  if (ALIUtils::debug >= 2)
989  std::cout << "$$ Fit::getFitQuality good " << fit_quality_cut << " <? " << theFitQualityCut << " || "
990  << rel_fit_quality << " <? " << theRelativeFitQualityCut << " GOOD " << canBeGood << std::endl;
991  fitQuality = FQsmallDistanceToMinimum;
992  if (ALIUtils::report >= 1) {
994  fileout << "STOP: SMALL IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut << " < "
995  << theFitQualityCut << " OR (RELATIVE) " << rel_fit_quality << " < " << theRelativeFitQualityCut
996  << std::endl;
997  }
998  if (ALIUtils::debug >= 4) {
999  std::cout << "STOP: SMALL IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut << " < "
1000  << theFitQualityCut << " OR (RELATIVE) " << rel_fit_quality << " < " << theRelativeFitQualityCut
1001  << std::endl;
1002  }
1003 
1004  //----- Big chi2 change: go to next iteration
1005  } else {
1006  if (ALIUtils::debug >= 2)
1007  std::cout << "$$ Fit::getFitQuality bad " << fit_quality_cut << " <? " << theFitQualityCut << " || "
1008  << rel_fit_quality << " <? " << theRelativeFitQualityCut << " GOOD " << canBeGood << std::endl;
1009  fitQuality = FQbigDistanceToMinimum;
1010  //----- set thePreviousIterationFitQuality for next iteration
1011  thePreviousIterationFitQuality = fit_quality;
1012 
1013  if (ALIUtils::report >= 2) {
1015  fileout << "CONTINUE: BIG IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut
1016  << " >= " << theFitQualityCut << " AND (RELATIVE) " << rel_fit_quality
1017  << " >= " << theRelativeFitQualityCut << std::endl;
1018  }
1019  if (ALIUtils::debug >= 4) {
1020  std::cout << "CONTINUE: BIG IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut
1021  << " >= " << theFitQualityCut << " AND (RELATIVE) " << rel_fit_quality
1022  << " >= " << theRelativeFitQualityCut << std::endl;
1023  }
1024  }
1025  }
1026 
1027  return fitQuality;
1028 }
1029 
1030 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1032  if (ALIUtils::debug >= 3)
1033  std::cout << "@@@ Fit::GetSChi2 useDa= " << useDa << std::endl;
1034 
1035  ALIMatrix* SMat = nullptr;
1036  if (useDa) {
1037  //----- Calculate variables to check quality of this set of parameters
1038 
1039  //----- Calculate Da = (At * W * A)-1 * At * W * (y-f)
1040  /*t DaMatrix = new ALIMatrix( *AtWAMatrix );
1041  *DaMatrix *= *AtMatrix * *WMatrix;
1042  if(ALIUtils::debug >= 5) DaMatrix->Dump("DaMatrix before yf ");
1043  *DaMatrix *= *yfMatrix;
1044  if(ALIUtils::debug >= 5) DaMatrix->Dump("DaMatrix");
1045  */
1046 
1047  DaMatrix = new ALIMatrix(0, 0);
1048  // if(ALIUtils::debug >= 5) AtWAMatrix->Dump("AtWAMatrix=0");
1049  *DaMatrix = (*AtWAMatrix * *AtMatrix * *WMatrix * *yfMatrix);
1050  if (ALIUtils::debug >= 5)
1051  DaMatrix->Dump("DaMatrix");
1052 
1053  //----- Calculate S = chi2 = Fit quality = r^T W r (r = residual = f + A*Da - y )
1054  //op ALIMatrix* tmpM = new ALIMatrix( *AMatrix * *DaMatrix + *PDMatrix );
1055  // ALIMatrix* tmpM = new ALIMatrix( *AMatrix * *DaMatrix + *yfMatrix );
1056 
1057  ALIMatrix* tmpM = new ALIMatrix(0, 0);
1058  *tmpM = *AMatrix * *DaMatrix - *yfMatrix;
1059  if (ALIUtils::debug >= 5)
1060  tmpM->Dump("A*Da + f - y Matrix ");
1061 
1062  ALIMatrix* tmptM = new ALIMatrix(*tmpM);
1063  tmptM->transpose();
1064  if (ALIUtils::debug >= 5)
1065  tmptM->Dump("tmptM after transpose");
1066  if (ALIUtils::debug >= 5)
1067  WMatrix->Dump("WMatrix");
1068 
1069  // std::cout << "smat " << std::endl;
1070  //o ALIMatrix* SMat = new ALIMatrix(*tmptM * *WMatrix * *tmpM);
1071  ALIMatrix* SMat1 = MatrixByMatrix(*tmptM, *WMatrix);
1072  // ALIMatrix* SMat1 = MatrixByMatrix(*AMatrix,*WMatrix);
1073  if (ALIUtils::debug >= 5)
1074  SMat1->Dump("(A*Da + f - y)^T * W Matrix");
1075  SMat = MatrixByMatrix(*SMat1, *tmpM);
1076  // std::cout << "smatc " << std::endl;
1077  delete tmpM;
1078  delete tmptM;
1079  if (ALIUtils::debug >= 5)
1080  SMat->Dump("SMatrix with Da");
1081  } else {
1082  ALIMatrix* yftMat = new ALIMatrix(*yfMatrix);
1083  yftMat->transpose();
1084  SMat = new ALIMatrix(*yftMat * *WMatrix * *yfMatrix);
1085  delete yftMat;
1086  if (ALIUtils::debug >= 5)
1087  SMat->Dump("SMatrix no Da");
1088  }
1089  ALIdouble fit_quality = (*SMat)(0, 0);
1090  delete SMat;
1091  if (ALIUtils::debug >= 5)
1092  std::cout << " GetSChi2 = " << fit_quality << std::endl;
1093 
1094  PrintChi2(fit_quality, !useDa);
1095 
1096  return fit_quality;
1097 }
1098 
1099 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1100 //@@ Correct entries with fitted values
1101 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1103  if (ALIUtils::debug >= 4) {
1104  std::cout << "@@ Adding correction (DaMatrix) to Entries " << std::endl;
1105  DaMatrix->Dump("DaMatrix =");
1106  }
1107 
1108  /*- Now there are other places where entries are changed
1109  if( ALIUtils::report >= 3 ) {
1110  ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
1111  fileout << std::endl << " CHANGE IN ENTRIES" << std::endl
1112  << " Optical Object Parameter xxxxxx " << std::endl;
1113  }*/
1114 
1115  ALIint nEnt = 0;
1116  std::vector<Entry*>::const_iterator vecite;
1117  for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
1118  //------------------ Number of parameters 'cal'
1119  // (= No parameters to be fitted - No parameters 'unk' )
1120  //- std::cout << "qual" << (*vecite)->quality() << theMinimumEntryQuality << std::endl;
1121  if ((*vecite)->quality() >= theMinimumEntryQuality) {
1122  if (ALIUtils::debug >= 5) {
1123  std::cout << std::endl
1124  << " @@@ PENTRY change " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " "
1125  << " change= " << (*DaMatrix)(nEnt, 0) << " value= " << (*vecite)->valueDisplacementByFitting()
1126  << std::endl;
1127  }
1128  /* if( ALIUtils::report >=3 ) {
1129  ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
1130  fileout << "dd" << std::setw(30) << (*vecite)->OptOCurrent()->name()
1131  << std::setw(8) << " " << (*vecite)->name() << " "
1132  << std::setw(8) << " " << (*DaMatrix)(nEnt,0) / (*vecite)->OutputValueDimensionFactor()
1133  << " " << (*vecite)->valueDisplacementByFitting() / (*vecite)->OutputValueDimensionFactor() << " fitpos " << (*vecite)->fitPos()
1134  << std::endl;
1135  }*/
1136 
1137  //----- Store this displacement
1138  if (ALIUtils::debug >= 4)
1139  std::cout << " old valueDisplacementByFitting " << (*vecite)->name() << " "
1140  << (*vecite)->valueDisplacementByFitting() << " original value " << (*vecite)->value() << std::endl;
1141 
1142  (*vecite)->addFittedDisplacementToValue((*DaMatrix)(nEnt, 0));
1143 
1144  if (ALIUtils::debug >= 4)
1145  std::cout << nEnt << " new valueDisplacementByFitting " << (*vecite)->OptOCurrent()->name() << " "
1146  << (*vecite)->name() << " = " << (*vecite)->valueDisplacementByFitting() << " "
1147  << (*DaMatrix)(nEnt, 0) << std::endl;
1148  nEnt++;
1149  }
1150  }
1151 }
1152 
1153 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1154 //@@ Delete the previous addition of fitted values (to try a new one daFactor times smaller )
1155 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1157  if (ALIUtils::debug >= 4) {
1158  std::cout << "@@ Fit::substractToHalfDaMatrixToEntries " << std::endl;
1159  }
1160 
1161  std::vector<Entry*>::const_iterator vecite;
1162  for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
1163  if ((*vecite)->quality() >= theMinimumEntryQuality) {
1164  //-- (*vecite)->addFittedDisplacementToValue( -(*DaMatrix)(nEnt,0) );!!! it is not substracting the new value of DaMatrix, but substracting the value that was added last iteration, with which the new value of DaMatrix has been calculated for this iteration
1165 
1166  ALIdouble lastadd = (*vecite)->lastAdditionToValueDisplacementByFitting() * factor;
1167  //- if( lastadd < 0 ) lastadd *= -1;
1168  (*vecite)->addFittedDisplacementToValue(-lastadd);
1169  (*vecite)->setLastAdditionToValueDisplacementByFitting(-(*vecite)->lastAdditionToValueDisplacementByFitting());
1170  // (*vecite)->substractToHalfFittedDisplacementToValue();
1171 
1172  if (ALIUtils::debug >= 4)
1173  std::cout << " new valueDisplacementByFitting " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name()
1174  << " = " << (*vecite)->valueDisplacementByFitting() << " " << std::endl;
1175  }
1176  }
1177 }
1178 
1179 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1180 //@@ Dump all the entries that have been fitted (those that were 'cal' or 'unk'
1181 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1182 void Fit::dumpFittedValues(ALIFileOut& fileout, ALIbool printErrors, ALIbool printOrig) {
1183  //---------- print
1184  if (ALIUtils::debug >= 0) {
1185  std::cout << "SRPRPOS "
1186  << " Optical Object "
1187  << " Parameter"
1188  << " Fit.Value "
1189  << " Orig.Value" << std::endl;
1190  }
1191  //---------- Dump header
1192  if (ALIUtils::debug >= 0)
1193  std::cout << std::endl << "FITTED VALUES " << std::endl;
1194  fileout << std::endl
1195  << "FITTED VALUES " << std::endl
1196  << "nEnt_unk"
1197  << " Optical Object"
1198  << " Parameter ";
1199  if (printErrors) {
1200  fileout << " value (+-error)"
1201  << " orig.val (+-error)";
1202  } else {
1203  fileout << " value "
1204  << " orig.val ";
1205  }
1206  fileout << " quality" << std::endl;
1207 
1208  //---------- Iterate over OptO list
1209  std::vector<Entry*> entries;
1210  // const Entry* entry;
1211  int ii, siz;
1212  std::vector<OpticalObject*>::const_iterator vocite;
1213  ALIstring sys = ALIstring("system");
1214  for (vocite = Model::OptOList().begin(); vocite != Model::OptOList().end(); ++vocite) {
1215  if ((*vocite)->type() == sys)
1216  continue;
1217 
1218  fileout << " %%%% Optical Object: " << (*vocite)->longName() << std::endl;
1219 
1220  entries = (*vocite)->CoordinateEntryList();
1221  siz = entries.size();
1222  if (siz != 6) {
1223  std::cerr << "!!! FATAL ERROR: strange number of coordinates = " << siz << std::endl;
1224  abort();
1225  }
1226 
1227  //----- Dump entry centre coordinates (centre in current coordinates of parent frame <> summ of displacements, as each displacement is done with a different rotation of parent frame)
1228  OpticalObject* opto = entries[0]->OptOCurrent();
1229  const OpticalObject* optoParent = opto->parent();
1230  printCentreInOptOFrame(opto, optoParent, fileout, printErrors, printOrig);
1231 
1232  //----- Dump entry angles coordinates
1233  printRotationAnglesInOptOFrame(opto, optoParent, fileout, printErrors, printOrig);
1234 
1235  //----- Dump extra entries
1236  entries = (*vocite)->ExtraEntryList();
1237  siz = entries.size();
1238  for (ii = 0; ii < siz; ii++) {
1239  double entryvalue = getEntryValue(entries[ii]);
1240  dumpEntryAfterFit(fileout, entries[ii], entryvalue, printErrors, printOrig);
1241  }
1242  }
1243 
1244  dumpEntryCorrelations(fileout);
1245 }
1246 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1247 //@@ Dump all the entries that have been fitted in reference frames of all ancestors
1248 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1249 void Fit::dumpFittedValuesInAllAncestorFrames(ALIFileOut& fileout, ALIbool printErrors, ALIbool printOrig) {
1250  //---------- print
1251  fileout << std::endl
1252  << "@@@@ FITTED VALUES IN ALL ANCESTORS " << std::endl
1253  << "nEnt_unk"
1254  << " Optical Object"
1255  << " Parameter ";
1256  if (printErrors) {
1257  fileout << " value (+-error)";
1258  if (printOrig) {
1259  fileout << " orig.val (+-error)";
1260  }
1261  } else {
1262  fileout << " value ";
1263  if (printOrig) {
1264  fileout << " orig.val ";
1265  }
1266  }
1267  fileout << " quality" << std::endl;
1268 
1269  //---------- Iterate over OptO list
1270  std::vector<Entry*> entries;
1271  std::vector<OpticalObject*>::const_iterator vocite;
1272  ALIstring sys = ALIstring("system");
1273  for (vocite = Model::OptOList().begin(); vocite != Model::OptOList().end(); ++vocite) {
1274  if ((*vocite)->type() == sys)
1275  continue;
1276 
1277  fileout << " %%%% Optical Object: " << (*vocite)->longName() << std::endl;
1278 
1279  entries = (*vocite)->CoordinateEntryList();
1280 
1281  //----- Dump entry centre coordinates (centre in current coordinates of parent frame <> summ of displacements, as each displacement is done with a different rotation of parent frame)
1282  OpticalObject* opto = *vocite;
1283  const OpticalObject* optoParent = opto->parent();
1284  do {
1285  fileout << " %% IN FRAME : " << optoParent->longName() << std::endl;
1286  printCentreInOptOFrame(opto, optoParent, fileout, printErrors, printOrig);
1287 
1288  //----- Dump entry angles coordinates
1289  printRotationAnglesInOptOFrame(opto, optoParent, fileout, printErrors, printOrig);
1290  optoParent = optoParent->parent();
1291  } while (optoParent);
1292  }
1293 }
1294 
1295 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1297  const OpticalObject* optoAncestor,
1298  ALIFileOut& fileout,
1299  ALIbool printErrors,
1300  ALIbool printOrig) {
1301  CLHEP::Hep3Vector centreLocal;
1302  if (optoAncestor->type() == "system") {
1303  centreLocal = opto->centreGlob();
1304  } else {
1305  centreLocal = opto->centreGlob() - optoAncestor->centreGlob();
1306  CLHEP::HepRotation parentRmGlobInv = inverseOf(optoAncestor->rmGlob());
1307  centreLocal = parentRmGlobInv * centreLocal;
1308  }
1309  if (ALIUtils::debug >= 2) {
1310  std::cout << "CENTRE LOCAL " << opto->name() << " " << centreLocal << " GLOBL " << opto->centreGlob()
1311  << " parent GLOB " << optoAncestor->centreGlob() << std::endl;
1312  ALIUtils::dumprm(optoAncestor->rmGlob(), " parent rm ");
1313  }
1314  std::vector<Entry*> entries = opto->CoordinateEntryList();
1315  for (ALIuint ii = 0; ii < 3; ii++) {
1316  /* double entryvalue = getEntryValue( entries[ii] );
1317  ALIdouble entryvalue;
1318  if( ii == 0 ) {
1319  entryvalue = centreLocal.x();
1320  }else if( ii == 1 ) {
1321  entryvalue = centreLocal.y();
1322  }else if( ii == 2 ) {
1323  entryvalue = centreLocal.z();
1324  }*/
1326  fileout, entries[ii], centreLocal[ii] / entries[ii]->OutputValueDimensionFactor(), printErrors, printOrig);
1327  }
1328 }
1329 
1330 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1332  const OpticalObject* optoAncestor,
1333  ALIFileOut& fileout,
1334  ALIbool printErrors,
1335  ALIbool printOrig) {
1336  std::vector<Entry*> entries = opto->CoordinateEntryList();
1337  std::vector<double> entryvalues = opto->getRotationAnglesInOptOFrame(optoAncestor, entries);
1338  //- std::cout << " after return entryvalues[0] " << entryvalues[0] << " entryvalues[1] " << entryvalues[1] << " entryvalues[2] " << entryvalues[2] << std::endl;
1339  for (ALIuint ii = 3; ii < entries.size(); ii++) {
1341  fileout, entries[ii], entryvalues[ii - 3] / entries[ii]->OutputValueDimensionFactor(), printErrors, printOrig);
1342  }
1343 }
1344 
1345 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1347  double entryvalue;
1348  if (entry->type() == "angles") {
1349  if (ALIUtils::debug >= 2)
1350  std::cout << "WARNING valueDisplacementByFitting has no sense for angles " << std::endl;
1351 
1352  // commenting out the following line as it is a dead assignment due to the
1353  // subsequent assignment below
1354  // -> silence static analyzer warnings, but leaving the commented line in
1355  // case someone wants to actively use this code again
1356 
1357  // entryvalue = -999;
1358  }
1359  entryvalue = (entry->value() + entry->valueDisplacementByFitting()) / entry->OutputValueDimensionFactor();
1360  return entryvalue;
1361 }
1362 
1363 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1365  ALIFileOut& fileout, const Entry* entry, double entryvalue, ALIbool printErrors, ALIbool printOrig) {
1366  //- std::cout << " Fit::dumpEntryAfterFit " << entryvalue << std::endl;
1367  ALIdouble dimv = entry->OutputValueDimensionFactor();
1368  ALIdouble dims = entry->OutputSigmaDimensionFactor();
1369  //----- Dump to std::cout
1370  if (ALIUtils::debug >= 3) {
1371  std::cout << "ENTRY: " << std::setw(30) << entry->OptOCurrent()->name() << std::setw(8) << " " << entry->name()
1372  << " " << std::setw(8) << (entry->value() + entry->valueDisplacementByFitting()) << " " << entry->value()
1373  << " Q" << entry->quality() << std::endl;
1374  }
1375 
1376  if (entry->quality() == 2) {
1377  fileout << "UNK: " << entry->fitPos() << " ";
1378  } else if (entry->quality() == 1) {
1379  fileout << "CAL: " << entry->fitPos() << " ";
1380  // fileout << "CAL: -1 ";
1381  } else {
1382  fileout << "FIX: -1 ";
1383  }
1384 
1385  fileout << std::setw(30) << entry->OptOCurrent()->name() << std::setw(8) << " " << entry->name() << " "
1386  << std::setw(8) << std::setprecision(8) << entryvalue;
1387  if (entry->quality() >= theMinimumEntryQuality) {
1388  if (printErrors)
1389  fileout << " +- " << std::setw(8) << sqrt(AtWAMatrix->Mat()->me[entry->fitPos()][entry->fitPos()]) / dims;
1390  } else {
1391  if (printErrors)
1392  fileout << " +- " << std::setw(8) << 0.;
1393  }
1394  if (printOrig) {
1395  fileout << std::setw(8) << " " << entry->value() / dimv;
1396  if (printErrors)
1397  fileout << " +- " << std::setw(8) << entry->sigma() / dims << " Q" << entry->quality();
1398 
1399  if (ALIUtils::report >= 2) {
1400  float dif = (entry->value() + entry->valueDisplacementByFitting()) / dimv - entry->value() / dimv;
1401  if (std::abs(dif) < 1.E-9)
1402  dif = 0.;
1403  fileout << " DIFF= " << dif;
1404  // << " == " << ( entry->value() + entry->valueDisplacementByFitting() ) / dimv - entryvalue << " @@ " << ( entry->value() + entry->valueDisplacementByFitting() ) / dimv << " @@ " << entryvalue;
1405  } else {
1406  // fileout << std::endl;
1407  }
1408  }
1409 
1410  fileout << std::endl;
1411 }
1412 
1413 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1415  //----- Only dump correlations bigger than a factor
1416  ALIdouble minCorrel = 1.E-6;
1417  //----- Dump correlations
1418  fileout << std::endl
1419  << "CORRELATION BETWEEN 'unk' ENTRIES: (>= " << minCorrel << " )" << std::endl
1420  << "No_1 No_2 correlation " << std::endl;
1421 
1422  ALIuint i1, i2;
1423  std::vector<Entry*>::iterator veite1, veite2;
1424  std::string E1, E2;
1425  for (veite1 = Model::EntryList().begin(); veite1 != Model::EntryList().end(); ++veite1) {
1426  if ((*veite1)->quality() == 0) {
1427  continue;
1428  } else if ((*veite1)->quality() == 1) {
1429  E1 = "C";
1430  } else if ((*veite1)->quality() == 2) {
1431  E1 = "U";
1432  }
1433  i1 = (*veite1)->fitPos();
1434 
1435  for (veite2 = veite1 + 1; veite2 != Model::EntryList().end(); ++veite2) {
1436  i2 = (*veite2)->fitPos();
1437  if ((*veite2)->quality() == 0) {
1438  continue;
1439  } else if ((*veite2)->quality() == 1) {
1440  E2 = "C";
1441  } else if ((*veite2)->quality() == 2) {
1442  E2 = "U";
1443  }
1444  ALIdouble corr = AtWAMatrix->Mat()->me[i1][i2];
1445  ALIdouble corrf = corr / sqrt(AtWAMatrix->Mat()->me[i1][i1]) / sqrt(AtWAMatrix->Mat()->me[i2][i2]);
1446  if (std::abs(corrf) >= minCorrel) {
1447  if (ALIUtils::debug >= 0) {
1448  std::cout << "CORR:" << E1 << "" << E2 << " (" << i1 << ")"
1449  << " (" << i2 << ")"
1450  << " " << corrf << std::endl;
1451  }
1452  fileout << "CORR:" << E1 << "" << E2 << " (" << i1 << ")"
1453  << " (" << i2 << ")"
1454  << " " << corrf << std::endl;
1455  }
1456  }
1457  }
1458  //------- Dump optical object list
1459  if (ALIUtils::debug >= 2)
1461 }
1462 
1463 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1464 //@@ Dump matrices used for the fit
1465 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1467  //----- Dump matrices for this iteration
1469  // ofstream matout("matrices.out");
1470  matout << std::endl << " @@@@@@@@@@@@@@@ Iteration No : " << theNoFitIterations << std::endl;
1471  AMatrix->ostrDump(matout, "Matrix A");
1472  AtMatrix->ostrDump(matout, "Matrix At");
1473  WMatrix->ostrDump(matout, "Matrix W");
1474  AtWAMatrix->ostrDump(matout, "Matrix AtWA");
1475  //op VaMatrix->ostrDump( matout, "Matrix Va" );
1476  DaMatrix->ostrDump(matout, "Matrix Da");
1477  yfMatrix->ostrDump(matout, "Matrix y");
1478  //op fMatrix->ostrDump( matout, "Matrix f" );
1479  //op thePropagationMatrix->ostrDump( matout, "propagation Matrix " );
1480  matout.close();
1481 }
1482 
1483 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1484 //@@ findEntryFitPosition
1485 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1486 ALIint Fit::findEntryFitPosition(const ALIstring& opto_name, const ALIstring& entry_name) {
1487  ALIint fitposi = -99;
1488 
1489  OpticalObject* opto = Model::getOptOByName(opto_name);
1490  //- std::cout << "OPTO = " << opto->name() << std::endl;
1491  std::vector<Entry*>::const_iterator vecite;
1492  for (vecite = opto->CoordinateEntryList().begin(); vecite < opto->CoordinateEntryList().end(); ++vecite) {
1493  //- std::cout << "ENTRYLIST" << (*vecite)->name() << std::endl;
1494  if ((*vecite)->name() == entry_name) {
1495  //- std::cout << "FOUND " << std::endl;
1496  fitposi = (*vecite)->fitPos();
1497  }
1498  }
1499  for (vecite = opto->ExtraEntryList().begin(); vecite < opto->ExtraEntryList().end(); ++vecite) {
1500  //- std::cout << "ENTRYLIST" << (*vecite)->name() << std::endl;
1501  if ((*vecite)->name() == entry_name) {
1502  //- std::cout << "FOUND " << std::endl;
1503  fitposi = (*vecite)->fitPos();
1504  }
1505  }
1506 
1507  if (fitposi == -99) {
1508  std::cerr << "!!EXITING: entry name not found: " << entry_name << std::endl;
1509  exit(2);
1510  } else {
1511  return fitposi;
1512  }
1513 }
1514 
1515 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1517  double chi2meas = 0;
1518  double chi2cal = 0;
1519  ALIint nMeas = 0, nUnk = 0;
1520 
1521  //----- Calculate the chi2 of measurements
1522  std::vector<Measurement*>::const_iterator vmcite;
1523  for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
1524  //--- Calculate Simulated Value Original
1525  for (ALIuint ii = 0; ii < ALIuint((*vmcite)->dim()); ii++) {
1526  nMeas++;
1527  double c2 = ((*vmcite)->value(ii) - (*vmcite)->valueSimulated(ii)) / (*vmcite)->sigma(ii);
1528  chi2meas += c2 * c2;
1529  if (ALIUtils::debug >= 2) {
1530  std::cout << c2 << " adding chi2meas " << chi2meas << " " << (*vmcite)->name() << ": " << ii
1531  << " (mm)R: " << (*vmcite)->value(ii) * 1000. << " S: " << (*vmcite)->valueSimulated(ii) * 1000.
1532  << " Diff= " << ((*vmcite)->value(ii) - (*vmcite)->valueSimulated(ii)) * 1000. << std::endl;
1533  }
1534  }
1535  }
1536 
1537  //----- Calculate the chi2 of calibrated parameters
1538  std::vector<Entry*>::iterator veite;
1539  for (veite = Model::EntryList().begin(); veite != Model::EntryList().end(); ++veite) {
1540  if ((*veite)->quality() == 2)
1541  nUnk++;
1542  if ((*veite)->quality() == 1) {
1543  double c2 = (*veite)->valueDisplacementByFitting() / (*veite)->sigma();
1544  //double c2 = (*veite)->value() / (*veite)->sigma();
1545  chi2cal += c2 * c2;
1546  if (ALIUtils::debug >= 2)
1547  std::cout << c2 << " adding chi2cal " << chi2cal << " " << (*veite)->OptOCurrent()->name() << " "
1548  << (*veite)->name() << std::endl;
1549  //- std::cout << " valueDisplacementByFitting " << (*veite)->valueDisplacementByFitting() << " sigma " << (*veite)->sigma() << std::endl;
1550  }
1551  }
1552 
1553  if (ALIUtils::report >= 1) {
1555  fileout << " Chi2= " << chi2meas + chi2cal << " / " << nMeas - nUnk << " dof "
1556  << " From measurements= " << chi2meas << " from calibrated parameters= " << chi2cal << std::endl;
1557  }
1558  if (ALIUtils::debug >= 3)
1559  std::cout << " quality Chi2 (no correlations) " << chi2meas + chi2cal << " " << chi2meas << " " << chi2cal
1560  << std::endl;
1561 
1562  if (!isFirst) {
1563  // double fit_quality_change = thePreviousIterationFitQuality - fit_quality;
1564 
1565  if (ALIUtils::debug >= 0) {
1566  std::cout << std::endl << "@@@@ Fit iteration " << theNoFitIterations << " ..." << std::endl;
1567  // std::cout << theNoFitIterations << " Chi2 improvement in this iteration = " << fit_quality_change << std::endl;
1568  }
1569  if (ALIUtils::report >= 1) {
1571  fileout << std::endl << "Fit iteration " << theNoFitIterations << " ..." << std::endl;
1572  // fileout << theNoFitIterations << " Chi2 improvement in this iteration = " << fit_quality_change << std::endl;
1573  }
1574  }
1575 
1576  //---- Print chi2
1577  if (ALIUtils::debug >= 0)
1578  std::cout << theNoFitIterations << " Chi2 after iteration = " << fit_quality << std::endl;
1579  if (ALIUtils::report >= 1) {
1580  //--------- Get report file handler
1582  fileout << theNoFitIterations << " Chi2 after iteration = " << fit_quality << std::endl;
1583  }
1584 }
1585 
1586 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1588  if (ALIUtils::debug >= 3)
1589  std::cout << "@@@ Fit::CheckIfFitPossible" << std::endl;
1590 
1591  //----- Check if there is an unknown parameter that is not affecting any measurement
1592  ALIint NolinMes = 0;
1593  std::vector<Measurement*>::const_iterator vmcite;
1594  for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
1595  NolinMes += (*vmcite)->dim();
1596  }
1597 
1598  std::vector<Entry*>::const_iterator vecite;
1599  for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
1600  if (ALIUtils::debug >= 4)
1601  std::cout << "Fit::CheckIfFitPossible looping for entry " << (*vecite)->longName() << std::endl;
1602  if ((*vecite)->quality() == 2) {
1603  ALIint nCol = (*vecite)->fitPos();
1604  //--- Check all measurements
1605  ALIbool noDepend = TRUE;
1606  if (ALIUtils::debug >= 4)
1607  std::cout << "Fit::CheckIfFitPossible looping for entry " << nCol << std::endl;
1608  for (ALIint ii = 0; ii < NolinMes; ii++) {
1609  if (ALIUtils::debug >= 5)
1610  std::cout << " Derivative= (" << ii << "," << nCol << ") = " << (*AMatrix)(ii, nCol) << std::endl;
1611 
1612  if (std::abs((*AMatrix)(ii, nCol)) > ALI_DBL_MIN) {
1613  if (ALIUtils::debug >= 5)
1614  std::cout << "Fit::CheckIfFitIsPossible " << nCol << " " << ii << " = " << (*AMatrix)(ii, nCol)
1615  << std::endl;
1616  noDepend = FALSE;
1617  break;
1618  }
1619  }
1620  if (noDepend) {
1621  throw cms::Exception("SDFError")
1622  << "!!!FATAL ERROR: Fit::CheckIfFitPossible() no measurement depends on unknown entry "
1623  << (*vecite)->OptOCurrent()->name() << "/" << (*vecite)->name() << std::endl
1624  << "!!! Fit will not be possible! " << std::endl;
1625  }
1626  }
1627  }
1628 
1629  //------ Check if there are two unknown entries that have the derivatives of all measurements w.r.t. then equal (or 100% proportional). In this case any value of the first entry can be fully compensated by another value in the second entry without change in any measurement ---> the two entries cannot be fitted!
1630 
1631  std::vector<Entry*>::const_iterator vecite1, vecite2;
1632  ALIint nLin = AMatrix->NoLines();
1633  ALIdouble derivPrec = ALI_DBL_MIN;
1634  //---------- Loop entries
1635  for (vecite1 = Model::EntryList().begin(); vecite1 != Model::EntryList().end(); ++vecite1) {
1636  if ((*vecite1)->quality() == 2) {
1637  vecite2 = vecite1;
1638  ++vecite2;
1639  for (; vecite2 != Model::EntryList().end(); ++vecite2) {
1640  if ((*vecite2)->quality() == 2) {
1641  ALIint fitpos1 = (*vecite1)->fitPos();
1642  ALIint fitpos2 = (*vecite2)->fitPos();
1643  if (ALIUtils::debug >= 5)
1644  std::cout << "Fit::CheckIfFitIsPossible checking " << (*vecite1)->longName() << " ( " << fitpos1 << " ) & "
1645  << (*vecite2)->longName() << " ( " << fitpos2 << " ) " << std::endl;
1646  ALIdouble prop = DBL_MAX;
1647  ALIbool isProp = TRUE;
1648  for (ALIint ii = 0; ii < nLin; ii++) {
1649  if (ALIUtils::debug >= 5)
1650  std::cout << "Fit::CheckIfFitIsPossible " << ii << " : " << (*AMatrix)(ii, fitpos1)
1651  << " ?= " << (*AMatrix)(ii, fitpos2) << std::endl;
1652  if (std::abs((*AMatrix)(ii, fitpos1)) < derivPrec) {
1653  if (std::abs((*AMatrix)(ii, fitpos2)) > derivPrec) {
1654  isProp = FALSE;
1655  break;
1656  }
1657  } else {
1658  ALIdouble propn = (*AMatrix)(ii, fitpos2) / (*AMatrix)(ii, fitpos1);
1659  if (prop != DBL_MAX && prop != propn) {
1660  isProp = FALSE;
1661  break;
1662  }
1663  prop = propn;
1664  }
1665  }
1666  if (isProp) {
1667  std::cerr << "!!!FATAL ERROR: Fit::CheckIfFitPossible() two entries for which the measurements have the "
1668  "same dependency (or 100% proportional) "
1669  << (*vecite1)->longName() << " & " << (*vecite2)->longName() << std::endl
1670  << "!!! Fit will not be possible! " << std::endl;
1671  throw cms::Exception("SDFError");
1672  }
1673  }
1674  }
1675  }
1676  }
1677 }
1678 
1679 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1681  int measProp = -1;
1682 
1683  std::set<ALIuint> columnsEqual;
1684  std::set<ALIuint> columnsEqualSave;
1685  ALIuint biggestColumn = 0;
1686  ALIdouble biggest = 0.;
1687  for (int ii = 0; ii < AMatrix->NoColumns(); ii++) {
1688  if (std::abs((*AMatrix)(measNo, ii)) > biggest) {
1689  biggest = std::abs((*AMatrix)(measNo, ii));
1690  biggestColumn = ii;
1691  }
1692  columnsEqualSave.insert(ii);
1693  }
1694 
1695  ALIdouble div;
1696 
1697  for (int jj = 0; jj < AMatrix->NoLines(); jj++) {
1698  if (jj == int(measNo))
1699  continue;
1700  columnsEqual = columnsEqualSave;
1701  // check if ratio of each column to 'biggestColumn' is the same as for the N measurement
1702  for (int ii = 0; ii < AMatrix->NoColumns(); ii++) {
1703  div = (*AMatrix)(measNo, ii) / (*AMatrix)(measNo, biggestColumn);
1704  if (std::abs((*AMatrix)(jj, ii)) > ALI_DBL_MIN &&
1705  std::abs(div - (*AMatrix)(jj, ii) / (*AMatrix)(jj, biggestColumn)) > ALI_DBL_MIN) {
1706  if (ALIUtils::debug >= 3)
1707  std::cout << "CheckIfMeasIsProportionalToAnother 2 columns = " << ii << " in " << measNo << " & " << jj
1708  << std::endl;
1709  } else {
1710  if (ALIUtils::debug >= 3)
1711  std::cout << "CheckIfMeasIsProportionalToAnother 2 columns != " << ii << " in " << measNo << " & " << jj
1712  << std::endl;
1713  // if it is not equal delete this column
1714  std::set<ALIuint>::iterator ite = columnsEqual.find(ii);
1715  if (ite != columnsEqual.end()) {
1716  columnsEqual.erase(ite);
1717  }
1718  }
1719  }
1720  // check if not all columns have been deleted from columnsEqual
1721  if (!columnsEqual.empty()) {
1722  if (ALIUtils::debug >= 3)
1723  std::cout << "CheckIfMeasIsProportionalToAnother " << measNo << " = " << jj << std::endl;
1724  measProp = jj;
1725  break;
1726  }
1727  }
1728 
1729  return measProp;
1730 }
1731 
1732 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1734  std::string measname = " ";
1735 
1736  std::cout << " imeas " << imeas << std::endl;
1737  int Aline = 0;
1738  std::vector<Measurement*>::const_iterator vmcite;
1739  for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
1740  for (ALIuint jj = 0; jj < ALIuint((*vmcite)->dim()); jj++) {
1741  if (Aline == imeas) {
1742  char ctmp[20];
1743  gcvt(jj, 10, ctmp);
1744  return ((*vmcite)->name()) + ":" + std::string(ctmp);
1745  }
1746  Aline++;
1747  }
1748  }
1749 
1750  std::cout << " return measname " << measname << std::endl;
1751  return measname;
1752 }
static void set_time_now(time_t now)
Definition: ALIUtils.h:40
static ALIbool only1
Definition: Measurement.h:212
static void multiplyMatrices()
Definition: Fit.cc:830
static ALIdouble theMinDaFactor
Definition: Fit.h:205
const CLHEP::Hep3Vector & centreGlob() const
Definition: OpticalObject.h:75
long double ALIdouble
Definition: CocoaGlobals.h:11
void AddFittedEntriesSet(FittedEntriesSet *fents)
static ALIMatrix * AMatrix
Definition: Fit.h:164
#define TRUE
Definition: scimark2.h:12
static void setCorrelationFromParamFitted(const pss &entry1, const pss &entry2, ALIdouble correl)
Definition: Fit.cc:808
static CocoaDaqReader * GetDaqReader()
static void setCorrelationsInWMatrix()
Definition: Fit.cc:786
std::pair< ALIstring, ALIstring > pss
Definition: Fit.h:25
const double ALI_DBL_MIN
Definition: CocoaGlobals.h:25
static cocoaStatus getCocoaStatus()
Definition: Model.h:63
Definition: Entry.h:18
static void setCocoaStatus(const cocoaStatus cs)
Definition: Model.h:64
static void printCentreInOptOFrame(const OpticalObject *opto, const OpticalObject *optoAncestor, ALIFileOut &fileout, ALIbool printErrors=true, ALIbool printOrig=true)
Definition: Fit.cc:1296
std::vector< double > getRotationAnglesInOptOFrame(const OpticalObject *optoAncestor, const std::vector< Entry *> &entries) const
ErrorCorrelation * getCorrelation(ALIint ii)
static ALIint theMinimumEntryQuality
Definition: Fit.h:188
static std::vector< OpticalObject * > & OptOList()
Definition: Model.h:84
static void dumpEntryAfterFit(ALIFileOut &fileout, const Entry *entry, double entryvalue, ALIbool printErrors=true, ALIbool printOrig=true)
Definition: Fit.cc:1364
bool DumpCocoaResults()
Definition: CocoaDBMgr.cc:49
static void dumprm(const CLHEP::HepRotation &rm, const std::string &msg, std::ostream &out=std::cout)
Definition: ALIUtils.cc:71
MatrixMeschach * MatrixByMatrix(const MatrixMeschach &mat1, const MatrixMeschach &mat2)
static ALIdouble getMaximumDeviationDerivative()
Definition: ALIUtils.h:89
int ALIint
Definition: CocoaGlobals.h:15
static OpticalObject * getOptOByName(const ALIstring &opto_name)
--— Find an OptO name in theOptOList and return a pointer to it
Definition: Model.cc:562
static ALIint debug
Definition: ALIUtils.h:34
static ALIMatrix * AtMatrix
Definition: Fit.h:165
static ALIstring & ReportFName()
the name of the report File
Definition: Model.h:99
static void setMaximumDeviationDerivative(ALIdouble val)
Definition: ALIUtils.h:90
static ErrorCorrelationMgr * getInstance()
const ALIstring & name() const
Definition: OpticalObject.h:58
#define FALSE
Definition: scimark2.h:16
static Fit & getInstance()
Definition: Fit.cc:69
static GlobalOptionMgr * getInstance()
static void setFittableEntries()
Definition: Fit.cc:382
void FillMeasurements()
MAT * MatNonConst() const
void dumpOptOs(std::ostream &out=std::cout) const
static Entry * getEntryByName(const ALIstring &opto_name, const ALIstring &entry_name)
--— Search an Entry name in the Entry* list and return a pointer to it
Definition: Model.cc:615
static double getEntryValue(const Entry *entry)
Definition: Fit.cc:1346
int nEvent
Definition: myFastSimVal.cc:44
static FittedEntriesReader * getFittedEntriesReader()
Definition: Model.h:254
ALIint NoLines() const
static ALIint _NoLinesA
Definition: Fit.h:177
const ALIstring longName() const
static OpticalObjectMgr * getInstance()
Get the only instance.
ALIint NoColumns() const
const std::vector< Entry * > & ExtraEntryList() const
Definition: OpticalObject.h:65
static ALIMatrix * AtWAMatrix
Definition: Fit.h:167
static ALIdouble theFitQualityCut
Definition: Fit.h:194
ALIint fitPos() const
Definition: Entry.h:58
int getGlobalOptionValue(const ALIstring &sstr, ALIdouble &val)
--— Search a string in theGlobalOptions and return 1 if found
bool ALIbool
Definition: CocoaGlobals.h:19
static void dumpFittedValuesInAllAncestorFrames(ALIFileOut &fileout, ALIbool printErrors=true, ALIbool printOrig=true)
Definition: Fit.cc:1249
dictionary corr
static FitQuality getFitQuality(const ALIbool canBeGood=TRUE)
Definition: Fit.cc:912
static void substractLastDisplacementToEntries(const ALIdouble factor)
Definition: Fit.cc:1156
Definition: Fit.h:32
static ALIint MaxNoFitIterations
Definition: Fit.h:202
T sqrt(T t)
Definition: SSEVec.h:23
static ALIint findEntryFitPosition(const ALIstring &opto_name, const ALIstring &entry_name)
Definition: Fit.cc:1486
static ALIMatrix * WMatrix
Definition: Fit.h:166
static std::string GetMeasurementName(int meas)
Definition: Fit.cc:1733
static void setFirstTime(ALIbool val)
Definition: ALIUtils.h:88
static ALIdouble thePreviousIterationFitQuality
Definition: Fit.h:191
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
static ALIMatrix * DaMatrix
Definition: Fit.h:169
static void PropagateErrors()
Definition: Fit.cc:457
virtual bool ReadNextEvent()=0
static ALIbool fitNextEvent(ALIuint &nEvent)
Definition: Fit.cc:144
static ALIMatrix * yfMatrix
Definition: Fit.h:174
static Fit * instance
Definition: Fit.h:162
static FittedEntriesManager * getInstance()
static ALIuint nEvent
Definition: Fit.h:208
void SetCorrelation(ALIint i1, ALIint i2, ALIdouble corr)
const std::vector< Entry * > & CoordinateEntryList() const
Definition: OpticalObject.h:63
isFirst
Definition: cuy.py:418
static void WriteVisualisationFiles()
Definition: Fit.cc:354
static ALIint _NoColumnsA
Definition: Fit.h:178
ii
Definition: cuy.py:589
MatrixMeschach ALIMatrix
static void setApply(ALIbool val)
const MAT * Mat() const
static void FillMatricesWithMeasurements()
Definition: Fit.cc:638
void ostrDump(std::ostream &fout, const ALIstring &mtext)
static void CreateMatrices()
Definition: Fit.cc:570
static void redoMatrices()
Definition: Fit.cc:442
static NtupleManager * getInstance()
static void dumpMatrices()
Definition: Fit.cc:1466
static void dumpFittedValues(ALIFileOut &fileout, ALIbool printErrors=true, ALIbool printOrig=true)
Definition: Fit.cc:1182
static void dumpEntryCorrelations(ALIFileOut &file)
Definition: Fit.cc:1414
const CLHEP::HepRotation & rmGlob() const
Definition: OpticalObject.h:83
static CocoaDBMgr * getInstance()
Definition: CocoaDBMgr.cc:38
void AddData(ALIuint col, ALIuint lin, ALIdouble data)
static void startFit()
Definition: Fit.cc:103
static ALIint report
Definition: ALIUtils.h:33
static ALIdouble cameraScaleFactor
Definition: Measurement.h:158
static ALIbool getFirstTime()
Definition: ALIUtils.h:87
static void deleteMatrices()
Definition: Fit.cc:553
static void CheckIfFitPossible()
Definition: Fit.cc:1587
std::map< ALIstring, ALIdouble, std::less< ALIstring > > & GlobalOptions()
static void printRotationAnglesInOptOFrame(const OpticalObject *opto, const OpticalObject *optoAncestor, ALIFileOut &fileout, ALIbool printErrors=true, ALIbool printOrig=true)
Definition: Fit.cc:1331
void FillOptObjects(MatrixMeschach *AtWAMatrix)
static int CheckIfMeasIsProportionalToAnother(ALIuint measNo)
Definition: Fit.cc:1680
Fit()
Definition: Fit.h:35
std::string ALIstring
Definition: CocoaGlobals.h:9
void FillFitParameters(MatrixMeschach *AtWAMatrix)
static ALIdouble theRelativeFitQualityCut
Definition: Fit.h:197
const OpticalObject * parent() const
Definition: OpticalObject.h:60
static void calculateSimulatedMeasurementsWithOriginalValues()
Definition: Fit.cc:532
static ALIstring & MatricesFName()
the name of the File for storing the matrices
Definition: Model.h:102
static void PrintChi2(ALIdouble fit_quality, ALIbool isFirst)
Definition: Fit.cc:1516
static FitQuality fitParameters(const double daFactor)
Definition: Fit.cc:407
static std::vector< Measurement * > & MeasurementList()
Definition: Model.h:88
FitQuality
Definition: Fit.h:30
static ALIFileOut & getInstance(const ALIstring &filename)
Definition: ALIFileOut.cc:18
static std::vector< Entry * > & EntryList()
Definition: Model.h:86
void Dump(const ALIstring &mtext)
static void addDaMatrixToEntries()
Definition: Fit.cc:1102
static ALIdouble GetSChi2(ALIbool useDa)
Definition: Fit.cc:1031
static ALIint theNoFitIterations
Definition: Fit.h:200
static void dumpDimensions(std::ofstream &fout)
Definition: ALIUtils.cc:308
void FillNtupleTree()
const ALIstring & type() const
Definition: OpticalObject.h:59
unsigned int ALIuint
Definition: CocoaGlobals.h:17
def exit(msg="")
static void FillMatricesWithCalibratedParameters()
Definition: Fit.cc:724
static time_t time_now()
Definition: ALIUtils.h:39