#ifndef ALIMILLEPEDE2_H
#define ALIMILLEPEDE2_H
#include <TObject.h>
#include <TString.h>
#include <TTree.h>
#include "AliMinResSolve.h"
#include "AliMillePedeRecord.h"
class TFile;
class AliMatrixSq;
class AliSymMatrix;
class AliRectMatrix;
class AliMatrixSparse;
class AliLog;
class TStopwatch;
class TArrayL;
class TArrayF;
class AliMillePede2: public TObject
{
public:
enum {kFailed,kInvert,kNoInversion};
enum {kFixParID=-1};
AliMillePede2();
AliMillePede2(const AliMillePede2& src);
virtual ~AliMillePede2();
AliMillePede2& operator=(const AliMillePede2& ) {printf("Dummy\n"); return *this;}
Int_t InitMille(int nGlo, int nLoc, Int_t lNStdDev=-1,double lResCut=-1., double lResCutInit=-1., const Int_t* regroup=0);
Int_t GetNGloPar() const {return fNGloPar;}
Int_t GetNGloParIni() const {return fNGloParIni;}
const Int_t* GetRegrouping() const {return fkReGroup;}
Int_t GetNLocPar() const {return fNLocPar;}
Long_t GetNLocalEquations() const {return fNLocEquations;}
Int_t GetCurrentIteration() const {return fIter;}
Int_t GetNMaxIterations() const {return fMaxIter;}
Int_t GetNStdDev() const {return fNStdDev;}
Int_t GetNGlobalConstraints() const {return fNGloConstraints;}
Int_t GetNLagrangeConstraints() const {return fNLagrangeConstraints;}
Long_t GetNLocalFits() const {return fNLocFits;}
Long_t GetNLocalFitsRejected() const {return fNLocFitsRejected;}
Int_t GetNGlobalsFixed() const {return fNGloFix;}
Int_t GetGlobalSolveStatus() const {return fGloSolveStatus;}
Float_t GetChi2CutFactor() const {return fChi2CutFactor;}
Float_t GetChi2CutRef() const {return fChi2CutRef;}
Float_t GetResCurInit() const {return fResCutInit;}
Float_t GetResCut() const {return fResCut;}
Int_t GetMinPntValid() const {return fMinPntValid;}
Int_t GetRGId(Int_t i) const {return fkReGroup ? (fkReGroup[i]<0 ? -1:fkReGroup[i]) : i;}
Int_t GetProcessedPoints(Int_t i) const {int ir=GetRGId(i); return ir<=0 ? 0:fProcPnt[ir];}
Int_t* GetProcessedPoints() const {return fProcPnt;}
Int_t GetParamGrID(Int_t i) const {int ir=GetRGId(i); return ir<=0 ? 0:fParamGrID[ir];}
AliMatrixSq* GetGlobalMatrix() const {return fMatCGlo;}
AliSymMatrix* GetLocalMatrix() const {return fMatCLoc;}
Double_t* GetGlobals() const {return fVecBGlo;}
Double_t* GetDeltaPars() const {return fDeltaPar;}
Double_t* GetInitPars() const {return fInitPar;}
Double_t* GetSigmaPars() const {return fSigmaPar;}
Bool_t* GetIsLinear() const {return fIsLinear;}
Double_t GetFinalParam(int i) const {int ir=GetRGId(i); return ir<0 ? 0:fDeltaPar[ir]+fInitPar[ir];}
Double_t GetFinalError(int i) const {return GetParError(i);}
Double_t GetPull(int i) const;
Double_t GetGlobal(Int_t i) const {int ir=GetRGId(i); return ir<0 ? 0:fVecBGlo[ir];}
Double_t GetInitPar(Int_t i) const {int ir=GetRGId(i); return ir<0 ? 0:fInitPar[ir];}
Double_t GetSigmaPar(Int_t i) const {int ir=GetRGId(i); return ir<0 ? 0:fSigmaPar[ir];}
Bool_t GetIsLinear(Int_t i) const {int ir=GetRGId(i); return ir<0 ? 0:fIsLinear[ir];}
static Bool_t IsGlobalMatSparse() {return fgIsMatGloSparse;}
static Bool_t IsWeightSigma() {return fgWeightSigma;}
void SetWghScale(Double_t wOdd=1,Double_t wEven=1) {fWghScl[0] = wOdd; fWghScl[1] = wEven;}
void SetUseRecordWeight(Bool_t v=kTRUE) {fUseRecordWeight=v;}
Bool_t GetUseRecordWeight() const {return fUseRecordWeight;}
void SetMinRecordLength(Int_t v=1) {fMinRecordLength = v;}
Int_t GetMinRecordLength() const {return fMinRecordLength;}
void SetParamGrID(Int_t grID,Int_t i) {int ir=GetRGId(i); if(ir<0) return; fParamGrID[ir] = grID; if(fNGroupsSet<grID)fNGroupsSet=grID;}
void SetNGloPar(Int_t n) {fNGloPar = n;}
void SetNLocPar(Int_t n) {fNLocPar = n;}
void SetNMaxIterations(Int_t n=10) {fMaxIter = n;}
void SetNStdDev(Int_t n) {fNStdDev = n;}
void SetChi2CutFactor(Float_t v) {fChi2CutFactor = v;}
void SetChi2CutRef(Float_t v) {fChi2CutRef = v;}
void SetResCurInit(Float_t v) {fResCutInit = v;}
void SetResCut(Float_t v) {fResCut = v;}
void SetMinPntValid(Int_t n) {fMinPntValid = n>0 ? n:1;}
static void SetGlobalMatSparse(Bool_t v=kTRUE) {fgIsMatGloSparse = v;}
static void SetWeightSigma(Bool_t v=kTRUE) {fgWeightSigma = v;}
void SetInitPars(const Double_t* par);
void SetSigmaPars(const Double_t* par);
void SetInitPar(Int_t i,Double_t par);
void SetSigmaPar(Int_t i,Double_t par);
Int_t GlobalFit(Double_t *par=0, Double_t *error=0, Double_t *pull=0);
Int_t GlobalFitIteration();
Int_t SolveGlobalMatEq();
static void SetInvChol(Bool_t v=kTRUE) {fgInvChol = v;}
static void SetMinResPrecondType(Int_t tp=0) {fgMinResCondType = tp;}
static void SetMinResTol(Double_t val=1e-12) {fgMinResTol = val;}
static void SetMinResMaxIter(Int_t val=2000) {fgMinResMaxIter = val;}
static void SetIterSolverType(Int_t val=AliMinResSolve::kSolMinRes) {fgIterSol = val;}
static void SetNKrylovV(Int_t val=60) {fgNKrylovV = val;}
static Bool_t GetInvChol() {return fgInvChol;}
static Int_t GetMinResPrecondType() {return fgMinResCondType;}
static Double_t GetMinResTol() {return fgMinResTol;}
static Int_t GetMinResMaxIter() {return fgMinResMaxIter;}
static Int_t GetIterSolverType() {return fgIterSol;}
static Int_t GetNKrylovV() {return fgNKrylovV;}
Double_t GetParError(int iPar) const;
Int_t PrintGlobalParameters() const;
void SetRejRunList(const UInt_t *runs, Int_t nruns);
void SetAccRunList(const UInt_t *runs, Int_t nruns, const Float_t* wghList=0);
Bool_t IsRecordAcceptable();
Int_t SetIterations(double lChi2CutFac);
void SetGlobalConstraint(const double *dergb, double val, double sigma=0);
void SetGlobalConstraint(const int *indgb, const double *dergb, int ngb, double val, double sigma=0);
void SetRecordRun(Int_t run);
void SetRecordWeight(double wgh);
void SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma);
void SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc,
double *derlc,int nlc,double lMeas,double lSigma);
void SetDataRecFName(const char* flname) {fDataRecFName = flname;}
const Char_t* GetDataRecFName() const {return fDataRecFName.Data();}
void SetConsRecFName(const char* flname) {fConstrRecFName = flname;}
const Char_t* GetConsRecFName() const {return fConstrRecFName.Data();}
void SetRecDataTreeName(const char* name=0) {fRecDataTreeName = name; if (fRecDataTreeName.IsNull()) fRecDataTreeName = "AliMillePedeRecords_Data";}
void SetRecConsTreeName(const char* name=0) {fRecConsTreeName = name; if (fRecConsTreeName.IsNull()) fRecConsTreeName = "AliMillePedeRecords_Consaints";}
void SetRecDataBranchName(const char* name=0) {fRecDataBranchName = name; if (fRecDataBranchName.IsNull()) fRecDataBranchName = "Record_Data";}
void SetRecConsBranchName(const char* name=0) {fRecConsBranchName = name; if (fRecConsBranchName.IsNull()) fRecConsBranchName = "Record_Consaints";}
const char* GetRecDataTreeName() const {return fRecDataTreeName.Data();}
const char* GetRecConsTreeName() const {return fRecConsTreeName.Data();}
const char* GetRecDataBranchName() const {return fRecDataBranchName.Data();}
const char* GetRecConsBranchName() const {return fRecConsBranchName.Data();}
Bool_t InitDataRecStorage(Bool_t read=kFALSE);
Bool_t InitConsRecStorage(Bool_t read=kFALSE);
Bool_t ImposeDataRecFile(const char* fname);
Bool_t ImposeConsRecFile(const char* fname);
void CloseDataRecStorage();
void CloseConsRecStorage();
void ReadRecordData(Long_t recID);
void ReadRecordConstraint(Long_t recID);
Bool_t ReadNextRecordData();
Bool_t ReadNextRecordConstraint();
void SaveRecordData();
void SaveRecordConstraint();
AliMillePedeRecord* GetRecord() const {return fRecord;}
Long_t GetSelFirst() const {return fSelFirst;}
Long_t GetSelLast() const {return fSelLast;}
void SetSelFirst(Long_t v) {fSelFirst = v;}
void SetSelLast(Long_t v) {fSelLast = v;}
Float_t Chi2DoFLim(int nSig, int nDoF) const;
void SetParSigma(Int_t i,Double_t par) {SetSigmaPar(i,par);}
void SetGlobalParameters(Double_t *par) {SetInitPars(par);}
void SetNonLinear(int index, Bool_t v=kTRUE) {int id = GetRGId(index); if (id<0) return; fIsLinear[id] = !v;}
protected:
Int_t LocalFit(double *localParams=0);
Bool_t IsZero(Double_t v,Double_t eps=1e-16) const {return TMath::Abs(v)<eps;}
protected:
Int_t fNLocPar;
Int_t fNGloPar;
Int_t fNGloParIni;
Int_t fNGloSize;
Long_t fNLocEquations;
Int_t fIter;
Int_t fMaxIter;
Int_t fNStdDev;
Int_t fNGloConstraints;
Int_t fNLagrangeConstraints;
Long_t fNLocFits;
Long_t fNLocFitsRejected;
Int_t fNGloFix;
Int_t fGloSolveStatus;
Float_t fChi2CutFactor;
Float_t fChi2CutRef;
Float_t fResCutInit;
Float_t fResCut;
Int_t fMinPntValid;
Int_t fNGroupsSet;
Int_t *fParamGrID;
Int_t *fProcPnt;
Double_t *fVecBLoc;
Double_t *fDiagCGlo;
Double_t *fVecBGlo;
Double_t *fInitPar;
Double_t *fDeltaPar;
Double_t *fSigmaPar;
Bool_t *fIsLinear;
Bool_t *fConstrUsed;
Int_t *fGlo2CGlo;
Int_t *fCGlo2Glo;
AliSymMatrix *fMatCLoc;
AliMatrixSq *fMatCGlo;
AliRectMatrix *fMatCGloLoc;
Int_t *fFillIndex;
Double_t *fFillValue;
TString fRecDataTreeName;
TString fRecConsTreeName;
TString fRecDataBranchName;
TString fRecConsBranchName;
TString fDataRecFName;
AliMillePedeRecord *fRecord;
TFile *fDataRecFile;
TTree *fTreeData;
Int_t fRecFileStatus;
TString fConstrRecFName;
TTree *fTreeConstr;
TFile *fConsRecFile;
Long_t fCurrRecDataID;
Long_t fCurrRecConstrID;
Bool_t fLocFitAdd;
Bool_t fUseRecordWeight;
Int_t fMinRecordLength;
Int_t fSelFirst;
Int_t fSelLast;
TArrayL* fRejRunList;
TArrayL* fAccRunList;
TArrayF* fAccRunListWgh;
Double_t fRunWgh;
Double_t fWghScl[2];
const Int_t* fkReGroup;
static Bool_t fgInvChol;
static Bool_t fgWeightSigma;
static Bool_t fgIsMatGloSparse;
static Int_t fgMinResCondType;
static Double_t fgMinResTol;
static Int_t fgMinResMaxIter;
static Int_t fgIterSol;
static Int_t fgNKrylovV;
ClassDef(AliMillePede2,1)
};
inline void AliMillePede2::ReadRecordData(Long_t recID) {fTreeData->GetEntry(recID); fCurrRecDataID=recID;}
inline void AliMillePede2::ReadRecordConstraint(Long_t recID) {fTreeConstr->GetEntry(recID); fCurrRecConstrID=recID;}
inline void AliMillePede2::SaveRecordData() {fTreeData->Fill(); fRecord->Reset(); fCurrRecDataID++;}
inline void AliMillePede2::SaveRecordConstraint() {fTreeConstr->Fill(); fRecord->Reset();fCurrRecConstrID++;}
#endif