Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:44:25

0001 #include <FWCore/Utilities/interface/Exception.h>
0002 //   COCOA class implementation file
0003 //Id:  Fit.cc
0004 //CAT: Fit
0005 //
0006 //   History: v1.0
0007 //   Pedro Arce
0008 
0009 #include <cstdlib>
0010 #include <iomanip>
0011 #include <cmath>  // among others include also floating-point std::abs functions
0012 #include <ctime>
0013 #include <set>
0014 
0015 #include "Alignment/CocoaModel/interface/Model.h"
0016 #include "Alignment/CocoaModel/interface/OpticalObject.h"
0017 #include "Alignment/CocoaFit/interface/Fit.h"
0018 
0019 #include "Alignment/CocoaModel/interface/Measurement.h"
0020 #include "Alignment/CocoaModel/interface/Entry.h"
0021 #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
0022 #include "Alignment/CocoaUtilities/interface/ALIFileOut.h"
0023 #include "Alignment/CocoaUtilities/interface/GlobalOptionMgr.h"
0024 #include "Alignment/CocoaModel/interface/DeviationsFromFileSensor2D.h"
0025 #include "Alignment/CocoaFit/interface/NtupleManager.h"
0026 #include "Alignment/CocoaFit/interface/FittedEntriesManager.h"
0027 #include "Alignment/CocoaFit/interface/FittedEntriesSet.h"
0028 #ifdef COCOA_VIS
0029 #include "Alignment/CocoaVisMgr/interface/ALIVRMLMgr.h"
0030 #include "Alignment/IgCocoaFileWriter/interface/IgCocoaFileMgr.h"
0031 #endif
0032 #include "Alignment/CocoaModel/interface/OpticalObjectMgr.h"
0033 #include "Alignment/CocoaModel/interface/ErrorCorrelationMgr.h"
0034 #include "Alignment/CocoaModel/interface/ErrorCorrelation.h"
0035 #include "Alignment/CocoaModel/interface/FittedEntriesReader.h"
0036 #include "Alignment/CocoaDaq/interface/CocoaDaqReader.h"
0037 #include "Alignment/CocoaFit/interface/CocoaDBMgr.h"
0038 
0039 Fit* Fit::instance = nullptr;
0040 
0041 ALIMatrix* Fit::AMatrix;
0042 ALIMatrix* Fit::AtMatrix;
0043 ALIMatrix* Fit::WMatrix;
0044 ALIMatrix* Fit::AtWAMatrix;
0045 //op ALIMatrix* Fit::VaMatrix;
0046 ALIMatrix* Fit::DaMatrix;
0047 //op ALIMatrix* Fit::PDMatrix;
0048 //-ALIMatrix* Fit::VyMatrix;
0049 ALIMatrix* Fit::yfMatrix;
0050 //ALIMatrix* Fit::fMatrix;
0051 
0052 ALIint Fit::_NoLinesA;
0053 ALIint Fit::_NoColumnsA;
0054 //op ALIMatrix* Fit::thePropagationMatrix;
0055 
0056 ALIint Fit::theMinimumEntryQuality;
0057 ALIdouble Fit::thePreviousIterationFitQuality = DBL_MAX;
0058 ALIdouble Fit::theFitQualityCut = -1;
0059 ALIdouble Fit::theRelativeFitQualityCut = -1;
0060 ALIint Fit::theNoFitIterations;
0061 ALIint Fit::MaxNoFitIterations = -1;
0062 ALIdouble Fit::theMinDaFactor = 1.e-8;
0063 
0064 ALIuint Fit::nEvent = 1;
0065 
0066 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0067 //@@  Gets the only instance of Model
0068 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0069 Fit& Fit::getInstance() {
0070   if (!instance) {
0071     instance = new Fit;
0072     ALIdouble go;
0073     GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
0074 
0075     gomgr->getGlobalOptionValue("maxDeviDerivative", go);
0076     ALIUtils::setMaximumDeviationDerivative(go);
0077     if (ALIUtils::debug >= 3)
0078       std::cout << " Fit::maximum_deviation_derivative " << ALIUtils::getMaximumDeviationDerivative() << std::endl;
0079 
0080     gomgr->getGlobalOptionValue("maxNoFitIterations", go);
0081     MaxNoFitIterations = int(go);
0082 
0083     gomgr->getGlobalOptionValue("fitQualityCut", go);
0084     theFitQualityCut = go;
0085     if (ALIUtils::debug >= 3)
0086       std::cout << " theFitQualityCut " << theFitQualityCut << std::endl;
0087 
0088     gomgr->getGlobalOptionValue("RelativeFitQualityCut", go);
0089     theRelativeFitQualityCut = go;
0090     if (ALIUtils::debug >= 3)
0091       std::cout << " theRelativeFitQualityCut " << theRelativeFitQualityCut << std::endl;
0092 
0093     gomgr->getGlobalOptionValue("minDaFactor", go);
0094     theMinDaFactor = go;
0095   }
0096 
0097   return *instance;
0098 }
0099 
0100 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0101 //@@  startFit: steering method to make the fit
0102 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0103 void Fit::startFit() {
0104   //  Model::setCocoaStatus( COCOA_InitFit );
0105   NtupleManager* NTmgr = NtupleManager::getInstance();
0106   if (GlobalOptionMgr::getInstance()->GlobalOptions()["rootResults"] > 0) {
0107     NTmgr->BookNtuple();
0108   }
0109 
0110   ALIUtils::setFirstTime(true);
0111 
0112   WriteVisualisationFiles();
0113 
0114   GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
0115   for (;;) {
0116     bool bend = fitNextEvent(nEvent);
0117     if (gomgr->GlobalOptions()["writeDBOptAlign"] > 0 || gomgr->GlobalOptions()["writeDBAlign"] > 0) {
0118       CocoaDBMgr::getInstance()->DumpCocoaResults();
0119     }
0120 
0121     if (!bend) {
0122       if (ALIUtils::debug >= 1)
0123         std::cout << "@@@ Fit::startFit  ended  n events = " << nEvent << std::endl;
0124       break;
0125     }
0126 
0127     //-    if ( ALIUtils::debug >= 0) std::cout << " FIT STATUS " << Model::printCocoaStatus( Model::getCocoaStatus() ) << std::endl;
0128 
0129     nEvent++;
0130   }
0131 
0132   //---------- Program ended, fill histograms of fitted entries
0133   if (gomgr->GlobalOptions()["histograms"] > 0) {
0134     FittedEntriesManager* FEmgr = FittedEntriesManager::getInstance();
0135     FEmgr->MakeHistos();
0136   }
0137 
0138   if (GlobalOptionMgr::getInstance()->GlobalOptions()["rootResults"] > 0) {
0139     NTmgr->WriteNtuple();
0140   }
0141 }
0142 
0143 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0144 ALIbool Fit::fitNextEvent(ALIuint& nEvent) {
0145   if (Model::getFittedEntriesReader() != nullptr)
0146     Model::getFittedEntriesReader()->readFittedEntriesFromFile();
0147 
0148   //----- Reset coordinates to those read at the start
0149   std::vector<OpticalObject*>::iterator voite;
0150   for (voite = Model::OptOList().begin(); voite != Model::OptOList().end(); ++voite) {
0151     (*voite)->resetOriginalOriginalCoordinates();
0152   }
0153 
0154   //----- Reset entries displacements to 0.
0155   std::vector<Entry*>::iterator veite;
0156   for (veite = Model::EntryList().begin(); veite != Model::EntryList().end(); ++veite) {
0157     (*veite)->resetValueDisplacementByFitting();
0158   }
0159 
0160   ALIbool lastEvent = false;
0161 
0162   //-    DeviationsFromFileSensor2D::setApply( 1 );
0163 
0164   //m  ALIbool moreDataSets = Model::readMeasurementsFromFile( Measurement::only1Date, Measurement::only1Time );
0165 
0166   //----- Check if there are more data sets
0167   ALIbool moreDataSets = true;
0168   if (CocoaDaqReader::GetDaqReader() != nullptr)
0169     moreDataSets = CocoaDaqReader::GetDaqReader()->ReadNextEvent();
0170 
0171   if (ALIUtils::debug >= 5)
0172     std::cout << CocoaDaqReader::GetDaqReader() << "$$$$$$$$$$$$$$$ More Data Sets to be processed: " << moreDataSets
0173               << std::endl;
0174 
0175   if (moreDataSets) {
0176     if (ALIUtils::debug >= 2)
0177       std::cout << std::endl << "@@@@@@@@@@@@@@@@@@ Starting data set fit : " << nEvent << std::endl;
0178 
0179     //----- Count entries to be fitted, and set their order in theFitPos
0180     setFittableEntries();
0181 
0182     //----- Dump dimensions of output in 'report.out' file
0183     ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
0184     fileout << std::endl << "@@@@@@@ NEW MEASUREMENT SET " << nEvent << std::endl;
0185     if (ALIUtils::report >= 1)
0186       ALIUtils::dumpDimensions(fileout);
0187 
0188     //----- reset Number of iterations of non linear fit
0189     theNoFitIterations = 0;
0190 
0191     GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
0192     ALIdouble dumpMat;
0193     gomgr->getGlobalOptionValue("save_matrices", dumpMat);
0194 
0195     //----- Fit parameters
0196     double daFactor = 1.;
0197     Model::setCocoaStatus(COCOA_FirstIterationInEvent);
0198     for (;;) {
0199       if (ALIUtils::debug >= 2) {
0200         std::cout << std::endl << "Fit iteration " << theNoFitIterations << " ..." << std::endl;
0201       }
0202 
0203       //---------- Calculate the original simulated values of each Measurement (when all entries have their read in values)
0204       calculateSimulatedMeasurementsWithOriginalValues();  //?? original changed atfer each iteration
0205 
0206       FitQuality fq = fitParameters(daFactor);
0207       if (dumpMat > 1)
0208         dumpMatrices();
0209 
0210       //-      evaluateFitQuality( fq, daFactor );
0211 
0212       if (ALIUtils::debug >= 2) {
0213         std::cout << std::endl << "@@@@ Check fit quality for iteration " << theNoFitIterations << std::endl;
0214       }
0215 
0216       //----- Check if new iteration must be done
0217       if (fq == FQsmallDistanceToMinimum) {
0218         if (ALIUtils::debug >= 2)
0219           std::cout << std::endl << "@@@@ Fit quality: distance SMALLER than mininum " << std::endl;
0220         addDaMatrixToEntries();
0221         if (ALIUtils::report >= 1)
0222           dumpFittedValues(ALIFileOut::getInstance(Model::ReportFName()), TRUE, TRUE);
0223         //--- Print entries in all ancestor frames
0224         ALIdouble go;
0225         gomgr->getGlobalOptionValue("dumpInAllFrames", go);
0226         if (go >= 1)
0227           dumpFittedValuesInAllAncestorFrames(ALIFileOut::getInstance(Model::ReportFName()), FALSE, FALSE);
0228 
0229         break;  // No more iterations
0230       } else if (fq == FQbigDistanceToMinimum) {
0231         if (ALIUtils::debug >= 2)
0232           std::cout << std::endl << "@@@@ Fit quality: distance BIGGER than mininum " << std::endl;
0233         addDaMatrixToEntries();
0234         if (ALIUtils::report >= 1)
0235           dumpFittedValues(ALIFileOut::getInstance(Model::ReportFName()), TRUE, TRUE);
0236 
0237         //----- Next iteration (if not too many already)
0238         theNoFitIterations++;
0239         daFactor = 1.;
0240 
0241         //----- Too many iterations: end event here
0242         if (theNoFitIterations >= MaxNoFitIterations) {
0243           if (ALIUtils::debug >= 1)
0244             std::cerr << "!!!! WARNING: Too many iterations " << theNoFitIterations << "  and fit DOES NOT CONVERGE "
0245                       << std::endl;
0246 
0247           if (ALIUtils::report >= 2) {
0248             ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
0249             fileout << "!!!! WARNING: Too many iterations " << theNoFitIterations << "  and fit DOES NOT CONVERGE "
0250                     << std::endl;
0251           }
0252           //      Model::setCocoaStatus( COCOA_FitCannotImprove );
0253           break;  // No more iterations
0254         }
0255 
0256       } else if (fq == FQchiSquareWorsened) {
0257         if (ALIUtils::debug >= 1) {
0258           //----- Recalculate fit quality with decreasing values of Da
0259           std::cerr << "!! WARNING: fit quality has worsened, Recalculate fit quality with decreasing values of Da "
0260                     << std::endl;
0261           std::cout << " quality daFactor= " << daFactor << " minimum= " << theMinDaFactor << std::endl;
0262         }
0263         daFactor *= 0.5;
0264         if (daFactor > theMinDaFactor) {
0265           substractLastDisplacementToEntries(0.5);
0266 
0267           if (ALIUtils::report >= 2) {
0268             ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
0269             fileout << " Redoing iteration with Da factor " << daFactor << std::endl;
0270           }
0271         } else {
0272           daFactor *= 2.;
0273           std::cerr << " !!!ERROR: not possible to get good fit quality even multiplying Da by " << daFactor
0274                     << std::endl;
0275           if (ALIUtils::report >= 2) {
0276             ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
0277             fileout << " !!!ERROR: not possible to get good fit quality even multiplying Da by " << daFactor
0278                     << std::endl;
0279           }
0280           //        Model::setCocoaStatus( COCOA_FitCannotImprove );
0281           //-    std::cout << "fdsaf FIT STATUS " << Model::printCocoaStatus( Model::getCocoaStatus() ) << std::endl;
0282           break;  // No more iterations
0283         }
0284       }
0285       Model::setCocoaStatus(COCOA_NextIterationInEvent);
0286     }
0287 
0288     //----- Iteration is finished: dump fitted entries
0289     if (ALIUtils::debug >= 1)
0290       calculateSimulatedMeasurementsWithOriginalValues();
0291     if (gomgr->GlobalOptions()["histograms"] > 0) {
0292       FittedEntriesManager::getInstance()->AddFittedEntriesSet(new FittedEntriesSet(AtWAMatrix));
0293     }
0294 
0295     if (GlobalOptionMgr::getInstance()->GlobalOptions()["rootResults"] > 0) {
0296       NtupleManager* ntupleMgr = NtupleManager::getInstance();
0297       ntupleMgr->InitNtuple();
0298       ntupleMgr->FillChi2();
0299       ntupleMgr->FillOptObjects(AtWAMatrix);
0300       ntupleMgr->FillMeasurements();
0301       ntupleMgr->FillFitParameters(AtWAMatrix);
0302       ntupleMgr->FillNtupleTree();
0303     }
0304 
0305     //- only if not stopped in worsening quality state        if(ALIUtils::report >= 0) dumpFittedValues( ALIFileOut::getInstance( Model::ReportFName() ));
0306 
0307     /*-      std::vector< OpticalObject* >::iterator voite;
0308       for( voite = Model::OptOList().begin(); voite !=  Model::OptOList().end(); voite++ ) {
0309       //-??         (*voite)->resetOriginalOriginalCoordinates();
0310     }*/
0311 
0312     //---- If no measurement file, break after looping once
0313     //-      std::cout << " Measurement::measurementsFileName() " << Measurement::measurementsFileName() << " Measurement::measurementsFileName()" <<std::endl;
0314     if (CocoaDaqReader::GetDaqReader() == nullptr) {
0315       //m    if( Measurement::measurementsFileName() == "" ) {
0316       if (ALIUtils::debug >= 1)
0317         std::cout << std::endl << "@@@@@@@@@@@@@@@@@@ Fit has ended : only one measurement " << nEvent << std::endl;
0318       lastEvent = true;
0319       return !lastEvent;
0320     }
0321 
0322     //-      std::cout << "  Measurement::only1" <<  Measurement::only1 << std::endl;
0323     if (Measurement::only1) {
0324       if (ALIUtils::debug >= 1)
0325         std::cout << std::endl << "@@@@@@@@@@@@@@@@@@ Fit has ended : 'Measurement::only1'  is set" << std::endl;
0326 
0327       lastEvent = true;
0328       return !lastEvent;
0329     }
0330 
0331     if (GlobalOptionMgr::getInstance()->GlobalOptions()["maxEvents"] <= nEvent) {
0332       if (ALIUtils::debug >= 1)
0333         std::cout << std::endl
0334                   << "@@@@@@@@@@@@@@@@@@ Fit has ended : 'Number of events exhausted " << nEvent << std::endl;
0335 
0336       lastEvent = true;
0337       return !lastEvent;
0338     }
0339 
0340   } else {
0341     lastEvent = true;
0342     if (ALIUtils::debug >= 1)
0343       std::cout << std::endl << "@@@@@@@@@@@@@@@@@@ Fit has ended : ??no more data sets' " << nEvent << std::endl;
0344     return !lastEvent;
0345   }
0346 
0347   if (ALIUtils::debug >= 1)
0348     std::cout << std::endl << "@@@@@@@@@@@@@@@@@@ Fit has ended : " << nEvent << std::endl;
0349 
0350   return !lastEvent;
0351 }
0352 
0353 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0354 void Fit::WriteVisualisationFiles() {
0355 #ifdef COCOA_VIS
0356   if (gomgr->GlobalOptions()["VisOnly"] == 1) {
0357     calculateSimulatedMeasurementsWithOriginalValues();  //?? original changed atfer each iteration
0358   }
0359 
0360   GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
0361   if (gomgr->GlobalOptions()["VisWriteVRML"] > 0) {
0362     if (ALIUtils::getFirstTime())
0363       ALIVRMLMgr::getInstance().writeFile();
0364   }
0365   if (gomgr->GlobalOptions()["VisWriteIguana"] > 0) {
0366     if (ALIUtils::getFirstTime())
0367       IgCocoaFileMgr::getInstance().writeFile();
0368   }
0369 
0370   if (gomgr->GlobalOptions()["VisOnly"] == 1) {
0371     if (ALIUtils::debug >= 1)
0372       std::cout << " Visualiation file(s) succesfully written. Ending.... " << std::endl;
0373     exit(1);
0374   }
0375 #endif
0376 }
0377 
0378 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0379 //@@  Count how many entries are going to be fitted (have quality >=  theMinimumEntryQuality)
0380 //@@ Set for this entries the value of theFitPos
0381 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0382 void Fit::setFittableEntries() {
0383   std::vector<Entry*>::const_iterator vecite;
0384 
0385   GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
0386   theMinimumEntryQuality = int(gomgr->GlobalOptions()[ALIstring("calcul_type")]) + 1;
0387   if (ALIUtils::debug >= 3)
0388     std::cout << "@@@ Fit::setFittableEntries: total Entry List size= " << Model::EntryList().size() << std::endl;
0389 
0390   int No_entry_to_fit = 0;
0391   for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
0392     // Number the parameters that are going to be fitted
0393     if ((*vecite)->quality() >= theMinimumEntryQuality) {
0394       (*vecite)->setFitPos(No_entry_to_fit);
0395       if (ALIUtils::debug >= 4)
0396         std::cout << " Entry To Fit= " << No_entry_to_fit << " " << (*vecite)->OptOCurrent()->name() << " "
0397                   << (*vecite)->name() << "   with quality= " << (*vecite)->quality() << std::endl;
0398       No_entry_to_fit++;
0399     }
0400   }
0401 }
0402 
0403 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0404 //@@ Main method in class Fit
0405 //@@ fitParameters: get the parameters through the chi square fit
0406 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0407 FitQuality Fit::fitParameters(const double daFactor) {
0408   if (ALIUtils::debug >= 3)
0409     std::cout << "@@@ Fit::fitParameters: Fit quality daFactor " << daFactor << std::endl;
0410 
0411   redoMatrices();
0412 
0413   //---- Get chi2 of first iteration
0414   if (Model::getCocoaStatus() == COCOA_FirstIterationInEvent) {
0415     thePreviousIterationFitQuality = GetSChi2(false);
0416 
0417     GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
0418     if (gomgr->GlobalOptions()[ALIstring("stopAfter1stIteration")] == 1) {
0419       std::cout << "@!! STOPPED by user after 1st iteration " << std::endl;
0420       exit(1);
0421     }
0422   }
0423 
0424   /*  //---------- Open output file
0425   if( ALIUtils::report >= 1 ) {
0426     ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
0427     //t    fileout << " REPORT.OUT " << std::endl;
0428     //t    ALIUtils::dumpDimensions( fileout );
0429     fileout << std::endl << "Fit iteration " << theNoFitIterations << " ..." << std::endl;
0430     }*/
0431 
0432   //-    std::cout << "2 FIT STATUS " << Model::printCocoaStatus( Model::getCocoaStatus() ) << std::endl;
0433 
0434   if (ALIUtils::debug >= 10) {
0435     std::cout << std::endl << " End fitParameters " << theNoFitIterations << " ..." << std::endl;
0436   }
0437 
0438   return getFitQuality();
0439 }
0440 
0441 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0442 void Fit::redoMatrices() {
0443   if (ALIUtils::debug >= 3)
0444     std::cout << "@@@ Fit::redoMatrices" << std::endl;
0445 
0446   deleteMatrices();
0447 
0448   calculateSimulatedMeasurementsWithOriginalValues();
0449 
0450   PropagateErrors();
0451 }
0452 
0453 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0454 //@@  Propagate the Errors from the entries to the measurements
0455 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0456 //cocoaStatus
0457 void Fit::PropagateErrors() {
0458   if (ALIUtils::debug >= 3)
0459     std::cout << "@@@ Fit::PropagateErrors" << std::endl;
0460 
0461   //----- Create empty matrices of appropiate size
0462   CreateMatrices();
0463 
0464   //---- count running time
0465   time_t now;
0466   now = clock();
0467   if (ALIUtils::debug >= 2)
0468     std::cout << "TIME:CREATE_MAT    : " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
0469   ALIUtils::set_time_now(now);
0470 
0471   //----- Fill the A, W & y matrices with the measurements
0472   FillMatricesWithMeasurements();
0473 
0474   //---- count running time
0475   now = clock();
0476   if (ALIUtils::debug >= 2)
0477     std::cout << "TIME:MAT_MEAS_FILLED: " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
0478   ALIUtils::set_time_now(now);
0479 
0480   //----- Fill the A, W & y matrices with the calibrated parameters
0481   GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
0482   if (gomgr->GlobalOptions()[ALIstring("calcul_type")] == 0) {
0483     FillMatricesWithCalibratedParameters();
0484 
0485     //---- count running time
0486     now = clock();
0487     if (ALIUtils::debug >= 0)
0488       std::cout << "TIME:MAT_CAL_FILLED : " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
0489     ALIUtils::set_time_now(now);
0490   }
0491 
0492   //----- Put by hand some correlations if known previously
0493   setCorrelationsInWMatrix();
0494 
0495   if (ALIUtils::debug >= 3)
0496     WMatrix->Dump("WMatrix before inverse");
0497 
0498   //----- Check first that matrix can be inverted
0499   if (m_norm1(WMatrix->MatNonConst()) == 0) {
0500     //    Model::setCocoaStatus( COCOA_FitMatrixNonInversable );
0501     return;  //  Model::getCocoaStatus();
0502   } else {
0503     WMatrix->inverse();
0504   }
0505 
0506   if (ALIUtils::debug >= 3)
0507     AMatrix->Dump("AMatrix");
0508   if (ALIUtils::debug >= 3)
0509     WMatrix->Dump("WMatrix");
0510   if (ALIUtils::debug >= 3)
0511     yfMatrix->Dump("yfMatrix");
0512 
0513   if (gomgr->GlobalOptions()["onlyDeriv"] >= 1) {
0514     std::cout << "ENDING after derivatives are calculated ('onlyDeriv' option set)" << std::endl;
0515     exit(1);
0516   }
0517 
0518   multiplyMatrices();
0519 
0520   now = clock();
0521   if (ALIUtils::debug >= 0)
0522     std::cout << "TIME:MAT_MULTIPLIED : " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
0523   ALIUtils::set_time_now(now);
0524 
0525   if (ALIUtils::getFirstTime() == 1)
0526     ALIUtils::setFirstTime(false);
0527 }
0528 
0529 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0530 //@@ Calculate the simulated value of each Measurement propagating the LightRay when all the entries have their original values
0531 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0532 void Fit::calculateSimulatedMeasurementsWithOriginalValues() {
0533   //  if( ALIUtils::debug >= 4) OpticalObjectMgr::getInstance()->dumpOptOs();
0534 
0535   //---------- Set DeviationsFromFileSensor2D::apply true
0536   DeviationsFromFileSensor2D::setApply(true);
0537 
0538   if (ALIUtils::debug >= 3)
0539     std::cout << "@@@ Fit::calculateSimulatedMeasurementsWithOriginalValues" << std::endl;
0540   //---------- Loop Measurements
0541   std::vector<Measurement*>::const_iterator vmcite;
0542   for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
0543     //----- Calculate Simulated Value Original
0544     (*vmcite)->calculateOriginalSimulatedValue();
0545   }
0546 
0547   //---------- Set DeviationsFromFileSensor2D::apply false
0548   // 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
0549   DeviationsFromFileSensor2D::setApply(false);
0550 }
0551 
0552 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0553 void Fit::deleteMatrices() {
0554   //delete matrices created in previous iteration
0555   delete DaMatrix;
0556   delete AMatrix;
0557   delete WMatrix;
0558   delete yfMatrix;
0559   //op  delete fMatrix;
0560   delete AtMatrix;
0561   delete AtWAMatrix;
0562   //op  delete VaMatrix;
0563   //-  delete VyMatrix;
0564   //op  delete PDMatrix;
0565 }
0566 
0567 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0568 //@@  Calculate the NoLines & NoColumns and create matrices
0569 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0570 void Fit::CreateMatrices() {
0571   //---------- Count number of measurements
0572   ALIint NoMeas = 0;
0573   std::vector<Measurement*>::const_iterator vmcite;
0574   for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
0575     NoMeas += (*vmcite)->dim();
0576   }
0577   if (ALIUtils::debug >= 9)
0578     std::cout << "NOMEAS" << NoMeas << std::endl;
0579 
0580   //-------- Count number of 'cal'ibrated parameters
0581   ALIint nEnt_cal = 0;
0582   ALIint noent = 0;
0583   //-  std::cout << Model::EntryList().size() << std::endl;
0584   GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
0585   if (gomgr->GlobalOptions()["calcul_type"] == 0) {  // fit also 'cal' parameters
0586     //-  if( ALIUtils::debug >= 9) std::cout << "NOENTCALLL " << nEnt_cal << std::endl;
0587     if (ALIUtils::debug >= 5)
0588       std::cout << " Count number of 'cal'ibrated parameters " << std::endl;
0589     std::vector<Entry*>::iterator veite;
0590     for (veite = Model::EntryList().begin(); veite != Model::EntryList().end(); ++veite) {
0591       if ((*veite)->quality() == 1)
0592         nEnt_cal++;
0593       noent++;
0594       if (ALIUtils::debug >= 6) {
0595         std::cout << (*veite)->quality() << " " << (*veite)->OptOCurrent()->name() << " " << (*veite)->name()
0596                   << " # ENT CAL " << nEnt_cal << " # ENT " << noent << std::endl;
0597       }
0598     }
0599   }
0600 
0601   //---------- Count number parameters to be fitted ('cal' + 'unk')
0602   ALIint NoParamFit = 0;
0603   std::vector<Entry*>::const_iterator vecite;
0604   for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
0605     if ((*vecite)->quality() >= theMinimumEntryQuality) {
0606       NoParamFit++;
0607       if (ALIUtils::debug >= 99)
0608         std::cout << (*vecite)->quality() << (*vecite)->OptOCurrent()->name() << (*vecite)->name() << "NoParamFit"
0609                   << NoParamFit << std::endl;
0610       //      break;
0611     }
0612   }
0613 
0614   //---------- Create Matrices
0615   ALIint NoLinesA = NoMeas + nEnt_cal;
0616   ALIint NoColumnsA = NoParamFit;
0617   AMatrix = new ALIMatrix(NoLinesA, NoColumnsA);
0618 
0619   ALIint NoLinesW = NoLinesA;
0620   ALIint NoColumnsW = NoLinesA;
0621   WMatrix = new ALIMatrix(NoLinesW, NoColumnsW);
0622 
0623   ALIint NoLinesY = NoLinesA;
0624   //op  yMatrix = new ALIMatrix( NoLinesY, 1 );
0625   yfMatrix = new ALIMatrix(NoLinesY, 1);
0626 
0627   //op  fMatrix = new ALIMatrix( NoLinesY, 1 );
0628 
0629   if (ALIUtils::debug >= 4)
0630     std::cout << "CreateMatrices: NoLinesA = " << NoLinesA << " NoColumnsA = " << NoColumnsA << std::endl;
0631 }
0632 
0633 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0634 //@@  Loop Measurements:
0635 //@@    Fill Matrix A with derivatives respect to affecting entries
0636 //@@    Fill Matrix W, y & f with values and sigmas of measurement coordinate
0637 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0638 void Fit::FillMatricesWithMeasurements() {
0639   if (ALIUtils::debug >= 3)
0640     std::cout << "@@@ Fit::FillMatricesWithMeasurements" << std::endl;
0641 
0642   int Aline = 0;
0643 
0644   //---------- Loop Measurements
0645   std::vector<Measurement*>::const_iterator vmcite;
0646   std::vector<Entry*>::const_iterator vecite;
0647   for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
0648     if (ALIUtils::debug >= 5)
0649       std::cout << "FillMatricesWithMeasurements: measurement " << (*vmcite)->name() << " # entries affecting "
0650                 << (*vmcite)->affectingEntryList().size() << std::endl;
0651 
0652     //-------- Array of derivatives with respect to entries
0653     ALIint measdim = (*vmcite)->dim();
0654     std::vector<ALIdouble> derivRE;
0655     //    derivRE = new ALIdouble[measdim];
0656 
0657     //-------- Fill matrix A:
0658     //------ Loop only Entries affecting this Measurement
0659     //-std::cout << "number affecting entries: " << (*vmcite)->affectingEntryList().size() << std::endl;
0660     for (vecite = (*vmcite)->affectingEntryList().begin(); vecite != (*vmcite)->affectingEntryList().end(); ++vecite) {
0661       //-------- If good quality, get derivative of measurement with respect to this Entry
0662       if ((*vecite)->quality() >= theMinimumEntryQuality) {
0663         if (ALIUtils::debug >= 4) {
0664           //      std::cout << "FillMatricesWithMeasurements: filling element ( " << Aline << " - " << Aline+measdim-1 << " , " << (*vecite)->fitPos() << std::endl;
0665           std::cout << "entry affecting: " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << std::endl;
0666         }
0667         derivRE = (*vmcite)->DerivativeRespectEntry(*vecite);
0668         //---------- Fill matrix A with derivatives
0669         for (ALIuint jj = 0; jj < ALIuint(measdim); jj++) {
0670           AMatrix->AddData(Aline + jj, (*vecite)->fitPos(), derivRE[jj]);
0671           if (ALIUtils::debug >= 5)
0672             std::cout << "FillMatricesWithMeasurements: AMATRIX (" << Aline + jj << "," << (*vecite)->fitPos() << " = "
0673                       << derivRE[jj] << std::endl;
0674           //---------- Reset Measurement simulated_value
0675           (*vmcite)->setValueSimulated(jj, (*vmcite)->valueSimulated_orig(jj));
0676         }
0677       }
0678     }
0679     //    delete[] derivRE;
0680 
0681     //---------- Fill matrices W, y and f:
0682     //------ Loop Measurement coordinates
0683     for (ALIuint jj = 0; jj < ALIuint((*vmcite)->dim()); jj++) {
0684       ALIdouble sigma = (*vmcite)->sigma()[jj];
0685       if (sigma == 0.) {
0686         std::cerr << "EXITING: Measurement number " << vmcite - Model::MeasurementList().begin() << "has 0 error!!"
0687                   << std::endl;
0688       } else {
0689         //----- Fill W Matrix with inverse of sigma squared
0690         // multiply error by cameraScaleFactor
0691         ALIdouble sigmanew = sigma * Measurement::cameraScaleFactor;
0692         //  std::cout << Aline+jj << " WMATRIX FILLING " << sigmanew << " * " << Measurement::cameraScaleFactor << std::endl;
0693         WMatrix->AddData(Aline + jj, Aline + jj, (sigmanew * sigmanew));
0694       }
0695       //op //----- Fill Matrices y with measurement value
0696       //op yMatrix->AddData( Aline+jj, 0, (*vmcite)->value()[jj] );
0697       //op //----- Fill f Matrix with simulated_value
0698       //op fMatrix->AddData( Aline+jj, 0, (*vmcite)->valueSimulated_orig(jj) );
0699       //----- Fill Matrix y - f with measurement value - simulated value
0700       yfMatrix->AddData(Aline + jj, 0, (*vmcite)->value()[jj] - (*vmcite)->valueSimulated_orig(jj));
0701       //      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;
0702     }
0703     if (ALIUtils::debug >= 99)
0704       std::cout << "change line" << Aline << std::endl;
0705     Aline += measdim;
0706     if (ALIUtils::debug >= 99)
0707       std::cout << "change line" << Aline << std::endl;
0708   }
0709 
0710   if (ALIUtils::debug >= 4)
0711     AMatrix->Dump("Matrix A with meas");
0712   if (ALIUtils::debug >= 4)
0713     WMatrix->Dump("Matrix W with meas");
0714   if (ALIUtils::debug >= 4)
0715     yfMatrix->Dump("Matrix y with meas");
0716 }
0717 
0718 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0719 //@@  Loop Measurements:
0720 //@@    Fill Matrix A with derivatives respect to affecting entries
0721 //@@    Fill Matrix W, y & f with values and sigmas of measurement coordinate
0722 //@@
0723 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0724 void Fit::FillMatricesWithCalibratedParameters() {
0725   if (ALIUtils::debug >= 3)
0726     std::cout << "@@@ Fit::FillMatricesWithCalibratedParameters" << std::endl;
0727 
0728   //---------- Count how many measurements
0729   ALIint NolinMes = 0;
0730   std::vector<Measurement*>::const_iterator vmcite;
0731   for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
0732     NolinMes += (*vmcite)->dim();
0733   }
0734   if (ALIUtils::debug >= 4)
0735     std::cout << "@@FillMatricesWithCalibratedParameters" << std::endl;
0736 
0737   std::vector<Entry*>::const_iterator vecite;
0738   ALIint nEntcal = 0;
0739   //---------- Loop entries
0740   for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
0741     //                  (= No parameters to be fitted - No parameters 'unk' )
0742     //-    std::cout << "entry" << (*veite) << std::endl;
0743     //----- Take entries of quality = 'cal'
0744     if ((*vecite)->quality() == 1) {
0745       //--- Matrix A: fill diagonals with 1. (derivatives of entry w.r.t itself)
0746       ALIint lineNo = NolinMes + nEntcal;
0747       ALIint columnNo = (*vecite)->fitPos();  //=? nEntcal
0748       AMatrix->AddData(lineNo, columnNo, 1.);
0749       if (ALIUtils::debug >= 4)
0750         std::cout << "Fit::FillMatricesWithCalibratedParameters:  AMatrix ( " << lineNo << " , " << columnNo
0751                   << ") = " << 1. << std::endl;
0752 
0753       //--- Matrix W: sigma*sigma
0754       ALIdouble entsig = (*vecite)->sigma();
0755       if (ALIUtils::debug >= 4)
0756         std::cout << "Fit::FillMatricesWithCalibratedParameters:  WMatrix ( " << lineNo << " , " << columnNo
0757                   << ") = " << entsig * entsig << std::endl;
0758       WMatrix->AddData(lineNo, lineNo, entsig * entsig);
0759 
0760       //--- Matrix y & f: fill it with 0.
0761       //op      yMatrix->AddData( lineNo, 0, (*vecite)->value());
0762       //op      yfMatrix->AddData( lineNo, 0, (*vecite)->value());
0763       ALIdouble calFit;
0764       GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
0765       gomgr->getGlobalOptionValue("calParamInyfMatrix", calFit);
0766       if (calFit) {
0767         yfMatrix->AddData(lineNo, 0, -(*vecite)->valueDisplacementByFitting());
0768         //- yfMatrix->AddData( lineNo, 0, (*vecite)->value() );
0769         //-         yfMatrix->AddData( lineNo, 0, (*vecite)->lastAdditionToValueDisplacementByFitting() );
0770         //- ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
0771         //  fileout << "cal to yf " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " " << (*vecite)->valueDisplacementByFitting() << endl;
0772         //  std::cout << "call to yf " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " " << (*vecite)->valueDisplacementByFitting() << std::endl;
0773 
0774       } else {
0775         yfMatrix->AddData(lineNo, 0, 0.);
0776       }
0777       //t      if(ALIUtils::debug >= 5) std::cout << "Fit::FillMatricesWithCalibratedParameters:  yfMatrix ( " << lineNo << " , " << columnNo  << ") = " << (*yfMatrix)(lineNo)(0) << std::endl;
0778       nEntcal++;
0779     }
0780   }
0781 }
0782 
0783 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0784 //@@  Gets the only instance of Model
0785 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0786 void Fit::setCorrelationsInWMatrix() {
0787   if (ALIUtils::debug >= 3)
0788     std::cout << "@@@ Fit::setCorrelationsInWMatrix" << std::endl;
0789 
0790   //----- Check if there are correlations to input
0791   ErrorCorrelationMgr* corrMgr = ErrorCorrelationMgr::getInstance();
0792   ALIint siz = corrMgr->getNumberOfCorrelations();
0793   if (siz == 0)
0794     return;
0795 
0796   //----- Set correlations
0797   ALIuint ii;
0798   for (ii = 0; ii < ALIuint(siz); ii++) {
0799     //t    if(ALIUtils::debug >= 5) std::cout << "globaloption cmslink fit" << Model::GlobalOptions()["cms_link"] << std::endl;
0800     ErrorCorrelation* corr = corrMgr->getCorrelation(ii);
0801     setCorrelationFromParamFitted(corr->getEntry1(), corr->getEntry2(), corr->getCorrelation());
0802   }
0803 }
0804 
0805 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0806 //@@  set correlation between two entries of two OptOs
0807 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0808 void Fit::setCorrelationFromParamFitted(const pss& entry1, const pss& entry2, ALIdouble correl) {
0809   ALIint pmsize = WMatrix->NoLines();
0810   ALIint fit_pos1 = Model::getEntryByName(entry1.first, entry1.second)->fitPos();
0811   ALIint fit_pos2 = Model::getEntryByName(entry2.first, entry2.second)->fitPos();
0812   std::cout << "CHECKsetCorrelatiFPF " << fit_pos1 << " " << fit_pos2 << std::endl;
0813 
0814   if (fit_pos1 >= 0 && fit_pos1 < pmsize && fit_pos2 >= 0 && fit_pos2 < pmsize) {
0815     setCorrelationFromParamFitted(fit_pos1, fit_pos2, correl);
0816   }
0817 }
0818 
0819 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0820 void Fit::setCorrelationFromParamFitted(const ALIint fit_pos1, const ALIint fit_pos2, ALIdouble correl) {
0821   //  ALIdouble error1 = sqrt( (*WMatrix)(fit_pos1, fit_pos1) );
0822   // ALIdouble error2 = sqrt( (*WMatrix)(fit_pos2, fit_pos2) );
0823   WMatrix->SetCorrelation(fit_pos1, fit_pos2, correl);
0824   std::cout << "setCorrelatiFPF " << fit_pos1 << " " << fit_pos2 << " " << correl << std::endl;
0825 }
0826 
0827 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0828 //@@  multiply matrices needed for fit
0829 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0830 void Fit::multiplyMatrices() {
0831   if (ALIUtils::debug >= 3)
0832     std::cout << "@@@ Fit::multiplyMatrices " << std::endl;
0833   //---------- Calculate transpose of A
0834   AtMatrix = new ALIMatrix(*AMatrix);
0835   if (ALIUtils::debug >= 5)
0836     AtMatrix->Dump("AtMatrix=A");
0837   //-  std::cout << "call transpose";
0838   AtMatrix->transpose();
0839   if (ALIUtils::debug >= 4)
0840     AtMatrix->Dump("AtMatrix");
0841 
0842   //---------- Calculate At * W * A
0843   AtWAMatrix = new ALIMatrix(0, 0);
0844   //  if(ALIUtils::debug >= 5) AtWAMatrix->Dump("AtWAMatrix=0");
0845   *AtWAMatrix = *AtMatrix * *WMatrix * *AMatrix;
0846   if (ALIUtils::debug >= 5)
0847     AtWAMatrix->Dump("AtWAMatrix");
0848 
0849   CheckIfFitPossible();
0850 
0851   //t  AtWAMatrix->EliminateLines(0,48);
0852   //t AtWAMatrix->EliminateColumns(0,48);
0853   time_t now;
0854   now = clock();
0855   if (ALIUtils::debug >= 0)
0856     std::cout << "TIME:BEFORE_INVERSE : " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
0857   ALIUtils::set_time_now(now);
0858 
0859   /*  std::cout << " DETERMINANT W " <<  m_norm1( AtWAMatrix->MatNonConst() ) << std::endl;
0860   if( m_norm1( AtWAMatrix->MatNonConst() ) == 0 ) {
0861     std::cout << " DETERMINANT W " <<  m_norm1( AtWAMatrix->MatNonConst() ) << std::endl;
0862     std::exception();
0863     } */
0864 
0865   AtWAMatrix->inverse();
0866   if (ALIUtils::debug >= 4)
0867     AtWAMatrix->Dump("inverse AtWAmatrix");
0868   now = clock();
0869   if (ALIUtils::debug >= 0)
0870     std::cout << "TIME:AFTER_INVERSE  : " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
0871   ALIUtils::set_time_now(now);
0872 
0873   //op  thePropagationMatrix = AtWAMatrix;
0874 
0875   //op  VaMatrix = new ALIMatrix( *AtWAMatrix );
0876 
0877   //----- Print out propagated errors of parameters (=AtWA diagonal elements)
0878   std::vector<Entry*>::const_iterator vecite;
0879 
0880   if (ALIUtils::debug >= 4) {
0881     std::cout << "PARAM"
0882               << "        Optical Object "
0883               << "   entry name "
0884               << "      Param.Value "
0885               << " Prog.Error"
0886               << " Orig.Error" << std::endl;
0887   }
0888 
0889   ALIint nEnt = 0;
0890   ALIint nEntUnk = 0;
0891   for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
0892     //------------------ Number of parameters 'cal'
0893     //                  (= No parameters to be fitted - No parameters 'unk' )
0894     if ((*vecite)->quality() >= theMinimumEntryQuality) {
0895       if (ALIUtils::debug >= 4) {
0896         std::cout << nEnt << "PARAM" << std::setw(26) << (*vecite)->OptOCurrent()->name().c_str() << std::setw(8) << " "
0897                   << (*vecite)->name().c_str() << " " << std::setw(8) << " " << (*vecite)->value() << " "
0898                   << std::setw(8) << sqrt(AtWAMatrix->Mat()->me[nEnt][nEnt]) / (*vecite)->OutputSigmaDimensionFactor()
0899                   << " " << (*vecite)->sigma() / (*vecite)->OutputSigmaDimensionFactor() << " Q" << (*vecite)->quality()
0900                   << std::endl;
0901       }
0902       nEnt++;
0903     }
0904     if ((*vecite)->quality() == 2)
0905       nEntUnk++;
0906   }
0907 
0908   if (ALIUtils::debug >= 5)
0909     yfMatrix->Dump("PD(y-f)Matrix final");
0910 }
0911 
0912 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0913 //@@ check also that the fit_quality = SMat(0,0) is smaller for each new iteration
0914 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0915 FitQuality Fit::getFitQuality(const ALIbool canBeGood) {
0916   if (ALIUtils::debug >= 3)
0917     std::cout << "@@@ Fit::getFitQuality" << std::endl;
0918 
0919   double fit_quality = GetSChi2(true);
0920 
0921   //---------- Calculate DS = Variable to recognize convergence (distance to minimum)
0922   ALIMatrix* DatMatrix = new ALIMatrix(*DaMatrix);
0923   //  delete DaMatrix; //op
0924   DatMatrix->transpose();
0925   if (ALIUtils::debug >= 5)
0926     DatMatrix->Dump("DatMatrix");
0927   //op  ALIMatrix* DSMat = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix * *PDMatrix);
0928   ALIMatrix* DSMat = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix * *yfMatrix);
0929   if (ALIUtils::debug >= 5) {
0930     ALIMatrix* DSMattemp = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix);
0931     DSMattemp->Dump("DSMattempMatrix=Dat*At*W");
0932     ALIMatrix* DSMattemp2 = new ALIMatrix(*AtMatrix * *WMatrix * *yfMatrix);
0933     DSMattemp2->Dump("DSMattempMatrix2=At*W*yf");
0934     ALIMatrix* DSMattemp3 = new ALIMatrix(*AtMatrix * *WMatrix);
0935     DSMattemp3->Dump("DSMattempMatrix3=At*W");
0936     AtMatrix->Dump("AtMatrix");
0937   }
0938 
0939   /*  for( int ii = 0; ii < DatMatrix->NoColumns(); ii++ ){
0940     std::cout << ii << " DS term " << (*DatMatrix)(0,ii) * (*DSMattemp2)(ii,0) << std::endl;
0941     }*/
0942   //  delete AtMatrix; //op
0943   //  delete WMatrix; //op
0944 
0945   //op  if(ALIUtils::debug >= 5) (*PDMatrix).Dump("PDMatrix");
0946   if (ALIUtils::debug >= 5)
0947     (*yfMatrix).Dump("yfMatrix");
0948   if (ALIUtils::debug >= 5)
0949     DSMat->Dump("DSMatrix final");
0950   //  delete yfMatrix; //op
0951 
0952   ALIdouble fit_quality_cut = (*DSMat)(0, 0);
0953   //-  ALIdouble fit_quality_cut =std::abs( (*DSMat)(0,0) );
0954   delete DSMat;
0955   if (ALIUtils::debug >= 0)
0956     std::cout << theNoFitIterations
0957               << " Fit quality predicted improvement in distance to minimum is = " << fit_quality_cut << std::endl;
0958   if (ALIUtils::report >= 2) {
0959     ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
0960     fileout << theNoFitIterations
0961             << " Fit quality predicted improvement in distance to minimum is = " << fit_quality_cut << std::endl;
0962   }
0963 
0964   //-  double fit_quality_cut = thePreviousIterationFitQuality - fit_quality;
0965   //- double fit_quality_cut = fit_quality;
0966   //-  std::cout << "  fit_quality_cut " <<  fit_quality_cut << " fit_quality " << fit_quality << std::endl;
0967 
0968   //----- Check quality
0969   time_t now;
0970   now = clock();
0971   if (ALIUtils::debug >= 0)
0972     std::cout << "TIME:QUALITY_CHECKED: " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
0973   ALIUtils::set_time_now(now);
0974 
0975   FitQuality fitQuality;
0976 
0977   //----- Chi2 is bigger, bad
0978   //    if( theNoFitIterations != 0 && fit_quality_cut > 0. ) {
0979   if (fit_quality_cut < 0.) {
0980     fitQuality = FQchiSquareWorsened;
0981     if (ALIUtils::debug >= 1)
0982       std::cerr << "!!WARNING: Fit quality has worsened: Fit Quality now = " << fit_quality << " before "
0983                 << thePreviousIterationFitQuality << " diff " << fit_quality - thePreviousIterationFitQuality
0984                 << std::endl;
0985 
0986     //----- Chi2 is smaller, check if we make another iteration
0987   } else {
0988     ALIdouble rel_fit_quality = std::abs(thePreviousIterationFitQuality - fit_quality) / fit_quality;
0989     //----- Small chi2 change: end
0990     if ((fit_quality_cut < theFitQualityCut || rel_fit_quality < theRelativeFitQualityCut) && canBeGood) {
0991       if (ALIUtils::debug >= 2)
0992         std::cout << "$$ Fit::getFitQuality good " << fit_quality_cut << " <? " << theFitQualityCut << " || "
0993                   << rel_fit_quality << " <? " << theRelativeFitQualityCut << " GOOD " << canBeGood << std::endl;
0994       fitQuality = FQsmallDistanceToMinimum;
0995       if (ALIUtils::report >= 1) {
0996         ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
0997         fileout << "STOP: SMALL IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut << " < "
0998                 << theFitQualityCut << " OR (RELATIVE) " << rel_fit_quality << " < " << theRelativeFitQualityCut
0999                 << std::endl;
1000       }
1001       if (ALIUtils::debug >= 4) {
1002         std::cout << "STOP: SMALL IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut << " < "
1003                   << theFitQualityCut << " OR (RELATIVE) " << rel_fit_quality << " < " << theRelativeFitQualityCut
1004                   << std::endl;
1005       }
1006 
1007       //----- Big chi2 change: go to next iteration
1008     } else {
1009       if (ALIUtils::debug >= 2)
1010         std::cout << "$$ Fit::getFitQuality bad " << fit_quality_cut << " <? " << theFitQualityCut << " || "
1011                   << rel_fit_quality << " <? " << theRelativeFitQualityCut << " GOOD " << canBeGood << std::endl;
1012       fitQuality = FQbigDistanceToMinimum;
1013       //----- set thePreviousIterationFitQuality for next iteration
1014       thePreviousIterationFitQuality = fit_quality;
1015 
1016       if (ALIUtils::report >= 2) {
1017         ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
1018         fileout << "CONTINUE: BIG IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut
1019                 << " >= " << theFitQualityCut << " AND (RELATIVE) " << rel_fit_quality
1020                 << " >= " << theRelativeFitQualityCut << std::endl;
1021       }
1022       if (ALIUtils::debug >= 4) {
1023         std::cout << "CONTINUE: BIG IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut
1024                   << " >= " << theFitQualityCut << " AND (RELATIVE) " << rel_fit_quality
1025                   << " >= " << theRelativeFitQualityCut << std::endl;
1026       }
1027     }
1028   }
1029 
1030   return fitQuality;
1031 }
1032 
1033 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1034 ALIdouble Fit::GetSChi2(ALIbool useDa) {
1035   if (ALIUtils::debug >= 3)
1036     std::cout << "@@@ Fit::GetSChi2  useDa= " << useDa << std::endl;
1037 
1038   ALIMatrix* SMat = nullptr;
1039   if (useDa) {
1040     //----- Calculate variables to check quality of this set of parameters
1041 
1042     //----- Calculate Da = (At * W * A)-1 * At * W * (y-f)
1043     /*t  DaMatrix = new ALIMatrix( *AtWAMatrix );
1044      *DaMatrix *= *AtMatrix * *WMatrix;
1045      if(ALIUtils::debug >= 5) DaMatrix->Dump("DaMatrix before yf ");
1046      *DaMatrix *= *yfMatrix;
1047      if(ALIUtils::debug >= 5) DaMatrix->Dump("DaMatrix");
1048     */
1049 
1050     DaMatrix = new ALIMatrix(0, 0);
1051     //  if(ALIUtils::debug >= 5) AtWAMatrix->Dump("AtWAMatrix=0");
1052     *DaMatrix = (*AtWAMatrix * *AtMatrix * *WMatrix * *yfMatrix);
1053     if (ALIUtils::debug >= 5)
1054       DaMatrix->Dump("DaMatrix");
1055 
1056     //----- Calculate S = chi2 = Fit quality = r^T W r (r = residual = f + A*Da - y )
1057     //op  ALIMatrix* tmpM = new ALIMatrix( *AMatrix * *DaMatrix + *PDMatrix );
1058     //  ALIMatrix* tmpM = new ALIMatrix( *AMatrix * *DaMatrix + *yfMatrix );
1059 
1060     ALIMatrix* tmpM = new ALIMatrix(0, 0);
1061     *tmpM = *AMatrix * *DaMatrix - *yfMatrix;
1062     if (ALIUtils::debug >= 5)
1063       tmpM->Dump("A*Da + f - y Matrix ");
1064 
1065     ALIMatrix* tmptM = new ALIMatrix(*tmpM);
1066     tmptM->transpose();
1067     if (ALIUtils::debug >= 5)
1068       tmptM->Dump("tmptM after transpose");
1069     if (ALIUtils::debug >= 5)
1070       WMatrix->Dump("WMatrix");
1071 
1072     //  std::cout << "smat " << std::endl;
1073     //o  ALIMatrix* SMat = new ALIMatrix(*tmptM * *WMatrix * *tmpM);
1074     ALIMatrix* SMat1 = MatrixByMatrix(*tmptM, *WMatrix);
1075     //  ALIMatrix* SMat1 = MatrixByMatrix(*AMatrix,*WMatrix);
1076     if (ALIUtils::debug >= 5)
1077       SMat1->Dump("(A*Da + f - y)^T * W  Matrix");
1078     SMat = MatrixByMatrix(*SMat1, *tmpM);
1079     //  std::cout << "smatc " << std::endl;
1080     delete tmpM;
1081     delete tmptM;
1082     if (ALIUtils::debug >= 5)
1083       SMat->Dump("SMatrix with Da");
1084   } else {
1085     ALIMatrix* yftMat = new ALIMatrix(*yfMatrix);
1086     yftMat->transpose();
1087     SMat = new ALIMatrix(*yftMat * *WMatrix * *yfMatrix);
1088     delete yftMat;
1089     if (ALIUtils::debug >= 5)
1090       SMat->Dump("SMatrix no Da");
1091   }
1092   ALIdouble fit_quality = (*SMat)(0, 0);
1093   delete SMat;
1094   if (ALIUtils::debug >= 5)
1095     std::cout << " GetSChi2 = " << fit_quality << std::endl;
1096 
1097   PrintChi2(fit_quality, !useDa);
1098 
1099   return fit_quality;
1100 }
1101 
1102 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1103 //@@ Correct entries with fitted values
1104 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1105 void Fit::addDaMatrixToEntries() {
1106   if (ALIUtils::debug >= 4) {
1107     std::cout << "@@ Adding correction (DaMatrix) to Entries " << std::endl;
1108     DaMatrix->Dump("DaMatrix =");
1109   }
1110 
1111   /*-  Now there are other places where entries are changed
1112     if( ALIUtils::report >= 3 ) {
1113     ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
1114     fileout << std::endl << " CHANGE IN ENTRIES" << std::endl
1115             << "          Optical Object       Parameter   xxxxxx " << std::endl;
1116         }*/
1117 
1118   ALIint nEnt = 0;
1119   std::vector<Entry*>::const_iterator vecite;
1120   for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
1121     //------------------ Number of parameters 'cal'
1122     //                  (= No parameters to be fitted - No parameters 'unk' )
1123     //-   std::cout << "qual" << (*vecite)->quality() << theMinimumEntryQuality << std::endl;
1124     if ((*vecite)->quality() >= theMinimumEntryQuality) {
1125       if (ALIUtils::debug >= 5) {
1126         std::cout << std::endl
1127                   << " @@@ PENTRY change " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " "
1128                   << " change= " << (*DaMatrix)(nEnt, 0) << " value= " << (*vecite)->valueDisplacementByFitting()
1129                   << std::endl;
1130       }
1131       /*      if( ALIUtils::report >=3 ) {
1132         ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
1133         fileout << "dd" << std::setw(30) << (*vecite)->OptOCurrent()->name() 
1134                 << std::setw(8) << " " << (*vecite)->name() << " " 
1135                 << std::setw(8) << " " << (*DaMatrix)(nEnt,0) / (*vecite)->OutputValueDimensionFactor()
1136             << " " << (*vecite)->valueDisplacementByFitting() / (*vecite)->OutputValueDimensionFactor() << " fitpos " << (*vecite)->fitPos()
1137                 << std::endl;
1138         }*/
1139 
1140       //----- Store this displacement
1141       if (ALIUtils::debug >= 4)
1142         std::cout << " old valueDisplacementByFitting " << (*vecite)->name() << " "
1143                   << (*vecite)->valueDisplacementByFitting() << " original value " << (*vecite)->value() << std::endl;
1144 
1145       (*vecite)->addFittedDisplacementToValue((*DaMatrix)(nEnt, 0));
1146 
1147       if (ALIUtils::debug >= 4)
1148         std::cout << nEnt << " new valueDisplacementByFitting " << (*vecite)->OptOCurrent()->name() << " "
1149                   << (*vecite)->name() << " = " << (*vecite)->valueDisplacementByFitting() << " "
1150                   << (*DaMatrix)(nEnt, 0) << std::endl;
1151       nEnt++;
1152     }
1153   }
1154 }
1155 
1156 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1157 //@@ Delete the previous addition of fitted values (to try a new one daFactor times smaller )
1158 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1159 void Fit::substractLastDisplacementToEntries(const ALIdouble factor) {
1160   if (ALIUtils::debug >= 4) {
1161     std::cout << "@@  Fit::substractToHalfDaMatrixToEntries " << std::endl;
1162   }
1163 
1164   std::vector<Entry*>::const_iterator vecite;
1165   for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
1166     if ((*vecite)->quality() >= theMinimumEntryQuality) {
1167       //--     (*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
1168 
1169       ALIdouble lastadd = (*vecite)->lastAdditionToValueDisplacementByFitting() * factor;
1170       //-      if( lastadd < 0 ) lastadd *= -1;
1171       (*vecite)->addFittedDisplacementToValue(-lastadd);
1172       (*vecite)->setLastAdditionToValueDisplacementByFitting(-(*vecite)->lastAdditionToValueDisplacementByFitting());
1173       //      (*vecite)->substractToHalfFittedDisplacementToValue();
1174 
1175       if (ALIUtils::debug >= 4)
1176         std::cout << " new valueDisplacementByFitting " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name()
1177                   << " = " << (*vecite)->valueDisplacementByFitting() << " " << std::endl;
1178     }
1179   }
1180 }
1181 
1182 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1183 //@@  Dump all the entries that have been fitted (those that were 'cal' or 'unk'
1184 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1185 void Fit::dumpFittedValues(ALIFileOut& fileout, ALIbool printErrors, ALIbool printOrig) {
1186   //---------- print
1187   if (ALIUtils::debug >= 0) {
1188     std::cout << "SRPRPOS "
1189               << "               Optical Object  "
1190               << "      Parameter"
1191               << " Fit.Value "
1192               << " Orig.Value" << std::endl;
1193   }
1194   //---------- Dump header
1195   if (ALIUtils::debug >= 0)
1196     std::cout << std::endl << "FITTED VALUES " << std::endl;
1197   fileout << std::endl
1198           << "FITTED VALUES " << std::endl
1199           << "nEnt_unk"
1200           << "             Optical Object"
1201           << "        Parameter  ";
1202   if (printErrors) {
1203     fileout << " value (+-error)"
1204             << " orig.val (+-error)";
1205   } else {
1206     fileout << " value "
1207             << " orig.val ";
1208   }
1209   fileout << " quality" << std::endl;
1210 
1211   //---------- Iterate over OptO list
1212   std::vector<Entry*> entries;
1213   //  const Entry* entry;
1214   int ii, siz;
1215   std::vector<OpticalObject*>::const_iterator vocite;
1216   for (vocite = Model::OptOList().begin(); vocite != Model::OptOList().end(); ++vocite) {
1217     if ((*vocite)->type() == ALIstring("system"))
1218       continue;
1219 
1220     fileout << " %%%% Optical Object: " << (*vocite)->longName() << std::endl;
1221 
1222     entries = (*vocite)->CoordinateEntryList();
1223     siz = entries.size();
1224     if (siz != 6) {
1225       std::cerr << "!!! FATAL ERROR: strange number of coordinates = " << siz << std::endl;
1226       abort();
1227     }
1228 
1229     //----- 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)
1230     OpticalObject* opto = entries[0]->OptOCurrent();
1231     const OpticalObject* optoParent = opto->parent();
1232     printCentreInOptOFrame(opto, optoParent, fileout, printErrors, printOrig);
1233 
1234     //----- Dump entry angles coordinates
1235     printRotationAnglesInOptOFrame(opto, optoParent, fileout, printErrors, printOrig);
1236 
1237     //----- Dump extra entries
1238     entries = (*vocite)->ExtraEntryList();
1239     siz = entries.size();
1240     for (ii = 0; ii < siz; ii++) {
1241       double entryvalue = getEntryValue(entries[ii]);
1242       dumpEntryAfterFit(fileout, entries[ii], entryvalue, printErrors, printOrig);
1243     }
1244   }
1245 
1246   dumpEntryCorrelations(fileout);
1247 }
1248 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1249 //@@  Dump all the entries that have been fitted in reference frames of all ancestors
1250 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1251 void Fit::dumpFittedValuesInAllAncestorFrames(ALIFileOut& fileout, ALIbool printErrors, ALIbool printOrig) {
1252   //---------- print
1253   fileout << std::endl
1254           << "@@@@ FITTED VALUES IN ALL ANCESTORS " << std::endl
1255           << "nEnt_unk"
1256           << "             Optical Object"
1257           << "        Parameter  ";
1258   if (printErrors) {
1259     fileout << " value (+-error)";
1260     if (printOrig) {
1261       fileout << " orig.val (+-error)";
1262     }
1263   } else {
1264     fileout << " value ";
1265     if (printOrig) {
1266       fileout << " orig.val ";
1267     }
1268   }
1269   fileout << " quality" << std::endl;
1270 
1271   //---------- Iterate over OptO list
1272   std::vector<Entry*> entries;
1273   std::vector<OpticalObject*>::const_iterator vocite;
1274   for (vocite = Model::OptOList().begin(); vocite != Model::OptOList().end(); ++vocite) {
1275     if ((*vocite)->type() == ALIstring("system"))
1276       continue;
1277 
1278     fileout << " %%%% Optical Object: " << (*vocite)->longName() << std::endl;
1279 
1280     entries = (*vocite)->CoordinateEntryList();
1281 
1282     //----- 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)
1283     OpticalObject* opto = *vocite;
1284     const OpticalObject* optoParent = opto->parent();
1285     do {
1286       fileout << " %% IN FRAME : " << optoParent->longName() << std::endl;
1287       printCentreInOptOFrame(opto, optoParent, fileout, printErrors, printOrig);
1288 
1289       //----- Dump entry angles coordinates
1290       printRotationAnglesInOptOFrame(opto, optoParent, fileout, printErrors, printOrig);
1291       optoParent = optoParent->parent();
1292     } while (optoParent);
1293   }
1294 }
1295 
1296 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1297 void Fit::printCentreInOptOFrame(const OpticalObject* opto,
1298                                  const OpticalObject* optoAncestor,
1299                                  ALIFileOut& fileout,
1300                                  ALIbool printErrors,
1301                                  ALIbool printOrig) {
1302   CLHEP::Hep3Vector centreLocal;
1303   if (optoAncestor->type() == "system") {
1304     centreLocal = opto->centreGlob();
1305   } else {
1306     centreLocal = opto->centreGlob() - optoAncestor->centreGlob();
1307     CLHEP::HepRotation parentRmGlobInv = inverseOf(optoAncestor->rmGlob());
1308     centreLocal = parentRmGlobInv * centreLocal;
1309   }
1310   if (ALIUtils::debug >= 2) {
1311     std::cout << "CENTRE LOCAL " << opto->name() << " " << centreLocal << " GLOBL " << opto->centreGlob()
1312               << " parent GLOB " << optoAncestor->centreGlob() << std::endl;
1313     ALIUtils::dumprm(optoAncestor->rmGlob(), " parent rm ");
1314   }
1315   std::vector<Entry*> entries = opto->CoordinateEntryList();
1316   for (ALIuint ii = 0; ii < 3; ii++) {
1317     /* double entryvalue = getEntryValue( entries[ii] );
1318        ALIdouble entryvalue;
1319        if( ii == 0 ) {
1320        entryvalue = centreLocal.x();
1321        }else if( ii == 1 ) {
1322        entryvalue = centreLocal.y();
1323        }else if( ii == 2 ) {
1324        entryvalue = centreLocal.z();
1325        }*/
1326     dumpEntryAfterFit(
1327         fileout, entries[ii], centreLocal[ii] / entries[ii]->OutputValueDimensionFactor(), printErrors, printOrig);
1328   }
1329 }
1330 
1331 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1332 void Fit::printRotationAnglesInOptOFrame(const OpticalObject* opto,
1333                                          const OpticalObject* optoAncestor,
1334                                          ALIFileOut& fileout,
1335                                          ALIbool printErrors,
1336                                          ALIbool printOrig) {
1337   std::vector<Entry*> entries = opto->CoordinateEntryList();
1338   std::vector<double> entryvalues = opto->getRotationAnglesInOptOFrame(optoAncestor, entries);
1339   //-    std::cout << " after return entryvalues[0] " << entryvalues[0] << " entryvalues[1] " << entryvalues[1] << " entryvalues[2] " << entryvalues[2] << std::endl;
1340   for (ALIuint ii = 3; ii < entries.size(); ii++) {
1341     dumpEntryAfterFit(
1342         fileout, entries[ii], entryvalues[ii - 3] / entries[ii]->OutputValueDimensionFactor(), printErrors, printOrig);
1343   }
1344 }
1345 
1346 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1347 double Fit::getEntryValue(const Entry* entry) {
1348   double entryvalue;
1349   if (entry->type() == "angles") {
1350     if (ALIUtils::debug >= 2)
1351       std::cout << "WARNING valueDisplacementByFitting has no sense for angles " << std::endl;
1352 
1353     // commenting out the following line as it is a dead assignment due to the
1354     // subsequent assignment below
1355     // -> silence static analyzer warnings, but leaving the commented line in
1356     //    case someone wants to actively use this code again
1357 
1358     // entryvalue = -999;
1359   }
1360   entryvalue = (entry->value() + entry->valueDisplacementByFitting()) / entry->OutputValueDimensionFactor();
1361   return entryvalue;
1362 }
1363 
1364 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1365 void Fit::dumpEntryAfterFit(
1366     ALIFileOut& fileout, const Entry* entry, double entryvalue, ALIbool printErrors, ALIbool printOrig) {
1367   //-  std::cout << " Fit::dumpEntryAfterFit " << entryvalue << std::endl;
1368   ALIdouble dimv = entry->OutputValueDimensionFactor();
1369   ALIdouble dims = entry->OutputSigmaDimensionFactor();
1370   //----- Dump to std::cout
1371   if (ALIUtils::debug >= 3) {
1372     std::cout << "ENTRY: " << std::setw(30) << entry->OptOCurrent()->name() << std::setw(8) << " " << entry->name()
1373               << " " << std::setw(8) << (entry->value() + entry->valueDisplacementByFitting()) << " " << entry->value()
1374               << " Q" << entry->quality() << std::endl;
1375   }
1376 
1377   if (entry->quality() == 2) {
1378     fileout << "UNK: " << entry->fitPos() << " ";
1379   } else if (entry->quality() == 1) {
1380     fileout << "CAL: " << entry->fitPos() << " ";
1381     //    fileout << "CAL: -1 ";
1382   } else {
1383     fileout << "FIX: -1 ";
1384   }
1385 
1386   fileout << std::setw(30) << entry->OptOCurrent()->name() << std::setw(8) << " " << entry->name() << " "
1387           << std::setw(8) << std::setprecision(8) << entryvalue;
1388   if (entry->quality() >= theMinimumEntryQuality) {
1389     if (printErrors)
1390       fileout << " +- " << std::setw(8) << sqrt(AtWAMatrix->Mat()->me[entry->fitPos()][entry->fitPos()]) / dims;
1391   } else {
1392     if (printErrors)
1393       fileout << " +- " << std::setw(8) << 0.;
1394   }
1395   if (printOrig) {
1396     fileout << std::setw(8) << " " << entry->value() / dimv;
1397     if (printErrors)
1398       fileout << " +- " << std::setw(8) << entry->sigma() / dims << " Q" << entry->quality();
1399 
1400     if (ALIUtils::report >= 2) {
1401       float dif = (entry->value() + entry->valueDisplacementByFitting()) / dimv - entry->value() / dimv;
1402       if (std::abs(dif) < 1.E-9)
1403         dif = 0.;
1404       fileout << " DIFF= " << dif;
1405       // << " == " << ( entry->value() + entry->valueDisplacementByFitting() )  / dimv - entryvalue << " @@ " << ( entry->value() + entry->valueDisplacementByFitting() ) / dimv << " @@ " <<  entryvalue;
1406     } else {
1407       //    fileout << std::endl;
1408     }
1409   }
1410 
1411   fileout << std::endl;
1412 }
1413 
1414 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1415 void Fit::dumpEntryCorrelations(ALIFileOut& fileout) {
1416   //----- Only dump correlations bigger than a factor
1417   ALIdouble minCorrel = 1.E-6;
1418   //----- Dump correlations
1419   fileout << std::endl
1420           << "CORRELATION BETWEEN 'unk' ENTRIES: (>= " << minCorrel << " )" << std::endl
1421           << "No_1  No_2   correlation " << std::endl;
1422 
1423   ALIuint i1, i2;
1424   std::vector<Entry*>::iterator veite1, veite2;
1425   std::string E1, E2;
1426   for (veite1 = Model::EntryList().begin(); veite1 != Model::EntryList().end(); ++veite1) {
1427     if ((*veite1)->quality() == 0) {
1428       continue;
1429     } else if ((*veite1)->quality() == 1) {
1430       E1 = "C";
1431     } else if ((*veite1)->quality() == 2) {
1432       E1 = "U";
1433     }
1434     i1 = (*veite1)->fitPos();
1435 
1436     for (veite2 = veite1 + 1; veite2 != Model::EntryList().end(); ++veite2) {
1437       i2 = (*veite2)->fitPos();
1438       if ((*veite2)->quality() == 0) {
1439         continue;
1440       } else if ((*veite2)->quality() == 1) {
1441         E2 = "C";
1442       } else if ((*veite2)->quality() == 2) {
1443         E2 = "U";
1444       }
1445       ALIdouble corr = AtWAMatrix->Mat()->me[i1][i2];
1446       ALIdouble corrf = corr / sqrt(AtWAMatrix->Mat()->me[i1][i1]) / sqrt(AtWAMatrix->Mat()->me[i2][i2]);
1447       if (std::abs(corrf) >= minCorrel) {
1448         if (ALIUtils::debug >= 0) {
1449           std::cout << "CORR:" << E1 << "" << E2 << " (" << i1 << ")"
1450                     << " (" << i2 << ")"
1451                     << " " << corrf << std::endl;
1452         }
1453         fileout << "CORR:" << E1 << "" << E2 << " (" << i1 << ")"
1454                 << " (" << i2 << ")"
1455                 << " " << corrf << std::endl;
1456       }
1457     }
1458   }
1459   //------- Dump optical object list
1460   if (ALIUtils::debug >= 2)
1461     OpticalObjectMgr::getInstance()->dumpOptOs();
1462 }
1463 
1464 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1465 //@@ Dump matrices used for the fit
1466 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1467 void Fit::dumpMatrices() {
1468   //----- Dump matrices for this iteration
1469   ALIFileOut& matout = ALIFileOut::getInstance(Model::MatricesFName());
1470   //  ofstream matout("matrices.out");
1471   matout << std::endl << " @@@@@@@@@@@@@@@  Iteration No : " << theNoFitIterations << std::endl;
1472   AMatrix->ostrDump(matout, "Matrix A");
1473   AtMatrix->ostrDump(matout, "Matrix At");
1474   WMatrix->ostrDump(matout, "Matrix W");
1475   AtWAMatrix->ostrDump(matout, "Matrix AtWA");
1476   //op  VaMatrix->ostrDump( matout, "Matrix Va" );
1477   DaMatrix->ostrDump(matout, "Matrix Da");
1478   yfMatrix->ostrDump(matout, "Matrix y");
1479   //op  fMatrix->ostrDump( matout, "Matrix f" );
1480   //op  thePropagationMatrix->ostrDump( matout, "propagation Matrix " );
1481   matout.close();
1482 }
1483 
1484 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1485 //@@ findEntryFitPosition
1486 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1487 ALIint Fit::findEntryFitPosition(const ALIstring& opto_name, const ALIstring& entry_name) {
1488   ALIint fitposi = -99;
1489 
1490   OpticalObject* opto = Model::getOptOByName(opto_name);
1491   //-  std::cout << "OPTO = " << opto->name() << std::endl;
1492   std::vector<Entry*>::const_iterator vecite;
1493   for (vecite = opto->CoordinateEntryList().begin(); vecite < opto->CoordinateEntryList().end(); ++vecite) {
1494     //-    std::cout << "ENTRYLIST" << (*vecite)->name() << std::endl;
1495     if ((*vecite)->name() == entry_name) {
1496       //-  std::cout << "FOUND " << std::endl;
1497       fitposi = (*vecite)->fitPos();
1498     }
1499   }
1500   for (vecite = opto->ExtraEntryList().begin(); vecite < opto->ExtraEntryList().end(); ++vecite) {
1501     //-    std::cout << "ENTRYLIST" << (*vecite)->name() << std::endl;
1502     if ((*vecite)->name() == entry_name) {
1503       //- std::cout << "FOUND " << std::endl;
1504       fitposi = (*vecite)->fitPos();
1505     }
1506   }
1507 
1508   if (fitposi == -99) {
1509     std::cerr << "!!EXITING: entry name not found: " << entry_name << std::endl;
1510     exit(2);
1511   } else {
1512     return fitposi;
1513   }
1514 }
1515 
1516 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1517 void Fit::PrintChi2(ALIdouble fit_quality, ALIbool isFirst) {
1518   double chi2meas = 0;
1519   double chi2cal = 0;
1520   ALIint nMeas = 0, nUnk = 0;
1521 
1522   //----- Calculate the chi2 of measurements
1523   std::vector<Measurement*>::const_iterator vmcite;
1524   for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
1525     //--- Calculate Simulated Value Original
1526     for (ALIuint ii = 0; ii < ALIuint((*vmcite)->dim()); ii++) {
1527       nMeas++;
1528       double c2 = ((*vmcite)->value(ii) - (*vmcite)->valueSimulated(ii)) / (*vmcite)->sigma(ii);
1529       chi2meas += c2 * c2;
1530       if (ALIUtils::debug >= 2) {
1531         std::cout << c2 << " adding chi2meas " << chi2meas << " " << (*vmcite)->name() << ": " << ii
1532                   << " (mm)R: " << (*vmcite)->value(ii) * 1000. << " S: " << (*vmcite)->valueSimulated(ii) * 1000.
1533                   << " Diff= " << ((*vmcite)->value(ii) - (*vmcite)->valueSimulated(ii)) * 1000. << std::endl;
1534       }
1535     }
1536   }
1537 
1538   //----- Calculate the chi2 of calibrated parameters
1539   std::vector<Entry*>::iterator veite;
1540   for (veite = Model::EntryList().begin(); veite != Model::EntryList().end(); ++veite) {
1541     if ((*veite)->quality() == 2)
1542       nUnk++;
1543     if ((*veite)->quality() == 1) {
1544       double c2 = (*veite)->valueDisplacementByFitting() / (*veite)->sigma();
1545       //double c2 = (*veite)->value() / (*veite)->sigma();
1546       chi2cal += c2 * c2;
1547       if (ALIUtils::debug >= 2)
1548         std::cout << c2 << " adding chi2cal " << chi2cal << " " << (*veite)->OptOCurrent()->name() << " "
1549                   << (*veite)->name() << std::endl;
1550       //-   std::cout << " valueDisplacementByFitting " << (*veite)->valueDisplacementByFitting() << " sigma " << (*veite)->sigma() << std::endl;
1551     }
1552   }
1553 
1554   if (ALIUtils::report >= 1) {
1555     ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
1556     fileout << " Chi2= " << chi2meas + chi2cal << " / " << nMeas - nUnk << " dof "
1557             << "  From measurements= " << chi2meas << " from calibrated parameters= " << chi2cal << std::endl;
1558   }
1559   if (ALIUtils::debug >= 3)
1560     std::cout << " quality Chi2 (no correlations) " << chi2meas + chi2cal << " " << chi2meas << " " << chi2cal
1561               << std::endl;
1562 
1563   if (!isFirst) {
1564     //    double fit_quality_change = thePreviousIterationFitQuality - fit_quality;
1565 
1566     if (ALIUtils::debug >= 0) {
1567       std::cout << std::endl << "@@@@ Fit iteration " << theNoFitIterations << " ..." << std::endl;
1568       //      std::cout << theNoFitIterations << " Chi2 improvement in this iteration = " << fit_quality_change << std::endl;
1569     }
1570     if (ALIUtils::report >= 1) {
1571       ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
1572       fileout << std::endl << "Fit iteration " << theNoFitIterations << " ..." << std::endl;
1573       //      fileout << theNoFitIterations << " Chi2 improvement in this iteration = " << fit_quality_change << std::endl;
1574     }
1575   }
1576 
1577   //---- Print chi2
1578   if (ALIUtils::debug >= 0)
1579     std::cout << theNoFitIterations << " Chi2 after iteration = " << fit_quality << std::endl;
1580   if (ALIUtils::report >= 1) {
1581     //--------- Get report file handler
1582     ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
1583     fileout << theNoFitIterations << " Chi2 after iteration = " << fit_quality << std::endl;
1584   }
1585 }
1586 
1587 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1588 void Fit::CheckIfFitPossible() {
1589   if (ALIUtils::debug >= 3)
1590     std::cout << "@@@ Fit::CheckIfFitPossible" << std::endl;
1591 
1592   //----- Check if there is an unknown parameter that is not affecting any measurement
1593   ALIint NolinMes = 0;
1594   std::vector<Measurement*>::const_iterator vmcite;
1595   for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
1596     NolinMes += (*vmcite)->dim();
1597   }
1598 
1599   std::vector<Entry*>::const_iterator vecite;
1600   for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
1601     if (ALIUtils::debug >= 4)
1602       std::cout << "Fit::CheckIfFitPossible looping for entry " << (*vecite)->longName() << std::endl;
1603     if ((*vecite)->quality() == 2) {
1604       ALIint nCol = (*vecite)->fitPos();
1605       //--- Check all measurements
1606       ALIbool noDepend = TRUE;
1607       if (ALIUtils::debug >= 4)
1608         std::cout << "Fit::CheckIfFitPossible looping for entry " << nCol << std::endl;
1609       for (ALIint ii = 0; ii < NolinMes; ii++) {
1610         if (ALIUtils::debug >= 5)
1611           std::cout << " Derivative= (" << ii << "," << nCol << ") = " << (*AMatrix)(ii, nCol) << std::endl;
1612 
1613         if (std::abs((*AMatrix)(ii, nCol)) > ALI_DBL_MIN) {
1614           if (ALIUtils::debug >= 5)
1615             std::cout << "Fit::CheckIfFitIsPossible " << nCol << " " << ii << " = " << (*AMatrix)(ii, nCol)
1616                       << std::endl;
1617           noDepend = FALSE;
1618           break;
1619         }
1620       }
1621       if (noDepend) {
1622         throw cms::Exception("SDFError")
1623             << "!!!FATAL ERROR: Fit::CheckIfFitPossible() no measurement depends on unknown entry "
1624             << (*vecite)->OptOCurrent()->name() << "/" << (*vecite)->name() << std::endl
1625             << "!!! Fit will not be possible! " << std::endl;
1626       }
1627     }
1628   }
1629 
1630   //------ 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!
1631 
1632   std::vector<Entry*>::const_iterator vecite1, vecite2;
1633   ALIint nLin = AMatrix->NoLines();
1634   ALIdouble derivPrec = ALI_DBL_MIN;
1635   //---------- Loop entries
1636   for (vecite1 = Model::EntryList().begin(); vecite1 != Model::EntryList().end(); ++vecite1) {
1637     if ((*vecite1)->quality() == 2) {
1638       vecite2 = vecite1;
1639       ++vecite2;
1640       for (; vecite2 != Model::EntryList().end(); ++vecite2) {
1641         if ((*vecite2)->quality() == 2) {
1642           ALIint fitpos1 = (*vecite1)->fitPos();
1643           ALIint fitpos2 = (*vecite2)->fitPos();
1644           if (ALIUtils::debug >= 5)
1645             std::cout << "Fit::CheckIfFitIsPossible checking " << (*vecite1)->longName() << " ( " << fitpos1 << " ) & "
1646                       << (*vecite2)->longName() << " ( " << fitpos2 << " ) " << std::endl;
1647           ALIdouble prop = DBL_MAX;
1648           ALIbool isProp = TRUE;
1649           for (ALIint ii = 0; ii < nLin; ii++) {
1650             if (ALIUtils::debug >= 5)
1651               std::cout << "Fit::CheckIfFitIsPossible " << ii << " : " << (*AMatrix)(ii, fitpos1)
1652                         << " ?= " << (*AMatrix)(ii, fitpos2) << std::endl;
1653             if (std::abs((*AMatrix)(ii, fitpos1)) < derivPrec) {
1654               if (std::abs((*AMatrix)(ii, fitpos2)) > derivPrec) {
1655                 isProp = FALSE;
1656                 break;
1657               }
1658             } else {
1659               ALIdouble propn = (*AMatrix)(ii, fitpos2) / (*AMatrix)(ii, fitpos1);
1660               if (prop != DBL_MAX && prop != propn) {
1661                 isProp = FALSE;
1662                 break;
1663               }
1664               prop = propn;
1665             }
1666           }
1667           if (isProp) {
1668             std::cerr << "!!!FATAL ERROR: Fit::CheckIfFitPossible() two entries for which the measurements have the "
1669                          "same dependency (or 100% proportional) "
1670                       << (*vecite1)->longName() << " & " << (*vecite2)->longName() << std::endl
1671                       << "!!! Fit will not be possible! " << std::endl;
1672             throw cms::Exception("SDFError");
1673           }
1674         }
1675       }
1676     }
1677   }
1678 }
1679 
1680 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1681 int Fit::CheckIfMeasIsProportionalToAnother(ALIuint measNo) {
1682   int measProp = -1;
1683 
1684   std::set<ALIuint> columnsEqual;
1685   std::set<ALIuint> columnsEqualSave;
1686   ALIuint biggestColumn = 0;
1687   ALIdouble biggest = 0.;
1688   for (int ii = 0; ii < AMatrix->NoColumns(); ii++) {
1689     if (std::abs((*AMatrix)(measNo, ii)) > biggest) {
1690       biggest = std::abs((*AMatrix)(measNo, ii));
1691       biggestColumn = ii;
1692     }
1693     columnsEqualSave.insert(ii);
1694   }
1695 
1696   ALIdouble div;
1697 
1698   for (int jj = 0; jj < AMatrix->NoLines(); jj++) {
1699     if (jj == int(measNo))
1700       continue;
1701     columnsEqual = columnsEqualSave;
1702     // check if ratio of each column to 'biggestColumn' is the same as for the N measurement
1703     for (int ii = 0; ii < AMatrix->NoColumns(); ii++) {
1704       div = (*AMatrix)(measNo, ii) / (*AMatrix)(measNo, biggestColumn);
1705       if (std::abs((*AMatrix)(jj, ii)) > ALI_DBL_MIN &&
1706           std::abs(div - (*AMatrix)(jj, ii) / (*AMatrix)(jj, biggestColumn)) > ALI_DBL_MIN) {
1707         if (ALIUtils::debug >= 3)
1708           std::cout << "CheckIfMeasIsProportionalToAnother 2 columns = " << ii << " in " << measNo << " & " << jj
1709                     << std::endl;
1710       } else {
1711         if (ALIUtils::debug >= 3)
1712           std::cout << "CheckIfMeasIsProportionalToAnother 2 columns != " << ii << " in " << measNo << " & " << jj
1713                     << std::endl;
1714         // if it is not equal delete this column
1715         std::set<ALIuint>::iterator ite = columnsEqual.find(ii);
1716         if (ite != columnsEqual.end()) {
1717           columnsEqual.erase(ite);
1718         }
1719       }
1720     }
1721     // check if not all columns have been deleted from columnsEqual
1722     if (!columnsEqual.empty()) {
1723       if (ALIUtils::debug >= 3)
1724         std::cout << "CheckIfMeasIsProportionalToAnother " << measNo << " = " << jj << std::endl;
1725       measProp = jj;
1726       break;
1727     }
1728   }
1729 
1730   return measProp;
1731 }
1732 
1733 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1734 std::string Fit::GetMeasurementName(int imeas) {
1735   std::string measname = " ";
1736 
1737   std::cout << " imeas " << imeas << std::endl;
1738   int Aline = 0;
1739   std::vector<Measurement*>::const_iterator vmcite;
1740   for (vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); ++vmcite) {
1741     for (ALIuint jj = 0; jj < ALIuint((*vmcite)->dim()); jj++) {
1742       if (Aline == imeas) {
1743         char ctmp[20];
1744         gcvt(jj, 10, ctmp);
1745         return ((*vmcite)->name()) + ":" + std::string(ctmp);
1746       }
1747       Aline++;
1748     }
1749   }
1750 
1751   std::cout << " return measname " << measname << std::endl;
1752   return measname;
1753 }