MillePede2.h

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#ifndef ALICEO2_FWDALIGN_MILLEPEDE2_H
#define ALICEO2_FWDALIGN_MILLEPEDE2_H
 
#include <vector>
#include <TString.h>
#include <TTree.h>
#include "ForwardAlign/MinResSolve.h"
#include "ForwardAlign/MillePedeRecord.h"
#include "ForwardAlign/SymMatrix.h"
#include "ForwardAlign/RectMatrix.h"
#include "ForwardAlign/MatrixSparse.h"
#include "ForwardAlign/MatrixSq.h"
#include "ForwardAlign/MilleRecordWriter.h"
#include "ForwardAlign/MilleRecordReader.h"
 
class TFile;
class TStopwatch;
class TArrayL;
class TArrayF;
 
namespace o2
{
namespace fwdalign
{
 
class MillePede2
{
 public:
  //
  enum { kFailed,
         kInvert,
         kNoInversion };   // used global matrix solution methods
  enum { kFixParID = -1 }; // dummy id for fixed param
 
  MillePede2();
  MillePede2(const MillePede2& src);
  virtual ~MillePede2();
  MillePede2& operator=(const MillePede2&)
  {
    printf("Dummy\n");
    return *this;
  }
 
  int InitMille(int nGlo, const int nLoc,
                const int lNStdDev = -1, const double lResCut = -1.,
                const double lResCutInit = -1., const std::vector<int>& regroup = {});
 
  int GetNGloPar() const { return fNGloPar; }
  int GetNGloParIni() const { return fNGloParIni; }
  std::vector<int> GetRegrouping() const { return fkReGroup; }
  int GetNLocPar() const { return fNLocPar; }
  long GetNLocalEquations() const { return fNLocEquations; }
  int GetCurrentIteration() const { return fIter; }
  int GetNMaxIterations() const { return fMaxIter; }
  int GetNStdDev() const { return fNStdDev; }
  int GetNGlobalConstraints() const { return fNGloConstraints; }
  int GetNLagrangeConstraints() const { return fNLagrangeConstraints; }
  int GetNLocalFits() const { return fNLocFits; }
  long GetNLocalFitsRejected() const { return fNLocFitsRejected; }
  int GetNGlobalsFixed() const { return fNGloFix; }
  int GetGlobalSolveStatus() const { return fGloSolveStatus; }
  float GetChi2CutFactor() const { return fChi2CutFactor; }
  float GetChi2CutRef() const { return fChi2CutRef; }
  float GetResCurInit() const { return fResCutInit; }
  float GetResCut() const { return fResCut; }
  int GetMinPntValid() const { return fMinPntValid; }
  int GetRGId(int i) const { return fkReGroup.size() ? (fkReGroup[i] < 0 ? -1 : fkReGroup[i]) : i; }
  int GetProcessedPoints(int i) const
  {
    int ir = GetRGId(i);
    return ir <= 0 ? 0 : fProcPnt[ir];
  }
  std::vector<int> GetProcessedPoints() const { return fProcPnt; }
  int GetParamGrID(int i) const
  {
    int ir = GetRGId(i);
    return ir <= 0 ? 0 : fParamGrID[ir];
  }
 
  MatrixSq* GetGlobalMatrix() const { return fMatCGlo; }
  SymMatrix* GetLocalMatrix() const { return fMatCLoc; }
  std::vector<double> GetGlobals() const { return fVecBGlo; }
  std::vector<double> GetDeltaPars() const { return fDeltaPar; }
  std::vector<double> GetInitPars() const { return fInitPar; }
  std::vector<double> GetSigmaPars() const { return fSigmaPar; }
  std::vector<bool> GetIsLinear() const { return fIsLinear; }
  double GetFinalParam(int i) const
  {
    int ir = GetRGId(i);
    return ir < 0 ? 0 : fDeltaPar[ir] + fInitPar[ir];
  }
  double GetFinalError(int i) const { return GetParError(i); }
 
  double GetPull(int i) const;
 
  double GetGlobal(int i) const
  {
    int ir = GetRGId(i);
    return ir < 0 ? 0 : fVecBGlo[ir];
  }
  double GetInitPar(int i) const
  {
    int ir = GetRGId(i);
    return ir < 0 ? 0 : fInitPar[ir];
  }
  double GetSigmaPar(int i) const
  {
    int ir = GetRGId(i);
    return ir < 0 ? 0 : fSigmaPar[ir];
  }
  bool GetIsLinear(int i) const
  {
    int ir = GetRGId(i);
    return ir < 0 ? 0 : fIsLinear[ir];
  }
  static bool IsGlobalMatSparse() { return fgIsMatGloSparse; }
  static bool IsWeightSigma() { return fgWeightSigma; }
  void SetWghScale(const double wOdd = 1, const double wEven = 1)
  {
    fWghScl[0] = wOdd;
    fWghScl[1] = wEven;
  }
  void SetUseRecordWeight(const bool v = true) { fUseRecordWeight = v; }
  bool GetUseRecordWeight() const { return fUseRecordWeight; }
  void SetMinRecordLength(const int v = 1) { fMinRecordLength = v; }
  int GetMinRecordLength() const { return fMinRecordLength; }
 
  void SetParamGrID(const int grID, int i)
  {
    int ir = GetRGId(i);
    if (ir < 0) {
      return;
    }
    fParamGrID[ir] = grID;
    if (fNGroupsSet < grID) {
      fNGroupsSet = grID;
    }
  }
  void ResetRecord() { fRecord->Reset(); }
  void SetNGloPar(const int n) { fNGloPar = n; }
  void SetNLocPar(const int n) { fNLocPar = n; }
  void SetNMaxIterations(const int n = 10) { fMaxIter = n; }
  void SetNStdDev(const int n) { fNStdDev = n; }
  void SetChi2CutFactor(const float v) { fChi2CutFactor = v; }
  void SetChi2CutRef(const float v) { fChi2CutRef = v; }
  void SetResCurInit(const float v) { fResCutInit = v; }
  void SetResCut(const float v) { fResCut = v; }
  void SetMinPntValid(const int n) { fMinPntValid = n > 0 ? n : 1; }
  static void SetGlobalMatSparse(const bool v = true) { fgIsMatGloSparse = v; }
  static void SetWeightSigma(const bool v = true) { fgWeightSigma = v; }
 
  void SetInitPars(const double* par);
 
  void SetSigmaPars(const double* par);
 
  void SetInitPar(int i, double par);
 
  void SetSigmaPar(int i, double par);
 
  int GlobalFit(std::vector<double>& par,
                std::vector<double>& error,
                std::vector<double>& pull);
 
  int GlobalFitIteration();
 
  int SolveGlobalMatEq();
 
  static void SetInvChol(const bool v = true) { fgInvChol = v; }
  static void SetMinResPrecondType(const int tp = 0) { fgMinResCondType = tp; }
  static void SetMinResTol(double val = 1e-12) { fgMinResTol = val; }
  static void SetMinResMaxIter(const int val = 2000) { fgMinResMaxIter = val; }
  static void SetIterSolverType(const int val = MinResSolve::kSolMinRes) { fgIterSol = val; }
  static void SetNKrylovV(const int val = 60) { fgNKrylovV = val; }
 
  static bool GetInvChol() { return fgInvChol; }
  static int GetMinResPrecondType() { return fgMinResCondType; }
  static double GetMinResTol() { return fgMinResTol; }
  static int GetMinResMaxIter() { return fgMinResMaxIter; }
  static int GetIterSolverType() { return fgIterSol; }
  static int GetNKrylovV() { return fgNKrylovV; }
 
  double GetParError(int iPar) const;
 
  int PrintGlobalParameters() const;
 
  void SetRejRunList(const int* runs, const int nruns);
 
  void SetAccRunList(const int* runs, const int nruns, const float* wghList = nullptr);
 
  bool IsRecordAcceptable();
 
  int SetIterations(const double lChi2CutFac);
 
  // constraints
 
  void SetGlobalConstraint(const std::vector<double>& dergb,
                           const double val, const double sigma = 0,
                           const bool doPrint = false);
 
  void SetGlobalConstraint(const std::vector<int>& indgb,
                           const std::vector<double>& dergb,
                           const int ngb, const double val,
                           double sigma = 0, const bool doPrint = false);
 
  void SetLocalEquation(std::vector<double>& dergb, std::vector<double>& derlc,
                        const double lMeas, const double lSigma);
 
  void SetLocalEquation(std::vector<int>& indgb, std::vector<double>& dergb,
                        int ngb, std::vector<int>& indlc,
                        std::vector<double>& derlc, const int nlc,
                        const double lMeas, const double lSigma);
 
  const char* GetRecChi2FName() const { return fRecChi2FName.Data(); }
 
  bool InitChi2Storage(const int nEntriesAutoSave = 10000);
 
  void EndChi2Storage();
 
  o2::fwdalign::MillePedeRecord* GetRecord() const { return fRecord; }
  long GetSelFirst() const { return fSelFirst; }
  long GetSelLast() const { return fSelLast; }
  void SetSelFirst(long v) { fSelFirst = v; }
  void SetSelLast(long v) { fSelLast = v; }
 
  void SetRecord(o2::fwdalign::MillePedeRecord* aRecord) { fRecord = aRecord; }
  void SetRecordWriter(o2::fwdalign::MilleRecordWriter* myP) { fRecordWriter = myP; }
  void SetConstraintsRecWriter(o2::fwdalign::MilleRecordWriter* myP) { fConstraintsRecWriter = myP; }
  void SetRecordReader(o2::fwdalign::MilleRecordReader* myP) { fRecordReader = myP; }
  void SetConstraintsRecReader(o2::fwdalign::MilleRecordReader* myP) { fConstraintsRecReader = myP; }
 
  float Chi2DoFLim(int nSig, int nDoF) const;
 
  // aliases for compatibility with millipede1
  void SetParSigma(int i, double par) { SetSigmaPar(i, par); }
  void SetGlobalParameters(double* par) { SetInitPars(par); }
  void SetNonLinear(int index, bool v = true)
  {
    int id = GetRGId(index);
    if (id < 0) {
      return;
    }
    fIsLinear[id] = !v;
  }
 
  void DisableRecordWriter() { fDisableRecordWriter = true; }
 
 protected:
  void ReadRecordData(const long recID, const bool doPrint = false);
 
  void ReadRecordConstraint(const long recID, const bool doPrint = false);
 
  int LocalFit(std::vector<double>& localParams);
 
  bool IsZero(const double v, const double eps = 1e-16) const { return TMath::Abs(v) < eps; }
 
 protected:
  int fNLocPar;    
  int fNGloPar;    
  int fNGloParIni; 
  int fNGloSize;   
 
  long fNLocEquations;       
  int fIter;                 
  int fMaxIter;              
  int fNStdDev;              
  int fNGloConstraints;      
  int fNLagrangeConstraints; 
  long fNLocFits;            
  long fNLocFitsRejected;    
  int fNGloFix;              
  int fGloSolveStatus;       
 
  float fChi2CutFactor; 
  float fChi2CutRef;    
  float fResCutInit;    
  float fResCut;        
  int fMinPntValid;     
 
  int fNGroupsSet;               
  std::vector<int> fParamGrID;   
  std::vector<int> fProcPnt;     
  std::vector<double> fVecBLoc;  
  std::vector<double> fDiagCGlo; 
  std::vector<double> fVecBGlo;  
 
  std::vector<double> fInitPar;  
  std::vector<double> fDeltaPar; 
  std::vector<double> fSigmaPar; 
 
  std::vector<bool> fIsLinear;   
  std::vector<bool> fConstrUsed; 
 
  std::vector<int> fGlo2CGlo; 
  std::vector<int> fCGlo2Glo; 
 
  // Matrices
  o2::fwdalign::SymMatrix* fMatCLoc;     
  o2::fwdalign::MatrixSq* fMatCGlo;      
  o2::fwdalign::RectMatrix* fMatCGloLoc; 
  std::vector<int> fFillIndex;           
  std::vector<double> fFillValue;        
 
  TFile* fRecChi2File;
  TString fRecChi2FName;
  TString fRecChi2TreeName; 
  TTree* fTreeChi2;
  float fSumChi2;
  bool fIsChi2BelowLimit;
  int fRecNDoF;
 
  o2::fwdalign::MillePedeRecord* fRecord; 
 
  long fCurrRecDataID;        
  long fCurrRecConstrID;      
  bool fLocFitAdd;            
  bool fUseRecordWeight;      
  bool fDisableRecordWriter;  
  int fMinRecordLength;       
  int fSelFirst;              
  int fSelLast;               
  TArrayL* fRejRunList;       
  TArrayL* fAccRunList;       
  TArrayF* fAccRunListWgh;    
  double fRunWgh;             
  double fWghScl[2];          
  std::vector<int> fkReGroup; 
 
  static bool fgInvChol;        
  static bool fgWeightSigma;    
  static bool fgIsMatGloSparse; 
  static int fgMinResCondType;  
  static double fgMinResTol;    
  static int fgMinResMaxIter;   
  static int fgIterSol;         
  static int fgNKrylovV;        
 
  // processed data record bufferization
  o2::fwdalign::MilleRecordWriter* fRecordWriter;         
  o2::fwdalign::MilleRecordWriter* fConstraintsRecWriter; 
  o2::fwdalign::MilleRecordReader* fRecordReader;         
  o2::fwdalign::MilleRecordReader* fConstraintsRecReader; 
 
  ClassDef(MillePede2, 0);
};
 
} // namespace fwdalign
} // namespace o2
 
#endif