CMSSW Reference Manual

00001 #include <FWCore/Utilities/interface/Exception.h>
00002 //   COCOA class implementation file
00003 //Id:  Fit.cc
00004 //CAT: Fit
00005 //
00006 //   History: v1.0 
00007 //   Pedro Arce
00008 
00009 #include "Alignment/CocoaModel/interface/Model.h"
00010 #include "Alignment/CocoaModel/interface/OpticalObject.h"
00011 #include "Alignment/CocoaFit/interface/Fit.h"
00012 
00013 #include "Alignment/CocoaModel/interface/Measurement.h"
00014 #include "Alignment/CocoaModel/interface/Entry.h"
00015 #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
00016 #include "Alignment/CocoaUtilities/interface/ALIFileOut.h"
00017 #include "Alignment/CocoaUtilities/interface/GlobalOptionMgr.h"
00018 #include "Alignment/CocoaModel/interface/DeviationsFromFileSensor2D.h"
00019 #include "Alignment/CocoaAnalysis/interface/NtupleManager.h"
00020 #include "Alignment/CocoaAnalysis/interface/FittedEntriesManager.h"
00021 #include "Alignment/CocoaAnalysis/interface/FittedEntriesSet.h"
00022 #ifdef COCOA_VIS
00023 #include "Alignment/CocoaVisMgr/interface/ALIVRMLMgr.h"
00024 #include "Alignment/IgCocoaFileWriter/interface/IgCocoaFileMgr.h"
00025 #endif
00026 #include "Alignment/CocoaModel/interface/OpticalObjectMgr.h"
00027 #include "Alignment/CocoaModel/interface/ErrorCorrelationMgr.h"
00028 #include "Alignment/CocoaModel/interface/ErrorCorrelation.h"
00029 #include "Alignment/CocoaModel/interface/FittedEntriesReader.h"
00030 #include "Alignment/CocoaDaq/interface/CocoaDaqReader.h"
00031 #include <stdlib.h>
00032 #include <iomanip>
00033 #include <math.h>
00034 #include <time.h>
00035 #include <set>
00036  
00037 Fit* Fit::instance = 0;
00038 
00039 ALIMatrix* Fit::AMatrix;
00040 ALIMatrix* Fit::AtMatrix;
00041 ALIMatrix* Fit::WMatrix;
00042 ALIMatrix* Fit::AtWAMatrix;
00043 //op ALIMatrix* Fit::VaMatrix;
00044 ALIMatrix* Fit::DaMatrix;
00045 //op ALIMatrix* Fit::PDMatrix;
00046 //-ALIMatrix* Fit::VyMatrix;
00047 ALIMatrix* Fit::yfMatrix;
00048 //ALIMatrix* Fit::fMatrix;
00049 
00050 ALIint Fit::_NoLinesA;
00051 ALIint Fit::_NoColumnsA;
00052 //op ALIMatrix* Fit::thePropagationMatrix;
00053 
00054 ALIint Fit::theMinimumEntryQuality;
00055 ALIdouble Fit::thePreviousIterationFitQuality = DBL_MAX;
00056 ALIdouble Fit::theFitQualityCut = -1;
00057 ALIdouble Fit::theRelativeFitQualityCut = -1;
00058 ALIint Fit::theNoFitIterations;
00059 ALIint Fit::MaxNoFitIterations = -1;
00060 ALIdouble Fit::theMinDaFactor = 1.e-8;
00061 
00062 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00063 //@@  Gets the only instance of Model
00064 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00065 Fit& Fit::getInstance()
00066 {
00067   if(!instance) {
00068     instance = new Fit;
00069     ALIdouble go;
00070     GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
00071 
00072     gomgr->getGlobalOptionValue("maxDeviDerivative", go );
00073     ALIUtils::setMaximumDeviationDerivative( go );
00074     if( ALIUtils::debug >= 3 ) std::cout << " Fit::maximum_deviation_derivative " << ALIUtils::getMaximumDeviationDerivative() << std::endl;
00075 
00076     gomgr->getGlobalOptionValue("maxNoFitIterations", go );
00077     MaxNoFitIterations = int(go);
00078 
00079     gomgr->getGlobalOptionValue("fitQualityCut", go );
00080     theFitQualityCut = go;
00081     if( ALIUtils::debug >= 3 ) std::cout << " theFitQualityCut " << theFitQualityCut << std::endl;
00082 
00083     gomgr->getGlobalOptionValue("RelativeFitQualityCut", go );
00084     theRelativeFitQualityCut = go;
00085     if( ALIUtils::debug >= 3 ) std::cout << " theRelativeFitQualityCut " << theRelativeFitQualityCut << std::endl;
00086 
00087     gomgr->getGlobalOptionValue("minDaFactor", go );
00088     theMinDaFactor = go;
00089 
00090   }
00091 
00092   return *instance;
00093 }
00094 
00095 
00096 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00097 //@@  startFit: steering method to make the fit
00098 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00099 void Fit::startFit()
00100 { 
00101   //  Model::setCocoaStatus( COCOA_InitFit );
00102 
00103   NtupleManager* NTmgr = NtupleManager::getInstance();
00104   NTmgr->BookNtuple();
00105 
00106   ALIuint nEvent = 0;
00107   ALIUtils::setFirstTime( 1 );
00108 
00109   WriteVisualisationFiles();
00110 
00111   for(;;) {
00112 
00113     if( !fitNextEvent( nEvent ) ) break;
00114 
00115     //-    if ( ALIUtils::debug >= 0) std::cout << " FIT STATUS " << Model::printCocoaStatus( Model::getCocoaStatus() ) << std::endl;
00116 
00117   }
00118 
00119   //---------- Program ended, fill histograms of fitted entries
00120   GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
00121   if(gomgr->GlobalOptions()["histograms"] > 0) {
00122     FittedEntriesManager* FEmgr = FittedEntriesManager::getInstance();
00123     FEmgr->MakeHistos();
00124   }
00125 
00126   NTmgr->WriteNtuple();
00127 
00128 }
00129 
00130 
00131 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00132 ALIbool Fit::fitNextEvent( ALIuint& nEvent )
00133 {
00134   if( Model::getFittedEntriesReader() != 0 ) Model::getFittedEntriesReader()->readFittedEntriesFromFile();
00135 
00136   //----- Reset coordinates to those read at the start
00137   std::vector< OpticalObject* >::iterator voite;
00138   for( voite = Model::OptOList().begin(); voite !=  Model::OptOList().end(); voite++ ) {
00139     (*voite)->resetOriginalOriginalCoordinates();
00140   }
00141   
00142   //----- Reset entries displacements to 0.
00143   std::vector< Entry* >::iterator veite;
00144   for( veite = Model::EntryList().begin(); veite !=  Model::EntryList().end(); veite++ ) {    
00145     (*veite)->resetValueDisplacementByFitting();
00146   }
00147   
00148   
00149   ALIbool lastEvent = 0;
00150   
00151   //-    DeviationsFromFileSensor2D::setApply( 1 );
00152   
00153   //m  ALIbool moreDataSets = Model::readMeasurementsFromFile( Measurement::only1Date, Measurement::only1Time );
00154   
00155   //----- Check if there are more data sets
00156   ALIbool moreDataSets = 1;
00157   if(CocoaDaqReader::GetDaqReader() != 0) moreDataSets = CocoaDaqReader::GetDaqReader()->ReadNextEvent();
00158   
00159   if(ALIUtils::debug >= 2)  std::cout << CocoaDaqReader::GetDaqReader() << "$$$$$$$$$$$$$$$ moreData Sets " << moreDataSets << std::endl;
00160   
00161   if( moreDataSets ) {
00162     if( ALIUtils::debug >= 2 ) std::cout << "@@@@@@@@@@@@@@@@@@ Starting data set fit ..." << nEvent << std::endl;
00163 
00164     //----- Count entries to be fitted, and set their order in theFitPos
00165     setFittableEntries();
00166     
00167     //----- Dump dimensions of output in 'report.out' file
00168     ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
00169     fileout << std::endl << "@@@@@@@ NEW MEASUREMENT SET " << nEvent << std::endl;
00170     if( ALIUtils::report >= 1 ) ALIUtils::dumpDimensions( fileout );
00171     
00172     //----- reset Number of iterations of non linear fit
00173     theNoFitIterations = 0;
00174     
00175     GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
00176     ALIdouble dumpMat;
00177     gomgr->getGlobalOptionValue("save_matrices", dumpMat );
00178 
00179     //----- Fit parameters
00180     double daFactor = 1.;
00181     Model::setCocoaStatus(COCOA_FirstIterationInEvent );
00182     for(;; ){
00183     //---------- Calculate the original simulated values of each Measurement (when all entries have their read in values)
00184       calculateSimulatedMeasurementsWithOriginalValues(); //?? original changed atfer each iteration
00185    
00186       FitQuality fq = fitParameters( daFactor );
00187       if( dumpMat > 1 ) dumpMatrices();
00188 
00189       //-      evaluateFitQuality( fq, daFactor );
00190 
00191       //----- Check if new iteration must be done
00192       if( fq == FQsmallDistanceToMinimum ) {
00193         addDaMatrixToEntries();
00194         if(ALIUtils::report >= 1) dumpFittedValues( ALIFileOut::getInstance( Model::ReportFName() ), TRUE, TRUE );
00195         //--- Print entries in all ancestor frames
00196         ALIdouble go;   
00197         gomgr->getGlobalOptionValue("dumpInAllFrames", go );
00198         if(go >= 1) dumpFittedValuesInAllAncestorFrames( ALIFileOut::getInstance( Model::ReportFName() ), FALSE, FALSE );
00199 
00200         break;  // No more iterations
00201       } else if( fq == FQbigDistanceToMinimum ) {
00202         addDaMatrixToEntries();
00203         if(ALIUtils::report >= 1) dumpFittedValues( ALIFileOut::getInstance( Model::ReportFName() ), TRUE, TRUE );
00204 
00205         //----- Next iteration (if not too many already)
00206         theNoFitIterations++;
00207         daFactor = 1.;
00208 
00209         //----- Too many iterations: end event here
00210         if( theNoFitIterations >= MaxNoFitIterations ) {
00211           std::cerr << "!!!! WARNING: Too many iterations " << theNoFitIterations << "  and fit DOES NOT CONVERGE " << std::endl;
00212  
00213           if(ALIUtils::report >= 2) {
00214             ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
00215             fileout << "!!!! WARNING: Too many iterations " << theNoFitIterations << "  and fit DOES NOT CONVERGE " << std::endl;
00216           }
00217           //      Model::setCocoaStatus( COCOA_FitCannotImprove );
00218           break;  // No more iterations
00219         }
00220 
00221       } else if( fq == FQchiSquareWorsened ) {
00222         //----- Recalculate fit quality with decreasing values of Da
00223         //      std::cout << " quality daFactor " << daFactor << " " << theMinDaFactor << std::endl;
00224         daFactor *= 0.5;      
00225         if( daFactor > theMinDaFactor ){
00226           substractLastDisplacementToEntries( 0.5 );
00227 
00228           if(ALIUtils::report >= 2) {
00229             ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
00230             fileout << " Redoing iteration with Da factor " << daFactor << std::endl;
00231           }
00232         } else {
00233           daFactor *= 2.;
00234           std::cerr << " !!!ERROR: not possible to get good fit quality even multiplying Da by " << daFactor << std::endl;
00235           if(ALIUtils::report >= 2) {
00236             ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
00237             fileout << " !!!ERROR: not possible to get good fit quality even multiplying Da by " << daFactor << std::endl;
00238           }
00239           //        Model::setCocoaStatus( COCOA_FitCannotImprove );
00240           //-    std::cout << "fdsaf FIT STATUS " << Model::printCocoaStatus( Model::getCocoaStatus() ) << std::endl;
00241           break;  // No more iterations
00242         }
00243       }
00244       Model::setCocoaStatus(COCOA_NextIterationInEvent );
00245 
00246     }
00247     
00248     //----- Iteration is finished: dump fitted entries
00249     if(ALIUtils::debug >= 1) calculateSimulatedMeasurementsWithOriginalValues();
00250     if(gomgr->GlobalOptions()["histograms"] > 0) {
00251        FittedEntriesManager::getInstance()->AddFittedEntriesSet( new FittedEntriesSet( AtWAMatrix ) );
00252     }
00253     
00254     NtupleManager* ntupleMgr = NtupleManager::getInstance();
00255     ntupleMgr->InitNtuple();
00256     ntupleMgr->FillChi2();
00257     ntupleMgr->FillOptObjects(AtWAMatrix);
00258     ntupleMgr->FillMeasurements();
00259     ntupleMgr->FillFitParameters(AtWAMatrix);
00260     ntupleMgr->FillNtupleTree();
00261 
00262     //- only if not stopped in worsening quality state        if(ALIUtils::report >= 0) dumpFittedValues( ALIFileOut::getInstance( Model::ReportFName() ));
00263     
00264     /*-      std::vector< OpticalObject* >::iterator voite;
00265       for( voite = Model::OptOList().begin(); voite !=  Model::OptOList().end(); voite++ ) {
00266       //-??             (*voite)->resetOriginalOriginalCoordinates();
00267         }*/
00268     
00269     //---- If no measurement file, break after looping once
00270     //-      std::cout << " Measurement::measurementsFileName() " << Measurement::measurementsFileName() << " Measurement::measurementsFileName()" <<std::endl;
00271     if( CocoaDaqReader::GetDaqReader() == 0 ) {
00272       //m    if( Measurement::measurementsFileName() == "" ) {
00273       lastEvent = 1;
00274       return !lastEvent;
00275     }
00276     
00277     //-      std::cout << "  Measurement::only1" <<  Measurement::only1 << std::endl;
00278     if( Measurement::only1 ) {
00279       lastEvent = 1;
00280       return !lastEvent;
00281     }
00282     nEvent++;
00283   } else {
00284     lastEvent = 1;
00285     return !lastEvent;
00286   }
00287 
00288   return !lastEvent;
00289 }
00290 
00291 
00292 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00293 void Fit::WriteVisualisationFiles()
00294 {
00295 #ifdef COCOA_VIS
00296   if(gomgr->GlobalOptions()["VisOnly"] == 1) {
00297       calculateSimulatedMeasurementsWithOriginalValues(); //?? original changed atfer each iteration
00298   }
00299 
00300   GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
00301   if(gomgr->GlobalOptions()["VisWriteVRML"] > 0) {
00302     if(ALIUtils::getFirstTime()) ALIVRMLMgr::getInstance().writeFile();
00303   }
00304   if(gomgr->GlobalOptions()["VisWriteIguana"] > 0) {
00305     if(ALIUtils::getFirstTime()) IgCocoaFileMgr::getInstance().writeFile();
00306   }
00307   
00308   if(gomgr->GlobalOptions()["VisOnly"] == 1) {
00309     if(ALIUtils::debug >= 1 )std::cout << " Visualiation file(s) succesfully written. Ending.... " << std::endl;
00310     exit(1);
00311   }
00312 #endif    
00313 }
00314 
00315 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00316 //@@  Count how many entries are going to be fitted (have quality >=  theMinimumEntryQuality)
00317 //@@ Set for this entries the value of theFitPos
00318 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00319 void Fit::setFittableEntries() 
00320 {
00321 
00322   std::vector< Entry* >::const_iterator vecite;
00323 
00324   GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
00325   theMinimumEntryQuality = int(gomgr->GlobalOptions()[ALIstring("calcul_type")]) + 1;
00326   if ( ALIUtils::debug >= 4) std::cout << "Fit::setFittableEntries: total Entry List size= " << Model::EntryList().size() << std::endl;
00327 
00328   int No_entry_to_fit = 0;
00329   for ( vecite = Model::EntryList().begin();
00330        vecite != Model::EntryList().end(); vecite++ ) {  
00331 
00332     // Number the parameters that are going to be fitted
00333       if ( (*vecite)->quality() >= theMinimumEntryQuality ) {
00334           (*vecite)->setFitPos( No_entry_to_fit );
00335           if( ALIUtils::debug >= 4 ) std::cout << " Entry To Fit= " << No_entry_to_fit << " " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << "   with quality= " << (*vecite)->quality() << std::endl; 
00336           No_entry_to_fit++;
00337       }
00338   }
00339 
00340 }
00341 
00342 
00343 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00344 //@@ Main method in class Fit 
00345 //@@ fitParameters: get the parameters through the chi square fit
00346 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00347 FitQuality Fit::fitParameters( const double daFactor )
00348 {
00349   redoMatrices();
00350 
00351   if(ALIUtils::debug >= 2) std::cout << " Fit quality daFactor " << daFactor << std::endl;
00352 
00353   if(ALIUtils::debug >= 0) {
00354     std::cout << std::endl << "Fit iteration " << theNoFitIterations << " ..." << std::endl;
00355   }
00356 
00357 
00358   //---- Get chi2 of first iteration
00359   if( Model::getCocoaStatus() == COCOA_FirstIterationInEvent ) {
00360     thePreviousIterationFitQuality = GetSChi2( 0 );
00361 
00362     GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
00363     if (gomgr->GlobalOptions()[ ALIstring("stopAfter1stIteration") ] == 1) {
00364       std::cout << "@!! STOPPED by user after 1st iteration " << std::endl;
00365       exit(1);
00366     }
00367   }
00368 
00369   /*  //---------- Open output file
00370   if( ALIUtils::report >= 1 ) {
00371     ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
00372     //t    fileout << " REPORT.OUT " << std::endl;
00373     //t    ALIUtils::dumpDimensions( fileout );
00374     fileout << std::endl << "Fit iteration " << theNoFitIterations << " ..." << std::endl;
00375     }*/
00376  
00377   //-    std::cout << "2 FIT STATUS " << Model::printCocoaStatus( Model::getCocoaStatus() ) << std::endl;
00378 
00379   if(ALIUtils::debug >= 10) {
00380     std::cout << std::endl << " End fitParameters " << theNoFitIterations << " ..." << std::endl;
00381   }
00382 
00383   return getFitQuality();
00384 
00385 }
00386 
00387 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00388 void Fit::redoMatrices()
00389 {
00390   deleteMatrices(); 
00391 
00392   calculateSimulatedMeasurementsWithOriginalValues();
00393 
00394   PropagateErrors();
00395 
00396 }
00397 
00398 
00399 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00400 //@@  Propagate the Errors from the entries to the measurements
00401 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00402 //cocoaStatus
00403 void Fit::PropagateErrors()
00404 {
00405 
00406   //----- Create empty matrices of appropiate size
00407   CreateMatrices();
00408 
00409   //---- count running time
00410   time_t now;
00411   now = clock();
00412   if(ALIUtils::debug >= 2) std::cout << "TIME:CREATE_MAT    : " << now << " " << difftime(now, ALIUtils::time_now())/1.E6 << std::endl;
00413   ALIUtils::set_time_now(now); 
00414 
00415   //----- Fill the A, W & y matrices with the measurements
00416   FillMatricesWithMeasurements();
00417 
00418   //---- count running time
00419   now = clock();
00420   if(ALIUtils::debug >= 0) std::cout << "TIME:MAT_MEAS_FILLED: " << now << " " << difftime(now, ALIUtils::time_now())/1.E6 << std::endl;
00421   ALIUtils::set_time_now(now); 
00422 
00423   //----- Fill the A, W & y matrices with the calibrated parameters
00424   GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
00425   if (gomgr->GlobalOptions()[ ALIstring("calcul_type") ] == 0) {
00426     FillMatricesWithCalibratedParameters();
00427 
00428     //---- count running time
00429     now = clock();
00430     if(ALIUtils::debug >= 0) std::cout << "TIME:MAT_CAL_FILLED : " << now << " " << difftime(now, ALIUtils::time_now())/1.E6 << std::endl;
00431     ALIUtils::set_time_now(now); 
00432 
00433   }
00434 
00435   //----- Put by hand some correlations if known previously
00436   setCorrelationsInWMatrix();
00437 
00438   if(ALIUtils::debug >= 3) WMatrix->Dump("WMatrix before inverse");
00439 
00440   //----- Check first that matrix can be inverted
00441   if( m_norm1( WMatrix->MatNonConst() ) == 0 ) {
00442     //    Model::setCocoaStatus( COCOA_FitMatrixNonInversable );
00443     return; //  Model::getCocoaStatus();
00444   } else {
00445     WMatrix->inverse();
00446   }
00447   
00448   if(ALIUtils::debug >= 3) AMatrix->Dump("AMatrix");
00449   if(ALIUtils::debug >= 3) WMatrix->Dump("WMatrix");
00450   if(ALIUtils::debug >= 3) yfMatrix->Dump("yfMatrix");
00451 
00452   if(gomgr->GlobalOptions()["onlyDeriv"] >= 1) {
00453     std::cout << "ENDING after derivatives are calculated ('onlyDeriv' option set)" << std::endl;
00454     exit(1);
00455   }
00456 
00457   multiplyMatrices();
00458 
00459   now = clock();
00460   if(ALIUtils::debug >= 0) std::cout << "TIME:MAT_MULTIPLIED : " << now << " " << difftime(now, ALIUtils::time_now())/1.E6 << std::endl;
00461   ALIUtils::set_time_now(now); 
00462 
00463   if( ALIUtils::getFirstTime() == 1) ALIUtils::setFirstTime( 0 );
00464 
00465 
00466 }
00467 
00468 
00469 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00470 //@@ Calculate the simulated value of each Measurement propagating the LightRay when all the entries have their original values
00471 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00472 void Fit::calculateSimulatedMeasurementsWithOriginalValues()
00473 {
00474   //  if( ALIUtils::debug >= 4) OpticalObjectMgr::getInstance()->dumpOptOs();
00475 
00476   //---------- Set DeviationsFromFileSensor2D::apply true
00477   DeviationsFromFileSensor2D::setApply( 1 );
00478 
00479   if(ALIUtils::debug >= 2) std::cout << "Fit::calculateSimulatedMeasurementsWithOriginalValues" <<std::endl;
00480   //---------- Loop Measurements
00481   std::vector< Measurement* >::const_iterator vmcite;
00482   for ( vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); vmcite++) {
00483     //----- Calculate Simulated Value Original
00484     (*vmcite)->calculateOriginalSimulatedValue();
00485   }
00486 
00487   //---------- Set DeviationsFromFileSensor2D::apply false
00488   // 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
00489   DeviationsFromFileSensor2D::setApply( 0 );
00490 
00491 }
00492 
00493 
00494 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00495 void Fit::deleteMatrices()
00496 {
00497  //delete matrices created in previous iteration
00498   delete DaMatrix;
00499   delete AMatrix;
00500   delete WMatrix;
00501   delete yfMatrix;
00502   //op  delete fMatrix;
00503   delete AtMatrix;
00504   delete AtWAMatrix;
00505   //op  delete VaMatrix;
00506   //-  delete VyMatrix;
00507   //op  delete PDMatrix;
00508 
00509 }
00510 
00511 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00512 //@@  Calculate the NoLines & NoColumns and create matrices 
00513 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00514 void Fit::CreateMatrices()
00515 {
00516 
00517   //---------- Count number of measurements
00518   ALIint NoMeas = 0;
00519   std::vector< Measurement* >::const_iterator vmcite;
00520   for ( vmcite = Model::MeasurementList().begin();
00521     vmcite != Model::MeasurementList().end(); vmcite++ ) {
00522     NoMeas += (*vmcite)->dim();
00523   }
00524   if( ALIUtils::debug >= 9) std::cout << "NOMEAS" << NoMeas << std::endl;
00525 
00526    //-------- Count number of 'cal'ibrated parameters
00527   ALIint nEnt_cal = 0;
00528   ALIint noent = 0;
00529   //-  std::cout << Model::EntryList().size() << std::endl;
00530   GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
00531   if ( gomgr->GlobalOptions()[ "calcul_type" ] == 0) { // fit also 'cal' parameters
00532     //-  if( ALIUtils::debug >= 9) std::cout << "NOENTCALLL " << nEnt_cal << std::endl;
00533     if( ALIUtils::debug >= 5 ) std::cout << " Count number of 'cal'ibrated parameters " << std::endl;
00534     std::vector< Entry* >::iterator veite;
00535     for ( veite = Model::EntryList().begin();
00536           veite != Model::EntryList().end(); veite++ ) {
00537       if ( (*veite)->quality() == 1 ) nEnt_cal++; 
00538       noent++;
00539       if( ALIUtils::debug >= 6) {
00540         std::cout <<(*veite)->quality() << " " << (*veite)->OptOCurrent()->name() << " "
00541              << (*veite)->name() << " # ENT CAL " << nEnt_cal << " # ENT " << noent << std::endl; 
00542       }
00543     }
00544   }
00545 
00546   //---------- Count number parameters to be fitted ('cal' + 'unk')
00547   ALIint NoParamFit = 0;
00548   std::vector<Entry*>::const_iterator vecite;    
00549   for ( vecite = Model::EntryList().begin();  
00550         vecite != Model::EntryList().end(); vecite++) { 
00551     if ( (*vecite)->quality() >= theMinimumEntryQuality ) {
00552       NoParamFit++;
00553       if( ALIUtils::debug >= 99) std::cout <<(*vecite)->quality() << (*vecite)->OptOCurrent()->name() << (*vecite)->name() << "NoParamFit" << NoParamFit << std::endl;
00554       //      break;
00555     }
00556   }  
00557  
00558   //---------- Create Matrices
00559   ALIint NoLinesA = NoMeas + nEnt_cal;
00560   ALIint NoColumnsA = NoParamFit;
00561   AMatrix = new ALIMatrix( NoLinesA, NoColumnsA );
00562 
00563   ALIint NoLinesW = NoLinesA;
00564   ALIint NoColumnsW = NoLinesA;
00565   WMatrix = new ALIMatrix( NoLinesW, NoColumnsW );
00566 
00567   ALIint NoLinesY = NoLinesA;
00568   //op  yMatrix = new ALIMatrix( NoLinesY, 1 );
00569   yfMatrix = new ALIMatrix( NoLinesY, 1 );
00570 
00571   //op  fMatrix = new ALIMatrix( NoLinesY, 1 );
00572 
00573   if ( ALIUtils::debug >= 4 ) std::cout << "CreateMatrices: NoLinesA = " << NoLinesA <<
00574   " NoColumnsA = " << NoColumnsA << std::endl;
00575 }
00576 
00577 
00578 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00579 //@@  Loop Measurements:
00580 //@@    Fill Matrix A with derivatives respect to affecting entries 
00581 //@@    Fill Matrix W, y & f with values and sigmas of measurement coordinate
00582 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00583 void Fit::FillMatricesWithMeasurements() 
00584 {
00585 
00586   int Aline = 0; 
00587 
00588   //---------- Loop Measurements
00589   std::vector<Measurement*>::const_iterator vmcite;
00590   std::vector<Entry*>::const_iterator vecite;
00591   for ( vmcite = Model::MeasurementList().begin();
00592         vmcite != Model::MeasurementList().end(); vmcite++) {
00593     if( ALIUtils::debug >= 5 ) std::cout << "FillMatricesWithMeasurements: measurement " << (*vmcite)->name() << " # entries affecting " <<(*vmcite)->affectingEntryList().size() << std::endl;
00594 
00595     //-------- Array of derivatives with respect to entries
00596     ALIint measdim = (*vmcite)->dim(); 
00597     ALIdouble* derivRE;
00598     derivRE = new ALIdouble[measdim];
00599 
00600     //-------- Fill matrix A:
00601     //------ Loop only Entries affecting this Measurement
00602     //-std::cout << "number affecting entries: " << (*vmcite)->affectingEntryList().size() << std::endl;
00603     for ( vecite = (*vmcite)->affectingEntryList().begin();  
00604           vecite != (*vmcite)->affectingEntryList().end(); vecite++) { 
00605       //-------- If good quality, get derivative of measurement with respect to this Entry
00606       if ( (*vecite)->quality() >= theMinimumEntryQuality ) {
00607         if ( ALIUtils::debug >= 4) {
00608           //      std::cout << "FillMatricesWithMeasurements: filling element ( " << Aline << " - " << Aline+measdim-1 << " , " << (*vecite)->fitPos() << std::endl;
00609           std::cout <<"entry affecting: " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() <<std::endl;
00610         }
00611         derivRE = (*vmcite)->DerivativeRespectEntry(*vecite);
00612         //---------- Fill matrix A with derivatives
00613         for ( ALIuint jj = 0; jj < ALIuint(measdim); jj++) {
00614           AMatrix->AddData( Aline+jj, (*vecite)->fitPos(), derivRE[jj] );
00615           if( ALIUtils::debug >= 5) std::cout << "FillMatricesWithMeasurements: AMATRIX (" << Aline+jj << "," << (*vecite)->fitPos() << " = " << derivRE[jj] << std::endl;
00616           //---------- Reset Measurement simulated_value
00617           (*vmcite)->setValueSimulated( jj, (*vmcite)->valueSimulated_orig(jj) );         
00618         }
00619       }
00620     }
00621     delete[] derivRE;
00622     
00623     //---------- Fill matrices W, y and f:
00624     //------ Loop Measurement coordinates
00625     for ( ALIuint jj=0; jj < ALIuint((*vmcite)->dim()); jj++) {
00626       ALIdouble sigma = (*vmcite)->sigma()[jj];
00627       if ( sigma == 0. ) {
00628         std::cerr << "EXITING: Measurement number " << 
00629           vmcite - Model::MeasurementList().begin() << 
00630           "has 0 error!!" << std::endl;
00631       } else {
00632         //----- Fill W Matrix with inverse of sigma squared
00633         // multiply error by cameraScaleFactor
00634         ALIdouble sigmanew = sigma * Measurement::cameraScaleFactor;
00635         //      std::cout << Aline+jj << " WMATRIX FILLING " << sigmanew << " * " << Measurement::cameraScaleFactor << std::endl;
00636         WMatrix->AddData( Aline+jj, Aline+jj, (sigmanew*sigmanew) );
00637       }
00638       //op //----- Fill Matrices y with measurement value 
00639       //op yMatrix->AddData( Aline+jj, 0, (*vmcite)->value()[jj] );
00640       //op //----- Fill f Matrix with simulated_value
00641       //op fMatrix->AddData( Aline+jj, 0, (*vmcite)->valueSimulated_orig(jj) );
00642       //----- Fill Matrix y - f with measurement value - simulated value
00643       yfMatrix->AddData( Aline+jj, 0, (*vmcite)->value()[jj] - (*vmcite)->valueSimulated_orig(jj) );
00644       //      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;
00645     }
00646     if ( ALIUtils::debug >= 99) std::cout << "change line" << Aline << std::endl;
00647     Aline += measdim;
00648     if ( ALIUtils::debug >= 99) std::cout << "change line" << Aline << std::endl;
00649     
00650   }
00651 
00652   if ( ALIUtils::debug >= 4) AMatrix->Dump("Matrix A with meas");
00653   if ( ALIUtils::debug >= 4) WMatrix->Dump("Matrix W with meas");
00654   if ( ALIUtils::debug >= 4) yfMatrix->Dump("Matrix y with meas");
00655 
00656 }
00657 
00658 
00659 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00660 //@@  Loop Measurements:
00661 //@@    Fill Matrix A with derivatives respect to affecting entries 
00662 //@@    Fill Matrix W, y & f with values and sigmas of measurement coordinate
00663 //@@ 
00664 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00665 void Fit::FillMatricesWithCalibratedParameters() 
00666 {
00667   //---------- Count how many measurements 
00668   ALIint NolinMes = 0;
00669   std::vector<Measurement*>::const_iterator vmcite;
00670   for ( vmcite = Model::MeasurementList().begin();
00671        vmcite != Model::MeasurementList().end(); vmcite++) {
00672       NolinMes += (*vmcite)->dim();
00673   }
00674   if(ALIUtils::debug>=4) std::cout << "@@FillMatricesWithCalibratedParameters" << std::endl;
00675 
00676   std::vector< Entry* >::const_iterator vecite;
00677   ALIint nEntcal = 0;
00678   //---------- Loop entries 
00679   for ( vecite = Model::EntryList().begin();
00680         vecite != Model::EntryList().end(); vecite++ ) {
00681 //                  (= No parameters to be fitted - No parameters 'unk' )
00682     //-    std::cout << "entry" << (*veite) << std::endl;
00683     //----- Take entries of quality = 'cal' 
00684     if ( (*vecite)->quality() == 1 ){
00685       //--- Matrix A: fill diagonals with 1. (derivatives of entry w.r.t itself)
00686       ALIint lineNo = NolinMes + nEntcal;  
00687       ALIint columnNo = (*vecite)->fitPos();  //=? nEntcal
00688       AMatrix->AddData( lineNo, columnNo, 1. );
00689       if(ALIUtils::debug >= 4) std::cout << "Fit::FillMatricesWithCalibratedParameters:  AMatrix ( " << lineNo << " , " << columnNo  << ") = " << 1. << std::endl;
00690 
00691       //--- Matrix W: sigma*sigma
00692       ALIdouble entsig = (*vecite)->sigma();
00693       if(ALIUtils::debug >= 4) std::cout << "Fit::FillMatricesWithCalibratedParameters:  WMatrix ( " << lineNo << " , " << columnNo  << ") = " << entsig*entsig << std::endl;
00694       WMatrix->AddData( lineNo, lineNo, entsig*entsig );
00695 
00696       //--- Matrix y & f: fill it with 0. 
00697       //op      yMatrix->AddData( lineNo, 0, (*vecite)->value());
00698       //op      yfMatrix->AddData( lineNo, 0, (*vecite)->value());
00699       ALIdouble calFit;
00700       GlobalOptionMgr* gomgr = GlobalOptionMgr::getInstance();
00701       gomgr->getGlobalOptionValue("calParamInyfMatrix", calFit );
00702       if( calFit ) {
00703         yfMatrix->AddData( lineNo, 0, -(*vecite)->valueDisplacementByFitting() );
00704         //-     yfMatrix->AddData( lineNo, 0, (*vecite)->value() );
00705         //-             yfMatrix->AddData( lineNo, 0, (*vecite)->lastAdditionToValueDisplacementByFitting() );
00706         //-     ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
00707         //      fileout << "cal to yf " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " " << (*vecite)->valueDisplacementByFitting() << endl;
00708         //      std::cout << "call to yf " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " " << (*vecite)->valueDisplacementByFitting() << std::endl;
00709 
00710       } else {
00711         yfMatrix->AddData( lineNo, 0, 0. );
00712       }
00713       //t      if(ALIUtils::debug >= 5) std::cout << "Fit::FillMatricesWithCalibratedParameters:  yfMatrix ( " << lineNo << " , " << columnNo  << ") = " << (*yfMatrix)(lineNo)(0) << std::endl;
00714       nEntcal++;
00715     }
00716   }
00717   
00718 }
00719 
00720 
00721 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00722 //@@  Gets the only instance of Model
00723 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00724 void Fit::setCorrelationsInWMatrix()
00725 {  
00726   //----- Check if there are correlations to input
00727   ErrorCorrelationMgr* corrMgr = ErrorCorrelationMgr::getInstance();
00728   ALIint siz = corrMgr->getNumberOfCorrelations();
00729   if( siz == 0 ) return;
00730 
00731   //----- Set correlations
00732   ALIuint ii;
00733   for( ii = 0; ii < ALIuint(siz); ii++ ){
00734   //t    if(ALIUtils::debug >= 5) std::cout << "globaloption cmslink fit" << Model::GlobalOptions()["cms_link"] << std::endl;
00735     ErrorCorrelation* corr = corrMgr->getCorrelation( ii );
00736     setCorrelationFromParamFitted( corr->getEntry1(), corr->getEntry2(), corr->getCorrelation() );
00737   }
00738 
00739 }
00740 
00741 
00742 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00743 //@@  set correlation between two entries of two OptOs
00744 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00745 void Fit::setCorrelationFromParamFitted( const pss& entry1, const pss& entry2,
00746  ALIdouble correl ) 
00747 {
00748 
00749   ALIint pmsize = WMatrix->NoLines();
00750   ALIint fit_pos1 = Model::getEntryByName(entry1.first, entry1.second)->fitPos();
00751   ALIint fit_pos2 = Model::getEntryByName(entry2.first, entry2.second)->fitPos();
00752   std::cout << "CHECKsetCorrelatiFPF " << fit_pos1 << " " << fit_pos2 << std::endl;
00753 
00754   if( fit_pos1 >= 0 && fit_pos1 < pmsize && fit_pos2 >= 0 && fit_pos2 < pmsize ) {
00755     setCorrelationFromParamFitted( fit_pos1, fit_pos2, correl );
00756   }
00757 }
00758 
00759 
00760 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00761 void Fit::setCorrelationFromParamFitted( const ALIint fit_pos1, const ALIint fit_pos2, ALIdouble correl ) 
00762 {
00763   //  ALIdouble error1 = sqrt( (*WMatrix)(fit_pos1, fit_pos1) );
00764   // ALIdouble error2 = sqrt( (*WMatrix)(fit_pos2, fit_pos2) );
00765   WMatrix->SetCorrelation( fit_pos1, fit_pos2, correl );
00766   std::cout << "setCorrelatiFPF " << fit_pos1 << " " << fit_pos2 << " " << correl << std::endl;
00767 }
00768 
00769 
00770 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00771 //@@  multiply matrices needed for fit
00772 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00773 void Fit::multiplyMatrices() 
00774 {
00775   if(ALIUtils::debug >= 4) std::cout << "@@Multiplying matrices " << std::endl;
00776   //---------- Calculate transpose of A
00777   AtMatrix = new ALIMatrix( *AMatrix );
00778   if(ALIUtils::debug >= 5) AtMatrix->Dump("AtMatrix=A");
00779   //-  std::cout << "call transpose";
00780   AtMatrix->transpose();
00781   if(ALIUtils::debug >= 4) AtMatrix->Dump("AtMatrix");
00782 
00783   //---------- Calculate At * W * A
00784   AtWAMatrix = new ALIMatrix(0, 0);   
00785   //  if(ALIUtils::debug >= 5) AtWAMatrix->Dump("AtWAMatrix=0");
00786   *AtWAMatrix = *AtMatrix * *WMatrix * *AMatrix;   
00787   if(ALIUtils::debug >= 5) AtWAMatrix->Dump("AtWAMatrix");
00788 
00789   CheckIfFitPossible();
00790 
00791   //t  AtWAMatrix->EliminateLines(0,48);
00792   //t AtWAMatrix->EliminateColumns(0,48);
00793   time_t now;
00794   now = clock();
00795   if(ALIUtils::debug >= 0) std::cout << "TIME:BEFORE_INVERSE : " << now << " " << difftime(now, ALIUtils::time_now())/1.E6 << std::endl;
00796   ALIUtils::set_time_now(now); 
00797 
00798   /*  std::cout << " DETERMINANT W " <<  m_norm1( AtWAMatrix->MatNonConst() ) << std::endl;
00799   if( m_norm1( AtWAMatrix->MatNonConst() ) == 0 ) {
00800     std::cout << " DETERMINANT W " <<  m_norm1( AtWAMatrix->MatNonConst() ) << std::endl;
00801     std::exception();
00802     } */
00803 
00804   AtWAMatrix->inverse();
00805   if(ALIUtils::debug >= 4) AtWAMatrix->Dump("inverse AtWAmatrix");
00806   now = clock();
00807   if(ALIUtils::debug >= 0) std::cout << "TIME:AFTER_INVERSE  : " << now << " " << difftime(now, ALIUtils::time_now())/1.E6 << std::endl;
00808   ALIUtils::set_time_now(now); 
00809 
00810   //op  thePropagationMatrix = AtWAMatrix;
00811 
00812   //op  VaMatrix = new ALIMatrix( *AtWAMatrix );
00813   
00814   //----- Print out propagated errors of parameters (=AtWA diagonal elements)
00815   std::vector< Entry* >::const_iterator vecite;
00816  
00817   if( ALIUtils::debug >= 4 ) {
00818     std::cout << "PARAM" << "        Optical Object " << "   entry name " << "      Param.Value " 
00819          << " Prog.Error" << " Orig.Error" << std::endl;
00820   }
00821 
00822   ALIint nEnt = 0;
00823   ALIint nEntUnk = 0;
00824   for ( vecite = Model::EntryList().begin();
00825     vecite != Model::EntryList().end(); vecite++ ) {
00826 //------------------ Number of parameters 'cal' 
00827 //                  (= No parameters to be fitted - No parameters 'unk' )
00828     if( (*vecite)->quality() >= theMinimumEntryQuality ){
00829       if( ALIUtils::debug >= 4) {
00830         std::cout << nEnt << "PARAM" << std::setw(26) 
00831              << (*vecite)->OptOCurrent()->name().c_str() 
00832              << std::setw(8) << " " << (*vecite)->name().c_str() << " " 
00833              << std::setw(8) << " " << (*vecite)->value() << " " 
00834              << std::setw(8) << sqrt(AtWAMatrix->Mat()->me[nEnt][nEnt]) / 
00835                 (*vecite)->OutputSigmaDimensionFactor() 
00836              << " " << (*vecite)->sigma() / (*vecite)->OutputSigmaDimensionFactor() 
00837              << " Q" << (*vecite)->quality() << std::endl;
00838       }
00839       nEnt++;
00840     }
00841     if ( (*vecite)->quality() == 2 ) nEntUnk++;
00842   }
00843 
00844   if(ALIUtils::debug >= 5) yfMatrix->Dump("PD(y-f)Matrix final");  
00845 
00846 }
00847 
00848 
00849 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00850 //@@ check also that the fit_quality = SMat(0,0) is smaller for each new iteration
00851 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00852 FitQuality Fit::getFitQuality( const ALIbool canBeGood ) 
00853 {
00854 
00855   double fit_quality = GetSChi2(1);
00856 
00857   //---------- Calculate DS = Variable to recognize convergence (distance to minimum)
00858   ALIMatrix* DatMatrix = new ALIMatrix( *DaMatrix );
00859   //  delete DaMatrix; //op
00860   DatMatrix->transpose();
00861   if(ALIUtils::debug >= 5) DatMatrix->Dump("DatMatrix");
00862   //op  ALIMatrix* DSMat = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix * *PDMatrix);
00863   ALIMatrix* DSMat = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix * *yfMatrix);
00864   ALIMatrix* DSMattemp = new ALIMatrix(*DatMatrix * *AtMatrix * *WMatrix);
00865   if(ALIUtils::debug >= 5) DSMattemp->Dump("DSMattempMatrix=Dat*At*W");
00866   ALIMatrix* DSMattemp2 = new ALIMatrix(*AtMatrix * *WMatrix * *yfMatrix);
00867   if(ALIUtils::debug >= 5) DSMattemp2->Dump("DSMattempMatrix2=At*W*yf");
00868   ALIMatrix* DSMattemp3 = new ALIMatrix(*AtMatrix * *WMatrix);
00869   if(ALIUtils::debug >= 5) DSMattemp3->Dump("DSMattempMatrix3=At*W");
00870   if(ALIUtils::debug >= 5) AtMatrix->Dump("AtMatrix");
00871   /*  for( int ii = 0; ii < DatMatrix->NoColumns(); ii++ ){
00872     std::cout << ii << " DS term " << (*DatMatrix)(0,ii) * (*DSMattemp2)(ii,0) << std::endl;
00873     }*/
00874   //  delete AtMatrix; //op
00875   //  delete WMatrix; //op
00876 
00877   //op  if(ALIUtils::debug >= 5) (*PDMatrix).Dump("PDMatrix");
00878   if(ALIUtils::debug >= 5) (*yfMatrix).Dump("yfMatrix");
00879   if(ALIUtils::debug >= 5) DSMat->Dump("DSMatrix final");
00880   //  delete yfMatrix; //op
00881 
00882   ALIdouble fit_quality_cut = (*DSMat)(0,0);  
00883   //-  ALIdouble fit_quality_cut =fabs( (*DSMat)(0,0) );  
00884   delete DSMat;
00885   if(ALIUtils::debug >= 0) std::cout << theNoFitIterations << " Fit quality predicted improvement in distance to minimum is = " << fit_quality_cut << std::endl;
00886   if( ALIUtils::report >= 2 ) {
00887     ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
00888     fileout << theNoFitIterations << " Fit quality predicted improvement in distance to minimum is = " << fit_quality_cut << std::endl;
00889   }
00890 
00891   //-  double fit_quality_cut = thePreviousIterationFitQuality - fit_quality;
00892   //- double fit_quality_cut = fit_quality;
00893   
00894   //----- Check quality 
00895   time_t now;
00896   now = clock();
00897   if(ALIUtils::debug >= 0) std::cout << "TIME:QUALITY_CHECKED: " << now << " " << difftime(now, ALIUtils::time_now())/1.E6 << std::endl;
00898   ALIUtils::set_time_now(now); 
00899 
00900   FitQuality fitQuality;
00901 
00902   //----- Chi2 is bigger, bad
00903   //    if( theNoFitIterations != 0 && fit_quality_cut > 0. ) {
00904   if( fit_quality_cut < 0. ) {
00905     fitQuality = FQchiSquareWorsened;
00906    if(ALIUtils::debug >= 0) std::cerr << "!!WARNING: Fit quality has worsened: Fit Quality now = " << fit_quality
00907               << " before " << thePreviousIterationFitQuality << " diff " << fit_quality - thePreviousIterationFitQuality << std::endl;
00908 
00909   //----- Chi2 is smaller, check if we make another iteration    
00910   } else {
00911     //----- Small chi2 change: end
00912     if( (fit_quality_cut < theFitQualityCut || fabs(thePreviousIterationFitQuality - fit_quality_cut )/fit_quality_cut < theRelativeFitQualityCut ) && canBeGood ) {
00913       fitQuality = FQsmallDistanceToMinimum;
00914       if(ALIUtils::report >= 1) {
00915         ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
00916         fileout << "STOP: SMALL IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut << " < " << theFitQualityCut << std::endl;
00917       }
00918       if(ALIUtils::debug >= 4) {
00919         std::cout << "STOP: SMALL IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut << " < " << theFitQualityCut << std::endl;
00920       }
00921       
00922       //----- Big chi2 change: go to next iteration
00923     } else {
00924       fitQuality = FQbigDistanceToMinimum;
00925       //----- set thePreviousIterationFitQuality for next iteration 
00926       thePreviousIterationFitQuality = fit_quality;     
00927       
00928       if(ALIUtils::report >= 2) {
00929         ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
00930         fileout << "CONTINUE: BIG IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut << " >= " << theFitQualityCut << std::endl;
00931       }
00932       if(ALIUtils::debug >= 4) {
00933         std::cout << "CONTINUE: BIG IMPROVEMENT IN ITERATION " << theNoFitIterations << " = " << fit_quality_cut << " >= " << theFitQualityCut << std::endl;
00934       } 
00935     }
00936   }
00937   
00938   return fitQuality;
00939 
00940 }
00941 
00942 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00943 ALIdouble Fit::GetSChi2( ALIbool useDa )
00944 {
00945   ALIMatrix* SMat = 0;
00946   if( useDa ){
00947     //----- Calculate variables to check quality of this set of parameters
00948     
00949     //----- Calculate Da = (At * W * A)-1 * At * W * (y-f)
00950     /*t  DaMatrix = new ALIMatrix( *AtWAMatrix );
00951      *DaMatrix *= *AtMatrix * *WMatrix;
00952      if(ALIUtils::debug >= 5) DaMatrix->Dump("DaMatrix before yf ");
00953      *DaMatrix *= *yfMatrix;
00954      if(ALIUtils::debug >= 5) DaMatrix->Dump("DaMatrix");
00955     */
00956     
00957     DaMatrix = new ALIMatrix(0, 0);   
00958     //  if(ALIUtils::debug >= 5) AtWAMatrix->Dump("AtWAMatrix=0");
00959     *DaMatrix = ( *AtWAMatrix * *AtMatrix * *WMatrix * *yfMatrix);
00960     if(ALIUtils::debug >= 5) DaMatrix->Dump("DaMatrix");
00961     
00962     //----- Calculate S = chi2 = Fit quality = r^T W r (r = residual = f + A*Da - y )
00963     //op  ALIMatrix* tmpM = new ALIMatrix( *AMatrix * *DaMatrix + *PDMatrix );
00964     //  ALIMatrix* tmpM = new ALIMatrix( *AMatrix * *DaMatrix + *yfMatrix );
00965     
00966     ALIMatrix* tmpM = new ALIMatrix( 0,0 );
00967     *tmpM  = *AMatrix * *DaMatrix - *yfMatrix;    
00968     if(ALIUtils::debug >= 5) tmpM->Dump("A*Da + f - y Matrix ");
00969 
00970     ALIMatrix* tmptM = new ALIMatrix( *tmpM );
00971     tmptM->transpose();
00972     if(ALIUtils::debug >= 5) tmptM->Dump("tmptM after transpose");
00973     if(ALIUtils::debug >= 5) WMatrix->Dump("WMatrix");
00974 
00975     //  std::cout << "smat " << std::endl;
00976     //o  ALIMatrix* SMat = new ALIMatrix(*tmptM * *WMatrix * *tmpM);
00977     ALIMatrix* SMat1 = MatrixByMatrix(*tmptM,*WMatrix);
00978     //  ALIMatrix* SMat1 = MatrixByMatrix(*AMatrix,*WMatrix);
00979     if(ALIUtils::debug >= 5) SMat1->Dump("(A*Da + f - y)^T * W  Matrix");
00980     SMat = MatrixByMatrix(*SMat1,*tmpM);
00981     //  std::cout << "smatc " << std::endl;
00982     delete tmpM;
00983     delete tmptM;
00984     if(ALIUtils::debug >= 5) SMat->Dump("SMatrix with Da");
00985   } else {
00986     ALIMatrix* yftMat = new ALIMatrix(*yfMatrix);
00987     yftMat->transpose();
00988     SMat = new ALIMatrix(*yftMat *  *WMatrix * *yfMatrix);
00989     delete yftMat;
00990     if(ALIUtils::debug >= 5) SMat->Dump("SMatrix no Da");
00991   }
00992   ALIdouble fit_quality = (*SMat)(0,0);
00993   delete SMat;
00994   if(ALIUtils::debug >= 5) std::cout << " GetSChi2 " << useDa << " = " << fit_quality << std::endl;
00995 
00996   PrintChi2( fit_quality, !useDa );
00997 
00998   return fit_quality;
00999   
01000 }
01001 
01002 
01003 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01004 //@@ Correct entries with fitted values  
01005 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01006 void Fit::addDaMatrixToEntries() 
01007 {
01008 
01009   if(ALIUtils::debug >= 4) {
01010     std::cout << "@@ Adding correction (DaMatrix) to Entries " << std::endl;
01011     DaMatrix->Dump("DaMatrix =");
01012   }
01013 
01014   /*-  Now there are other places where entries are changed
01015     if( ALIUtils::report >= 3 ) {
01016     ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
01017     fileout << std::endl << " CHANGE IN ENTRIES" << std::endl
01018             << "          Optical Object       Parameter   xxxxxx " << std::endl;
01019             }*/
01020 
01021   ALIint nEnt = 0;
01022   std::vector<Entry*>::const_iterator vecite; 
01023   for ( vecite = Model::EntryList().begin();
01024     vecite != Model::EntryList().end(); vecite++ ) {
01025 //------------------ Number of parameters 'cal' 
01026 //                  (= No parameters to be fitted - No parameters 'unk' )
01027     //-   std::cout << "qual" << (*vecite)->quality() << theMinimumEntryQuality << std::endl;
01028     if ( (*vecite)->quality() >= theMinimumEntryQuality ){
01029       if ( ALIUtils::debug >= 5) {
01030         std::cout << std::endl << " @@@ PENTRY change " 
01031                   << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " " 
01032                   << " change= " << (*DaMatrix)(nEnt,0)
01033                   << " value= " << (*vecite)->valueDisplacementByFitting()
01034                   << std::endl;
01035       }
01036       /*      if( ALIUtils::report >=3 ) {
01037         ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
01038         fileout << "dd" << std::setw(30) << (*vecite)->OptOCurrent()->name() 
01039                 << std::setw(8) << " " << (*vecite)->name() << " " 
01040                 << std::setw(8) << " " << (*DaMatrix)(nEnt,0) / (*vecite)->OutputValueDimensionFactor()
01041                 << " " << (*vecite)->valueDisplacementByFitting() / (*vecite)->OutputValueDimensionFactor() << " fitpos " << (*vecite)->fitPos()
01042                 << std::endl;
01043                 }*/
01044 
01045       //----- Store this displacement 
01046       if(ALIUtils::debug >= 4) std::cout << " old valueDisplacementByFitting " << (*vecite)->name() << " " << (*vecite)->valueDisplacementByFitting() << " original value " <<  (*vecite)->value() <<std::endl; 
01047 
01048       (*vecite)->addFittedDisplacementToValue( (*DaMatrix)(nEnt,0) );
01049 
01050       if(ALIUtils::debug >= 4) std::cout << nEnt << " new valueDisplacementByFitting " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " = " << (*vecite)->valueDisplacementByFitting() << " " << (*DaMatrix)(nEnt,0) << std::endl ; 
01051       nEnt++;
01052     }
01053   }
01054 
01055 }
01056 
01057 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01058 //@@ Delete the previous addition of fitted values (to try a new one daFactor times smaller )
01059 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01060 void Fit::substractLastDisplacementToEntries( const ALIdouble factor ) 
01061 {
01062 
01063   if(ALIUtils::debug >= 4) {
01064     std::cout << "@@  Fit::substractToHalfDaMatrixToEntries " << std::endl;
01065   }
01066 
01067   std::vector<Entry*>::const_iterator vecite; 
01068   for ( vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); vecite++ ) {
01069     if ( (*vecite)->quality() >= theMinimumEntryQuality ){
01070       //--     (*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
01071 
01072       ALIdouble lastadd = (*vecite)->lastAdditionToValueDisplacementByFitting() * factor;
01073       //-      if( lastadd < 0 ) lastadd *= -1;
01074       (*vecite)->addFittedDisplacementToValue( -lastadd );
01075       (*vecite)->setLastAdditionToValueDisplacementByFitting( - (*vecite)->lastAdditionToValueDisplacementByFitting() );
01076       //      (*vecite)->substractToHalfFittedDisplacementToValue();
01077 
01078         if(ALIUtils::debug >= 4) std::cout << " new valueDisplacementByFitting " << (*vecite)->OptOCurrent()->name() << " " << (*vecite)->name() << " = " << (*vecite)->valueDisplacementByFitting() << " " << std::endl ; 
01079     }
01080   }
01081 
01082 }
01083 
01084 
01085 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01086 //@@  Dump all the entries that have been fitted (those that were 'cal' or 'unk'
01087 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01088 void Fit::dumpFittedValues( ALIFileOut& fileout, ALIbool printErrors, ALIbool printOrig )
01089 {
01090   //---------- print
01091   if(ALIUtils::debug >= 0) {
01092     std::cout << "SRPARPOS " << "               Optical Object  " 
01093          << "      Parameter" <<  " Fit.Value " << " Orig.Value" << std::endl;
01094   }
01095   //---------- Dump header
01096   if(ALIUtils::debug >= 0) std::cout << std::endl << "FITTED VALUES " << std::endl; 
01097   fileout << std::endl << "FITTED VALUES " << std::endl 
01098           << "nEnt_unk"
01099           << "             Optical Object" 
01100           << "        Parameter  ";
01101   if( printErrors ) {
01102     fileout << " value (+-error)" 
01103             << " orig.val (+-error)";
01104   } else {
01105     fileout << " value " 
01106             << " orig.val ";
01107   }
01108   fileout << " quality" 
01109           << std::endl; 
01110 
01111   //---------- Iterate over OptO list
01112   std::vector< Entry* > entries;
01113   //  const Entry* entry;
01114   int ii, siz;
01115   std::vector< OpticalObject* >::const_iterator vocite;
01116   for( vocite = Model::OptOList().begin(); vocite != Model::OptOList().end(); vocite++ ) {
01117     if( (*vocite)->type() == ALIstring("system") ) continue;
01118 
01119     fileout << " %%%% Optical Object: " << (*vocite)->longName() << std::endl;
01120 
01121     entries = (*vocite)->CoordinateEntryList();
01122     siz = entries.size();
01123     if( siz != 6 ) {
01124       std::cerr << "!!! FATAL ERROR: strange number of coordinates = " << siz << std::endl;
01125       abort();
01126     }
01127 
01128     //----- 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)
01129     OpticalObject* opto = entries[0]->OptOCurrent();
01130     const OpticalObject* optoParent = opto->parent();
01131     printCentreInOptOFrame( opto, optoParent, fileout, printErrors, printOrig );
01132 
01133     //----- Dump entry angles coordinates
01134     printRotationAnglesInOptOFrame( opto, optoParent, fileout, printErrors, printOrig );
01135 
01136     //----- Dump extra entries
01137     entries = (*vocite)->ExtraEntryList();
01138     siz = entries.size();
01139     for( ii = 0; ii < siz; ii++ ){
01140       double entryvalue = getEntryValue( entries[ii] );
01141       dumpEntryAfterFit( fileout, entries[ii], entryvalue, printErrors, printOrig );
01142     }
01143   }
01144 
01145   dumpEntryCorrelations( fileout );
01146 
01147 }
01148 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01149 //@@  Dump all the entries that have been fitted in reference frames of all ancestors
01150 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01151 void Fit::dumpFittedValuesInAllAncestorFrames( ALIFileOut& fileout, ALIbool printErrors, ALIbool printOrig )
01152 {
01153   //---------- print
01154   fileout << std::endl << "@@@@ FITTED VALUES IN ALL ANCESTORS " << std::endl 
01155           << "nEnt_unk"
01156           << "             Optical Object" 
01157           << "        Parameter  ";
01158   if( printErrors ) {
01159     fileout << " value (+-error)"; 
01160       if( printOrig ){
01161         fileout << " orig.val (+-error)";
01162       }
01163   } else {
01164     fileout << " value ";
01165       if( printOrig ){
01166         fileout << " orig.val ";
01167       }
01168   }
01169   fileout << " quality" 
01170           << std::endl; 
01171 
01172   //---------- Iterate over OptO list
01173   std::vector< Entry* > entries;
01174   std::vector< OpticalObject* >::const_iterator vocite;
01175   for( vocite = Model::OptOList().begin(); vocite != Model::OptOList().end(); vocite++ ) {
01176     if( (*vocite)->type() == ALIstring("system") ) continue;
01177 
01178     fileout << " %%%% Optical Object: " << (*vocite)->longName() << std::endl;
01179 
01180     entries = (*vocite)->CoordinateEntryList();
01181 
01182     //----- 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)
01183     OpticalObject* opto = *vocite;
01184     const OpticalObject* optoParent = opto->parent();
01185     do {
01186       fileout << " %% IN FRAME : " << optoParent->longName() << std::endl;
01187       printCentreInOptOFrame( opto, optoParent, fileout, printErrors, printOrig );
01188 
01189       //----- Dump entry angles coordinates
01190       printRotationAnglesInOptOFrame( opto, optoParent, fileout, printErrors, printOrig );
01191       optoParent = optoParent->parent();
01192     }while( optoParent );
01193   }
01194 
01195 }
01196 
01197 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01198 void Fit::printCentreInOptOFrame( const OpticalObject* opto, const OpticalObject* optoAncestor, ALIFileOut& fileout, ALIbool printErrors, ALIbool printOrig )
01199 {
01200   Hep3Vector centreLocal;
01201   if( optoAncestor->type() == "system" ) {
01202     centreLocal = opto->centreGlob();
01203   } else {
01204     centreLocal = opto->centreGlob() - optoAncestor->centreGlob();
01205     CLHEP::HepRotation parentRmGlobInv = inverseOf( optoAncestor->rmGlob() );
01206     centreLocal = parentRmGlobInv * centreLocal;
01207   }
01208   if(ALIUtils::debug >= 2 ) {
01209     std::cout << "CENTRE LOCAL "<< opto->name() << " " << centreLocal << " GLOBL " << opto->centreGlob() << " parent GLOB " << optoAncestor->centreGlob() << std::endl;
01210     ALIUtils::dumprm( optoAncestor->rmGlob(), " parent rm " );
01211   }
01212   std::vector< Entry* > entries = opto->CoordinateEntryList();  
01213   for( uint ii = 0; ii < 3; ii++ ){
01214     /* double entryvalue = getEntryValue( entries[ii] );
01215        ALIdouble entryvalue;
01216        if( ii == 0 ) {
01217        entryvalue = centreLocal.x();
01218        }else if( ii == 1 ) {
01219        entryvalue = centreLocal.y();
01220        }else if( ii == 2 ) {
01221        entryvalue = centreLocal.z();
01222        }*/
01223     dumpEntryAfterFit( fileout, entries[ii], centreLocal[ii] / entries[ii]->OutputValueDimensionFactor(), printErrors, printOrig );
01224   }
01225 
01226 }
01227 
01228 
01229 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01230 void Fit::printRotationAnglesInOptOFrame( const OpticalObject* opto, const OpticalObject* optoAncestor, ALIFileOut& fileout, ALIbool printErrors, ALIbool printOrig )
01231 {
01232   std::vector< Entry* > entries = opto->CoordinateEntryList();  
01233   std::vector<double> entryvalues = opto->getRotationAnglesInOptOFrame( optoAncestor, entries );
01234   //-    std::cout << " after return entryvalues[0] " << entryvalues[0] << " entryvalues[1] " << entryvalues[1] << " entryvalues[2] " << entryvalues[2] << std::endl;
01235   for( uint ii = 3; ii < entries.size(); ii++ ){
01236     dumpEntryAfterFit( fileout, entries[ii], entryvalues[ii-3]/entries[ii]->OutputValueDimensionFactor(), printErrors, printOrig );
01237   }
01238 
01239 }
01240 
01241 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01242 double Fit::getEntryValue( const Entry* entry )
01243 {
01244   double entryvalue;
01245   if( entry->type() == "angles") {
01246     if(ALIUtils::debug >= 2 ) std::cout << "WARNING valueDisplacementByFitting has no sense for angles "  << std::endl;
01247     entryvalue = -999;
01248   }
01249   entryvalue = ( entry->value() + entry->valueDisplacementByFitting() ) / entry->OutputValueDimensionFactor();
01250   return entryvalue;
01251 }
01252 
01253 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01254 void Fit::dumpEntryAfterFit( ALIFileOut& fileout, const Entry* entry, double entryvalue, ALIbool printErrors, ALIbool printOrig )
01255 {
01256   //-  std::cout << " Fit::dumpEntryAfterFit " << entryvalue << std::endl;
01257   ALIdouble dimv = entry->OutputValueDimensionFactor();
01258   ALIdouble dims = entry->OutputSigmaDimensionFactor();
01259   //----- Dump to std::cout
01260   if(ALIUtils::debug >= 3) {
01261     std::cout << "ENTRY: " 
01262          << std::setw(30) << entry->OptOCurrent()->name()
01263          << std::setw(8) << " " << entry->name() << " " 
01264          << std::setw(8) << ( entry->value() + entry->valueDisplacementByFitting() ) 
01265          <<" " << entry->value() 
01266          << " Q" << entry->quality() << std::endl;
01267   }
01268   
01269   if ( entry->quality() == 2 ) {
01270     fileout << "UNK: " << entry->fitPos() << " ";   
01271   } else if ( entry->quality() == 1 ) {
01272     fileout << "CAL: " << entry->fitPos() << " "; 
01273     //    fileout << "CAL: -1 "; 
01274   } else {
01275     fileout << "FIX: -1 ";
01276   }
01277   
01278   fileout << std::setw(30)  << entry->OptOCurrent()->name()
01279           << std::setw(8) << " " << entry->name() << " " 
01280           << std::setw(8) << std::setprecision(8) << entryvalue;
01281   if ( entry->quality() >= theMinimumEntryQuality ) {
01282     if( printErrors ) fileout << " +- " << std::setw(8) << sqrt(AtWAMatrix->Mat()->me[entry->fitPos()][entry->fitPos()]) / dims;
01283   } else { 
01284     if( printErrors ) fileout << " +- " << std::setw(8) << 0.;
01285   }
01286   if( printOrig ) {
01287     fileout << std::setw(8) << " " << entry->value() / dimv;
01288     if( printErrors ) fileout << " +- " << std::setw(8) << entry->sigma() /dims << " Q" << entry->quality();
01289 
01290     if( ALIUtils::report >= 2) {
01291       float dif = ( entry->value() + entry->valueDisplacementByFitting() ) / dimv - entry->value() / dimv;
01292       if( fabs(dif) < 1.E-9 ) dif = 0.;
01293       fileout << " DIFF= " << dif;
01294       // << " == " << ( entry->value() + entry->valueDisplacementByFitting() )  / dimv - entryvalue << " @@ " << ( entry->value() + entry->valueDisplacementByFitting() ) / dimv << " @@ " <<  entryvalue;
01295     } else {
01296       //        fileout << std::endl;
01297     }
01298   }
01299 
01300   fileout << std::endl;
01301   
01302 }
01303 
01304 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01305 void Fit::dumpEntryCorrelations( ALIFileOut& fileout )
01306 {
01307   //----- Only dump correlations bigger than a factor
01308   ALIdouble minCorrel = 1.E-6;
01309   //----- Dump correlations
01310   fileout << std::endl << "CORRELATION BETWEEN 'unk' ENTRIES: (>= " << minCorrel<< " )" << std::endl
01311           << "No_1  No_2   correlation " << std::endl;
01312 
01313   ALIuint i1,i2;
01314   std::vector< Entry* >::iterator veite1, veite2;
01315   std::string E1, E2;
01316   for ( veite1 = Model::EntryList().begin();
01317         veite1 != Model::EntryList().end(); veite1++ ) {
01318     if( (*veite1)->quality() == 0 ) {
01319       continue;
01320     } else if( (*veite1)->quality() == 1 ) {
01321       E1 = "C";
01322     } else if( (*veite1)->quality() == 2 ) {
01323       E1 = "U";
01324     }
01325     i1 = (*veite1)->fitPos();
01326 
01327     for ( veite2 = veite1+1;
01328         veite2 != Model::EntryList().end(); veite2++ ) {
01329       i2 = (*veite2)->fitPos(); 
01330       if( (*veite2)->quality() == 0 ) {
01331         continue;
01332       } else if( (*veite2)->quality() == 1 ) {
01333         E2 = "C";
01334       } else if( (*veite2)->quality() == 2 ) {
01335         E2 = "U";
01336       }
01337       ALIdouble corr = AtWAMatrix->Mat()->me[i1][i2];
01338       ALIdouble corrf = corr / sqrt(AtWAMatrix->Mat()->me[i1][i1])
01339           / sqrt(AtWAMatrix->Mat()->me[i2][i2]);
01340       if (fabs(corrf) >= minCorrel ) {
01341         if(ALIUtils::debug >= 0) {
01342           std::cout  << "CORR:" << E1 << "" << E2 << " (" << i1 << ")" <<  " (" << i2 << ")" << " " << corrf << std::endl;
01343         }
01344         fileout  << "CORR:" << E1 << "" << E2 << " (" << i1 << ")" <<  " (" << i2 << ")" << " " << corrf << std::endl;
01345       }
01346     }
01347   }
01348   //------- Dump optical object list 
01349   if( ALIUtils::debug >= 2) OpticalObjectMgr::getInstance()->dumpOptOs();
01350  
01351 }
01352 
01353 
01354 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01355 //@@ Dump matrices used for the fit
01356 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01357 void Fit::dumpMatrices() 
01358 {
01359   //----- Dump matrices for this iteration
01360   ALIFileOut& matout = ALIFileOut::getInstance( Model::MatricesFName() );
01361   //  ofstream matout("matrices.out");
01362   matout << std::endl << " @@@@@@@@@@@@@@@  Iteration No : " << theNoFitIterations << std::endl;
01363   AMatrix->ostrDump( matout, "Matrix A" );
01364   AtMatrix->ostrDump( matout, "Matrix At" );
01365   WMatrix->ostrDump( matout, "Matrix W" );
01366   AtWAMatrix->ostrDump( matout, "Matrix AtWA" );
01367   //op  VaMatrix->ostrDump( matout, "Matrix Va" );
01368   DaMatrix->ostrDump( matout, "Matrix Da" );
01369   yfMatrix->ostrDump( matout, "Matrix y" );
01370   //op  fMatrix->ostrDump( matout, "Matrix f" );
01371   //op  thePropagationMatrix->ostrDump( matout, "propagation Matrix " );
01372   matout.close();
01373 
01374 }
01375 
01376 
01377 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01378 //@@ findEntryFitPosition 
01379 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01380 ALIint Fit::findEntryFitPosition( const ALIstring& opto_name, const ALIstring& entry_name ) 
01381 {
01382   ALIint fitposi = -99;
01383 
01384   OpticalObject* opto = Model::getOptOByName( opto_name );
01385   //-  std::cout << "OPTO = " << opto->name() << std::endl;
01386   std::vector<Entry*>::const_iterator vecite;
01387   for (vecite = opto->CoordinateEntryList().begin(); 
01388        vecite < opto->CoordinateEntryList().end(); vecite++) {
01389     //-    std::cout << "ENTRYLIST" << (*vecite)->name() << std::endl;
01390     if ((*vecite)->name() == entry_name ) {
01391       //-  std::cout << "FOUND " << std::endl;
01392       fitposi = (*vecite)->fitPos();
01393     }
01394   }
01395   for (vecite = opto->ExtraEntryList().begin(); 
01396        vecite < opto->ExtraEntryList().end(); vecite++) {
01397     //-    std::cout << "ENTRYLIST" << (*vecite)->name() << std::endl;
01398     if ((*vecite)->name() == entry_name ) {
01399       //- std::cout << "FOUND " << std::endl;
01400       fitposi = (*vecite)->fitPos();
01401     }
01402   }
01403 
01404   if(fitposi == -99) {
01405     std::cerr << "!!EXITING: entry name not found: " << entry_name << std::endl;
01406     exit(2);
01407   } else {
01408     return fitposi;
01409   }
01410 }
01411 
01412 
01413 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01414 void Fit::PrintChi2( ALIdouble fit_quality, ALIbool isFirst )
01415 {
01416   double chi2meas = 0; 
01417   double chi2cal = 0;
01418   ALIint nMeas = 0, nUnk = 0;
01419 
01420   //----- Calculate the chi2 of measurements
01421   std::vector< Measurement* >::const_iterator vmcite;
01422   for ( vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); vmcite++) {
01423     //--- Calculate Simulated Value Original
01424     for ( ALIuint ii = 0; ii < ALIuint((*vmcite)->dim()); ii++ ){
01425       nMeas++;
01426       double c2 = ( (*vmcite)->value(ii) - (*vmcite)->valueSimulated(ii) ) / (*vmcite)->sigma(ii);
01427       chi2meas += c2*c2; 
01428       if( ALIUtils::debug >= 0) {
01429         std::cout << c2 << " adding chi2meas "  << chi2meas << " " << (*vmcite)->name() << ": " << ii << " (mm)R: " << (*vmcite)->value(ii)*1000. << " S: " << (*vmcite)->valueSimulated(ii)*1000. << " Diff= " << ((*vmcite)->value(ii) - (*vmcite)->valueSimulated(ii))*1000. << std::endl;
01430       }
01431     }
01432   }
01433 
01434   //----- Calculate the chi2 of calibrated parameters
01435   std::vector< Entry* >::iterator veite;
01436   for ( veite = Model::EntryList().begin();
01437         veite != Model::EntryList().end(); veite++ ) {
01438     if ( (*veite)->quality() == 2 ) nUnk++;
01439     if ( (*veite)->quality() == 1 ) {
01440       double c2 = (*veite)->valueDisplacementByFitting() / (*veite)->sigma();
01441       //double c2 = (*veite)->value() / (*veite)->sigma();
01442       chi2cal += c2*c2;
01443       if( ALIUtils::debug >= 0) std::cout << c2 << " adding chi2cal "  << chi2cal << " " << (*veite)->OptOCurrent()->name() << " " << (*veite)->name() << std::endl;
01444       //-       std::cout << " valueDisplacementByFitting " << (*veite)->valueDisplacementByFitting() << " sigma " << (*veite)->sigma() << std::endl;
01445     }
01446   }
01447   
01448   if( ALIUtils::report >= 1) {
01449     ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
01450     fileout << " Chi2= " << chi2meas+chi2cal << " / " << nMeas-nUnk << " dof " << "  From measurements= " << chi2meas << " from calibrated parameters= " << chi2cal << std::endl;
01451   }
01452   if( ALIUtils::debug >= 3) std::cout << " quality Chi2 (no correlations) " << chi2meas+chi2cal << " " << chi2meas << " " << chi2cal << std::endl;
01453 
01454 
01455   if( !isFirst ) {
01456     //    double fit_quality_change = thePreviousIterationFitQuality - fit_quality;
01457     
01458     if(ALIUtils::debug >= 0) { 
01459       std::cout << std::endl << "Fit iteration " << theNoFitIterations << " ..." << std::endl;
01460       //      std::cout << theNoFitIterations << " Chi2 improvement in this iteration = " << fit_quality_change << std::endl;
01461     }
01462     if( ALIUtils::report >= 1 ) {
01463       ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
01464       fileout << std::endl << "Fit iteration " << theNoFitIterations << " ..." << std::endl;
01465       //      fileout << theNoFitIterations << " Chi2 improvement in this iteration = " << fit_quality_change << std::endl;
01466     }
01467   }
01468 
01469   //---- Print chi2
01470   if(ALIUtils::debug >= 0) std::cout << theNoFitIterations << " Chi2 after iteration = " << fit_quality << std::endl;
01471   if( ALIUtils::report >= 1 ) {
01472     //--------- Get report file handler
01473     ALIFileOut& fileout = ALIFileOut::getInstance( Model::ReportFName() );
01474     fileout << theNoFitIterations << " Chi2 after iteration = " << fit_quality << std::endl;
01475   }
01476   
01477 }
01478 
01479 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01480 void Fit::CheckIfFitPossible()
01481 {   
01482   //----- Check if there is an unknown parameter that is not affecting any measurement
01483   ALIint NolinMes = 0;
01484   std::vector<Measurement*>::const_iterator vmcite;
01485   for ( vmcite = Model::MeasurementList().begin();
01486         vmcite != Model::MeasurementList().end(); vmcite++) {
01487     NolinMes += (*vmcite)->dim();
01488   }
01489   
01490   std::vector< Entry* >::const_iterator vecite;
01491   for ( vecite = Model::EntryList().begin(); vecite != Model::EntryList().end(); vecite++ ) {
01492     if( ALIUtils::debug >= 3 ) std::cout << "Fit::CheckFitIsPossible looping for entry " << (*vecite)->longName() << std::endl;
01493     if ( (*vecite)->quality() == 2 ) {
01494       ALIint nCol =  (*vecite)->fitPos();
01495       //--- Check all measurements
01496       ALIbool noDepend = TRUE;
01497       if( ALIUtils::debug >= 3 ) std::cout << "Fit::CheckFitIsPossible looping for entry " << nCol << std::endl;
01498       for( ALIint ii = 0; ii < NolinMes; ii++ ) {
01499         if( ALIUtils::debug >= 4 ) std::cout << ii << "Derivative " <<  (*AMatrix)(ii,nCol)  << std::endl;
01500 
01501         if( fabs((*AMatrix)(ii,nCol)) > ALI_DBL_MIN ) {
01502           if( ALIUtils::debug >= 3 ) std::cout << "Fit::CheckIfFitIsPossible " << nCol << " " << ii << " = " << (*AMatrix)(ii,nCol) << std::endl;
01503           noDepend = FALSE;
01504           break;
01505         }
01506       }
01507       if( noDepend ){
01508         throw cms::Exception("SDFError") << "!!!FATAL ERROR: Fit::CheckIfFitPossible() no measurement depends on unknown entry " << (*vecite)->OptOCurrent()->name() << "/" << (*vecite)->name() << std::endl
01509                   << "!!! Fit will not be possible! " << std::endl;
01510       }
01511     }
01512   }
01513   
01514   //------ 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!
01515 
01516   std::vector< Entry* >::const_iterator vecite1,vecite2;
01517   ALIint nLin = AMatrix->NoLines();
01518   ALIdouble derivPrec = ALI_DBL_MIN;
01519   //---------- Loop entries 
01520   for ( vecite1 = Model::EntryList().begin();
01521         vecite1 != Model::EntryList().end(); vecite1++ ) {
01522     if( (*vecite1)->quality() == 2 ) {
01523       vecite2 = vecite1; vecite2++;
01524       for ( ;vecite2 != Model::EntryList().end(); vecite2++ ) {
01525         if( (*vecite2)->quality() == 2 ) { 
01526           ALIint fitpos1 =  (*vecite1)->fitPos();
01527           ALIint fitpos2 =  (*vecite2)->fitPos();
01528           if( ALIUtils::debug >= 5 ) std::cout << "Fit::CheckIfFitIsPossible checking " << (*vecite1)->longName() << " ( " << fitpos1 << " ) & " << (*vecite2)->longName() << " ( " << fitpos2 << " ) "<< std::endl;
01529           ALIdouble prop = DBL_MAX;
01530           ALIbool isProp = TRUE;
01531           for( ALIint ii = 0; ii < nLin; ii++ ) {
01532             if( ALIUtils::debug >= 5 ) std::cout << "Fit::CheckIfFitIsPossible " << ii << " : " << (*AMatrix)(ii,fitpos1) << " ?= " << (*AMatrix)(ii,fitpos2) << std::endl;
01533             if( fabs((*AMatrix)(ii,fitpos1)) < derivPrec ) {
01534               if( fabs((*AMatrix)(ii,fitpos2)) > derivPrec ) {
01535                 isProp = FALSE;
01536                 break;
01537               }
01538             } else {
01539               ALIdouble propn = (*AMatrix)(ii,fitpos2) / (*AMatrix)(ii,fitpos1);
01540               if( prop != DBL_MAX && prop != propn ) {
01541                 isProp = FALSE;
01542                 break;
01543               }
01544               prop = propn;
01545             }
01546           }
01547           if( isProp ) {
01548               std::cerr << "!!!FATAL ERROR: Fit::CheckIfFitPossible() two entries for which the measurements have the same dependency (or 100% proportional) " <<  (*vecite1)->longName() << " & " <<  (*vecite2)->longName() << std::endl
01549                         << "!!! Fit will not be possible! " << std::endl;
01550             throw cms::Exception("SDFError");
01551           }
01552         }
01553       }
01554     }
01555   }
01556 }
01557 
01558 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01559 int Fit::CheckIfMeasIsProportionalToAnother( uint measNo )
01560 {
01561   int measProp = -1;
01562 
01563   std::set<uint> columnsEqual;
01564   std::set<uint> columnsEqualSave;
01565   uint biggestColumn = 0;
01566   ALIdouble biggest = 0.;
01567   for (int ii = 0; ii < AMatrix->NoColumns(); ii++ ){
01568     if( fabs((*AMatrix)(measNo,ii)) > biggest ) {
01569       biggest = fabs((*AMatrix)(measNo,ii));
01570       biggestColumn = ii;
01571     }
01572     columnsEqualSave.insert(ii); 
01573   }
01574   
01575   ALIdouble div;
01576 
01577   for( int jj = 0; jj < AMatrix->NoLines(); jj++ ){
01578     if( jj == int(measNo) ) continue;
01579     columnsEqual = columnsEqualSave;
01580     // check if ratio of each column to 'biggestColumn' is the same as for the N measurement
01581     for (int ii = 0; ii < AMatrix->NoColumns(); ii++ ){
01582       div = (*AMatrix)(measNo,ii)/(*AMatrix)(measNo,biggestColumn);
01583       if( fabs((*AMatrix)(jj,ii))  > ALI_DBL_MIN && fabs(div - (*AMatrix)(jj,ii)/(*AMatrix)(jj,biggestColumn) ) > ALI_DBL_MIN ) {
01584         if( ALIUtils::debug >= 3 ) std::cout << "CheckIfMeasIsProportionalToAnother 2 columns = " << ii << " in " << measNo << " & " << jj << std::endl;
01585       } else {
01586         if( ALIUtils::debug >= 3 ) std::cout << "CheckIfMeasIsProportionalToAnother 2 columns != " << ii << " in " << measNo << " & " << jj << std::endl;
01587         // if it is not equal delete this column 
01588         std::set<uint>::iterator ite = columnsEqual.find( ii );
01589         if( ite != columnsEqual.end() ){
01590           columnsEqual.erase(ite);
01591         }
01592       }  
01593     }
01594     // check if not all columns have been deleted from columnsEqual
01595     if( columnsEqual.size() != 0 ) {
01596       if( ALIUtils::debug >= 3 ) std::cout << "CheckIfMeasIsProportionalToAnother " <<  measNo << " = " << jj << std::endl;
01597       measProp = jj;
01598       break;
01599     }
01600   }
01601 
01602   return measProp;
01603 }
01604 
01605 
01606 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
01607 std::string Fit::GetMeasurementName( int imeas )
01608 {
01609   std::string measname = " ";
01610   
01611   std::cout << " imeas " << imeas << std::endl;
01612   int Aline = 0;  
01613   std::vector< Measurement* >::const_iterator vmcite;
01614   for ( vmcite = Model::MeasurementList().begin(); vmcite != Model::MeasurementList().end(); vmcite++) {
01615     for ( ALIuint jj = 0; jj < ALIuint((*vmcite)->dim()); jj++) {
01616       if( Aline == imeas ) {
01617         char ctmp[20];
01618         gcvt( jj, 10, ctmp );
01619         return ((*vmcite)->name()) + ":" + std::string(ctmp);
01620       }
01621       Aline++;
01622     }
01623   }
01624 
01625   std::cout << " return measname " << measname << std::endl;
01626   return measname;
01627 }
01628