Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:55:59 

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   for (vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); ++vecite) {
0891     //------------------ Number of parameters 'cal'
0892     //                  (= No parameters to be fitted - No parameters 'unk' )
0893     if ((*vecite)->quality() >= theMinimumEntryQuality) {
0894       if (ALIUtils::debug >= 4) {
0895         std::cout << nEnt << "PARAM" << std::setw(26) << (*vecite)->OptOCurrent()->name().c_str() << std::setw(8) << " "
0896                   << (*vecite)->name().c_str() << " " << std::setw(8) << " " << (*vecite)->value() << " "
0897                   << std::setw(8) << sqrt(AtWAMatrix->Mat()->me[nEnt][nEnt]) / (*vecite)->OutputSigmaDimensionFactor()
0898                   << " " << (*vecite)->sigma() / (*vecite)->OutputSigmaDimensionFactor() << " Q" << (*vecite)->quality()
0899                   << std::endl;
0900       }
0901       nEnt++;
0902     }
0903   }
0904 
0905   if (ALIUtils::debug >= 5)
0906     yfMatrix->Dump("PD(y-f)Matrix final");
0907 }
0908 
0909 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0910 //@@ check also that the fit_quality = SMat(0,0) is smaller for each new iteration
0911 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0912 FitQuality Fit::getFitQuality(const ALIbool canBeGood) {
0913   if (ALIUtils::debug >= 3)
0914     std::cout << "@@@ Fit::getFitQuality" << std::endl;
0915 
0916   double fit_quality = GetSChi2(true);
0917 
0918   //---------- Calculate DS = Variable to recognize convergence (distance to minimum)
0919   ALIMatrix* DatMatrix = new ALIMatrix(*DaMatrix);
0920   //  delete DaMatrix; //op
0921   DatMatrix->transpose();
0922   if (ALIUtils::debug >= 5)
0923     DatMatrix->Dump("DatMatrix");
0924   //op  ALIMatrix* DSMat = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix * *PDMatrix);
0925   ALIMatrix* DSMat = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix * *yfMatrix);
0926   if (ALIUtils::debug >= 5) {
0927     ALIMatrix* DSMattemp = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix);
0928     DSMattemp->Dump("DSMattempMatrix=Dat*At*W");
0929     ALIMatrix* DSMattemp2 = new ALIMatrix(*AtMatrix * *WMatrix * *yfMatrix);
0930     DSMattemp2->Dump("DSMattempMatrix2=At*W*yf");
0931     ALIMatrix* DSMattemp3 = new ALIMatrix(*AtMatrix * *WMatrix);
0932     DSMattemp3->Dump("DSMattempMatrix3=At*W");
0933     AtMatrix->Dump("AtMatrix");
0934   }
0935 
0936   /*  for( int ii = 0; ii < DatMatrix->NoColumns(); ii++ ){
0937     std::cout << ii << " DS term " << (*DatMatrix)(0,ii) * (*DSMattemp2)(ii,0) << std::endl;
0938     }*/
0939   //  delete AtMatrix; //op
0940   //  delete WMatrix; //op
0941 
0942   //op  if(ALIUtils::debug >= 5) (*PDMatrix).Dump("PDMatrix");
0943   if (ALIUtils::debug >= 5)
0944     (*yfMatrix).Dump("yfMatrix");
0945   if (ALIUtils::debug >= 5)
0946     DSMat->Dump("DSMatrix final");
0947   //  delete yfMatrix; //op
0948 
0949   ALIdouble fit_quality_cut = (*DSMat)(0, 0);
0950   //-  ALIdouble fit_quality_cut =std::abs( (*DSMat)(0,0) );
0951   delete DSMat;
0952   if (ALIUtils::debug >= 0)
0953     std::cout << theNoFitIterations
0954               << " Fit quality predicted improvement in distance to minimum is = " << fit_quality_cut << std::endl;
0955   if (ALIUtils::report >= 2) {
0956     ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
0957     fileout << theNoFitIterations
0958             << " Fit quality predicted improvement in distance to minimum is = " << fit_quality_cut << std::endl;
0959   }
0960 
0961   //-  double fit_quality_cut = thePreviousIterationFitQuality - fit_quality;
0962   //- double fit_quality_cut = fit_quality;
0963   //-  std::cout << "  fit_quality_cut " <<  fit_quality_cut << " fit_quality " << fit_quality << std::endl;
0964 
0965   //----- Check quality
0966   time_t now;
0967   now = clock();
0968   if (ALIUtils::debug >= 0)
0969     std::cout << "TIME:QUALITY_CHECKED: " << now << " " << difftime(now, ALIUtils::time_now()) / 1.E6 << std::endl;
0970   ALIUtils::set_time_now(now);
0971 
0972   FitQuality fitQuality;
0973 
0974   //----- Chi2 is bigger, bad
0975   //    if( theNoFitIterations != 0 && fit_quality_cut > 0. ) {
0976   if (fit_quality_cut < 0.) {
0977     fitQuality = FQchiSquareWorsened;
0978     if (ALIUtils::debug >= 1)
0979       std::cerr << "!!WARNING: Fit quality has worsened: Fit Quality now = " << fit_quality << " before "
0980                 << thePreviousIterationFitQuality << " diff " << fit_quality - thePreviousIterationFitQuality
0981                 << std::endl;
0982 
0983     //----- Chi2 is smaller, check if we make another iteration
0984   } else {
0985     ALIdouble rel_fit_quality = std::abs(thePreviousIterationFitQuality - fit_quality) / fit_quality;
0986     //----- Small chi2 change: end
0987     if ((fit_quality_cut < theFitQualityCut || rel_fit_quality < theRelativeFitQualityCut) && canBeGood) {
0988       if (ALIUtils::debug >= 2)
0989         std::cout << "$$ Fit::getFitQuality good " << fit_quality_cut << " <? " << theFitQualityCut << " || "
0990                   << rel_fit_quality << " <? " << theRelativeFitQualityCut << " GOOD " << canBeGood << std::endl;
0991       fitQuality = FQsmallDistanceToMinimum;
0992       if (ALIUtils::report >= 1) {
0993         ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
0994         fileout << "STOP: SMALL IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut << " < "
0995                 << theFitQualityCut << " OR (RELATIVE) " << rel_fit_quality << " < " << theRelativeFitQualityCut
0996                 << std::endl;
0997       }
0998       if (ALIUtils::debug >= 4) {
0999         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) {
1014         ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
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 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1031 ALIdouble Fit::GetSChi2(ALIbool useDa) {
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 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1102 void Fit::addDaMatrixToEntries() {
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 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1156 void Fit::substractLastDisplacementToEntries(const ALIdouble factor) {
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 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1296 void Fit::printCentreInOptOFrame(const OpticalObject* opto,
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        }*/
1325     dumpEntryAfterFit(
1326         fileout, entries[ii], centreLocal[ii] / entries[ii]->OutputValueDimensionFactor(), printErrors, printOrig);
1327   }
1328 }
1329 
1330 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1331 void Fit::printRotationAnglesInOptOFrame(const OpticalObject* opto,
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++) {
1340     dumpEntryAfterFit(
1341         fileout, entries[ii], entryvalues[ii - 3] / entries[ii]->OutputValueDimensionFactor(), printErrors, printOrig);
1342   }
1343 }
1344 
1345 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1346 double Fit::getEntryValue(const Entry* entry) {
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 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1364 void Fit::dumpEntryAfterFit(
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 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1414 void Fit::dumpEntryCorrelations(ALIFileOut& fileout) {
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)
1460     OpticalObjectMgr::getInstance()->dumpOptOs();
1461 }
1462 
1463 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1464 //@@ Dump matrices used for the fit
1465 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1466 void Fit::dumpMatrices() {
1467   //----- Dump matrices for this iteration
1468   ALIFileOut& matout = ALIFileOut::getInstance(Model::MatricesFName());
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 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1516 void Fit::PrintChi2(ALIdouble fit_quality, ALIbool isFirst) {
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) {
1554     ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
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) {
1570       ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
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
1581     ALIFileOut& fileout = ALIFileOut::getInstance(Model::ReportFName());
1582     fileout << theNoFitIterations << " Chi2 after iteration = " << fit_quality << std::endl;
1583   }
1584 }
1585 
1586 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1587 void Fit::CheckIfFitPossible() {
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 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1680 int Fit::CheckIfMeasIsProportionalToAnother(ALIuint measNo) {
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 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
1733 std::string Fit::GetMeasurementName(int imeas) {
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!