#include <TMath.h>
#include <TVectorD.h>
#include <TArrayL.h>
+#include <TArrayF.h>
#include <TSystem.h>
#include "AliMatrixSq.h"
#include "AliSymMatrix.h"
#include <fcntl.h>
#include <fstream>
-using namespace std;
+//#define _DUMP_EQ_BEFORE_
+//#define _DUMP_EQ_AFTER_
+//#define _DUMPEQ_BEFORE_
+//#define _DUMPEQ_AFTER_
+using std::ifstream;
ClassImp(AliMillePede2)
Bool_t AliMillePede2::fgInvChol = kTRUE; // Invert global matrix with Cholesky solver
AliMillePede2::AliMillePede2()
: fNLocPar(0),
fNGloPar(0),
+ fNGloParIni(0),
fNGloSize(0),
//
fNLocEquations(0),
fFillIndex(0),
fFillValue(0),
//
+ fRecDataTreeName("AliMillePedeRecords_Data"),
+ fRecConsTreeName("AliMillePedeRecords_Consaints"),
+ fRecDataBranchName("Record_Data"),
+ fRecConsBranchName("Record_Consaints"),
+ //
fDataRecFName("/tmp/mp2_data_records.root"),
fRecord(0),
fDataRecFile(0),
fCurrRecDataID(0),
fCurrRecConstrID(0),
fLocFitAdd(kTRUE),
+ fUseRecordWeight(kTRUE),
+ fMinRecordLength(1),
fSelFirst(1),
fSelLast(-1),
fRejRunList(0),
- fAccRunList(0)
-{}
+ fAccRunList(0),
+ fAccRunListWgh(0),
+ fRunWgh(1),
+ fkReGroup(0)
+{
+ fWghScl[0] = fWghScl[1] = -1;
+}
//_____________________________________________________________________________________________
AliMillePede2::AliMillePede2(const AliMillePede2& src) :
- TObject(src),fNLocPar(0),fNGloPar(0),fNGloSize(0),fNLocEquations(0),fIter(0),
+ TObject(src),fNLocPar(0),fNGloPar(0),fNGloParIni(0),fNGloSize(0),fNLocEquations(0),fIter(0),
fMaxIter(10),fNStdDev(3),fNGloConstraints(0),fNLagrangeConstraints(0),
fNLocFits(0),fNLocFitsRejected(0),
fNGloFix(0),fGloSolveStatus(0),fChi2CutFactor(0),fChi2CutRef(0),fResCutInit(0),
fResCut(0),fMinPntValid(1),fNGroupsSet(0),fParamGrID(0),fProcPnt(0),fVecBLoc(0),fDiagCGlo(0),fVecBGlo(0),
fInitPar(0),fDeltaPar(0),fSigmaPar(0),fIsLinear(0),fConstrUsed(0),fGlo2CGlo(0),fCGlo2Glo(0),
fMatCLoc(0),fMatCGlo(0),fMatCGloLoc(0),fFillIndex(0),fFillValue(0),
+ fRecDataTreeName(0),fRecConsTreeName(0),fRecDataBranchName(0),fRecConsBranchName(0),
fDataRecFName(0),fRecord(0),fDataRecFile(0),
fTreeData(0),fRecFileStatus(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),
fCurrRecConstrID(0),fLocFitAdd(kTRUE),
+ fUseRecordWeight(kTRUE),
+ fMinRecordLength(1),
fSelFirst(1),
fSelLast(-1),
fRejRunList(0),
- fAccRunList(0)
-{printf("Dummy\n");}
+ fAccRunList(0),
+ fAccRunListWgh(0),
+ fRunWgh(1),
+ fkReGroup(0)
+{
+ fWghScl[0] = src.fWghScl[0];
+ fWghScl[1] = src.fWghScl[1];
+ printf("Dummy\n");
+}
//_____________________________________________________________________________________________
AliMillePede2::~AliMillePede2()
delete fMatCGloLoc;
delete fRejRunList;
delete fAccRunList;
+ delete fAccRunListWgh;
}
//_____________________________________________________________________________________________
-Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut, double lResCutInit)
+Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut, double lResCutInit, const Int_t* regroup)
{
// init all
+ //
+ fNGloParIni = nGlo;
+ if (regroup) { // regrouping is requested
+ fkReGroup = regroup;
+ int ng = 0; // recalculate N globals
+ int maxPID = -1;
+ for (int i=0;i<nGlo;i++) if (regroup[i]>=0) {ng++; if (regroup[i]>maxPID) maxPID = regroup[i];}
+ maxPID++;
+ AliInfo(Form("Regrouping is requested: from %d raw to %d formal globals grouped to %d real globals",nGlo,ng,maxPID));
+ nGlo = maxPID;
+ }
if (nLoc>0) fNLocPar = nLoc;
if (nGlo>0) fNGloPar = nGlo;
if (lResCutInit>0) fResCutInit = lResCutInit;
if (lResCut>0) fResCut = lResCut;
if (lNStdDev>0) fNStdDev = lNStdDev;
//
+ AliInfo(Form("NLoc: %d NGlo: %d",fNLocPar,fNGloPar));
+
fNGloSize = fNGloPar;
//
- try {
- //
- if (fgIsMatGloSparse) {fMatCGlo = new AliMatrixSparse(fNGloPar); fMatCGlo->SetSymmetric(kTRUE);}
- else fMatCGlo = new AliSymMatrix(fNGloPar);
- //
- fFillIndex = new Int_t[fNGloPar];
- fFillValue = new Double_t[fNGloPar];
- //
- fMatCLoc = new AliSymMatrix(fNLocPar);
- fMatCGloLoc = new AliRectMatrix(fNGloPar,fNLocPar);
- //
- fParamGrID = new Int_t[fNGloPar];
- fProcPnt = new Int_t[fNGloPar];
- fVecBLoc = new Double_t[fNLocPar];
- fDiagCGlo = new Double_t[fNGloPar];
- //
- fInitPar = new Double_t[fNGloPar];
- fDeltaPar = new Double_t[fNGloPar];
- fSigmaPar = new Double_t[fNGloPar];
- fIsLinear = new Bool_t[fNGloPar];
- //
- fGlo2CGlo = new Int_t[fNGloPar];
- fCGlo2Glo = new Int_t[fNGloPar];
- }
- catch(bad_alloc&) {
- AliInfo(Form("Failed to allocate the memory for %d global and %d local parameters",fNGloPar,fNLocPar));
- return 0;
- }
+ if (fgIsMatGloSparse) {fMatCGlo = new AliMatrixSparse(fNGloPar); fMatCGlo->SetSymmetric(kTRUE);}
+ else fMatCGlo = new AliSymMatrix(fNGloPar);
+ //
+ fFillIndex = new Int_t[fNGloPar];
+ fFillValue = new Double_t[fNGloPar];
+ //
+ fMatCLoc = new AliSymMatrix(fNLocPar);
+ fMatCGloLoc = new AliRectMatrix(fNGloPar,fNLocPar);
+ //
+ fParamGrID = new Int_t[fNGloPar];
+ fProcPnt = new Int_t[fNGloPar];
+ fVecBLoc = new Double_t[fNLocPar];
+ fDiagCGlo = new Double_t[fNGloPar];
+ //
+ fInitPar = new Double_t[fNGloPar];
+ fDeltaPar = new Double_t[fNGloPar];
+ fSigmaPar = new Double_t[fNGloPar];
+ fIsLinear = new Bool_t[fNGloPar];
+ //
+ fGlo2CGlo = new Int_t[fNGloPar];
+ fCGlo2Glo = new Int_t[fNGloPar];
//
memset(fVecBLoc ,0,fNLocPar*sizeof(Double_t));
memset(fDiagCGlo ,0,fNGloPar*sizeof(Double_t));
fParamGrID[i] = -1;
}
//
+ fWghScl[0] = -1;
+ fWghScl[1] = -1;
return 1;
}
fDataRecFile = TFile::Open(GetDataRecFName(),"recreate");
if (!fDataRecFile) {AliFatal(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;}
AliInfo(Form("File %s used for derivatives records",GetDataRecFName()));
- fTreeData = new TTree("AliMillePedeRecords_Data","Data Records for AliMillePede2");
- fTreeData->Branch("Record_Data","AliMillePedeRecord",&fRecord,32000,99);
+ fTreeData = new TTree(GetRecDataTreeName(),"Data Records for AliMillePede2");
+ fTreeData->Branch(GetRecDataBranchName(),"AliMillePedeRecord",&fRecord,32000,99);
}
else { // use chain
- TChain* ch = new TChain("AliMillePedeRecords_Data");
+ TChain* ch = new TChain(GetRecDataTreeName());
//
if (fDataRecFName.EndsWith(".root")) ch->AddFile(fDataRecFName);
else { // assume text file with list of filenames
Long64_t nent = ch->GetEntries();
if (nent<1) { AliInfo("Obtained chain is empty"); return kFALSE;}
fTreeData = ch;
- fTreeData->SetBranchAddress("Record_Data",&fRecord);
+ fTreeData->SetBranchAddress(GetRecDataBranchName(),&fRecord);
AliInfo(Form("Found %lld derivatives records",nent));
}
fCurrRecDataID = -1;
//
AliInfo(Form("File %s used for constraints records",GetConsRecFName()));
if (read) {
- fTreeConstr = (TTree*)fConsRecFile->Get("AliMillePedeRecords_Constraints");
+ fTreeConstr = (TTree*)fConsRecFile->Get(GetRecConsTreeName());
if (!fTreeConstr) {AliInfo(Form("Did not find constraints records tree in %s",GetConsRecFName())); return kFALSE;}
- fTreeConstr->SetBranchAddress("Record_Constraints",&fRecord);
+ fTreeConstr->SetBranchAddress(GetRecConsBranchName(),&fRecord);
AliInfo(Form("Found %lld constraints records",fTreeConstr->GetEntries()));
//
}
else {
//
- fTreeConstr = new TTree("AliMillePedeRecords_Constraints","Constraints Records for AliMillePede2");
- fTreeConstr->Branch("Record_Constraints","AliMillePedeRecord",&fRecord,32000,99);
+ fTreeConstr = new TTree(GetRecConsTreeName(),"Constraints Records for AliMillePede2");
+ fTreeConstr->Branch(GetRecConsBranchName(),"AliMillePedeRecord",&fRecord,32000,99);
}
fCurrRecConstrID = -1;
//
// write data of single measurement
if (lSigma<=0.0) { // If parameter is fixed, then no equation
for (int i=fNLocPar; i--;) derlc[i] = 0.0;
- for (int i=fNGloPar; i--;) dergb[i] = 0.0;
+ for (int i=fNGloParIni; i--;) dergb[i] = 0.0;
return;
}
//
fRecord->AddWeight( 1.0/lSigma/lSigma );
//
// Idem for global parameters
- for (int i=0;i<fNGloPar;i++) if (!IsZero(dergb[i])) {
+ for (int i=0;i<fNGloParIni;i++) if (!IsZero(dergb[i])) {
fRecord->AddIndexValue(i,dergb[i]); dergb[i] = 0.0;
- fRecord->MarkGroup(fParamGrID[i]);
+ int idrg = GetRGId(i);
+ fRecord->MarkGroup(idrg<0 ? -1 : fParamGrID[i]);
}
+ // fRecord->Print();
//
}
void AliMillePede2::SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc,
double *derlc,int nlc,double lMeas,double lSigma)
{
- // write data of single measurement
+ // write data of single measurement. Note: the records ignore regrouping, store direct parameters
if (lSigma<=0.0) { // If parameter is fixed, then no equation
for (int i=nlc;i--;) derlc[i] = 0.0;
for (int i=ngb;i--;) dergb[i] = 0.0;
fRecord->Reset();
fRecord->AddResidual(val);
fRecord->AddWeight(sigma);
- for (int i=0; i<fNGloPar; i++) if (!IsZero(dergb[i])) fRecord->AddIndexValue(i,dergb[i]);
+ for (int i=0; i<fNGloParIni; i++) if (!IsZero(dergb[i])) fRecord->AddIndexValue(i,dergb[i]);
fNGloConstraints++;
if (IsZero(sigma)) fNLagrangeConstraints++;
// printf("NewConstraint:\n"); fRecord->Print(); //RRR
//
nPoints++;
}
+ if (fMinRecordLength>0 && nPoints < fMinRecordLength) return 0; // ignore
+ //
double vl;
//
- double gloWgh = fRecord->GetWeight(); // global weight for this set
+ double gloWgh = fRunWgh;
+ if (fUseRecordWeight) gloWgh *= fRecord->GetWeight(); // global weight for this set
Int_t maxLocUsed = 0;
//
for (int ip=nPoints;ip--;) { // Transfer the measurement records to matrices
double resid = fRecord->GetValue( refLoc[ip]-1 );
- double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
+ double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
+ int odd = (ip&0x1);
+ if (fWghScl[odd]>0) weight *= fWghScl[odd];
double *derLoc = fRecord->GetValue()+refLoc[ip];
double *derGlo = fRecord->GetValue()+refGlo[ip];
int *indLoc = fRecord->GetIndex()+refLoc[ip];
int *indGlo = fRecord->GetIndex()+refGlo[ip];
//
for (int i=nrefGlo[ip];i--;) { // suppress the global part (only relevant with iterations)
+ //
+ // if regrouping was requested, do it here
+ if (fkReGroup) {
+ int idtmp = fkReGroup[ indGlo[i] ];
+ if (idtmp == kFixParID) indGlo[i] = kFixParID; // fixed param in regrouping
+ else indGlo[i] = idtmp;
+ }
+ //
int iID = indGlo[i]; // Global param indice
- if (fSigmaPar[iID]<=0.) continue; // fixed parameter RRRCheck
+ if (iID<0 || fSigmaPar[iID]<=0.) continue; // fixed parameter RRRCheck
if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter
else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter
}
//
for (int ip=nPoints;ip--;) { // Calculate residuals
double resid = fRecord->GetValue( refLoc[ip]-1 );
- double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
+ double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
+ int odd = (ip&0x1);
+ if (fWghScl[odd]>0) weight *= fWghScl[odd];
double *derLoc = fRecord->GetValue()+refLoc[ip];
double *derGlo = fRecord->GetValue()+refGlo[ip];
int *indLoc = fRecord->GetIndex()+refLoc[ip];
//
for (int i=nrefGlo[ip];i--;) { // global part
int iID = indGlo[i];
- if ( fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck
+ if ( iID<0 || fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck
if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter
else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter
}
//
for (int ip=nPoints;ip--;) { // Update matrices
double resid = fRecord->GetValue( refLoc[ip]-1 );
- double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
+ double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
+ int odd = (ip&0x1);
+ if (fWghScl[odd]>0) weight *= fWghScl[odd];
double *derLoc = fRecord->GetValue()+refLoc[ip];
double *derGlo = fRecord->GetValue()+refGlo[ip];
int *indLoc = fRecord->GetIndex()+refLoc[ip];
int *indGlo = fRecord->GetIndex()+refGlo[ip];
//
- for (int i=nrefGlo[ip];i--;) { // suppress the global part
+ for (int i=nrefGlo[ip];i--;) { // suppress the global part
int iID = indGlo[i]; // Global param indice
- if ( fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck
+ if ( iID<0 || fSigmaPar[iID] <= 0.) continue; // fixed parameter RRRCheck
if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter
else resid -= derGlo[i]*fDeltaPar[iID]; // nonlinear parameter
}
//
for (int ig=nrefGlo[ip];ig--;) {
int iIDg = indGlo[ig]; // Global param indice (the matrix line)
- if ( fSigmaPar[iIDg] <= 0.) continue; // fixed parameter RRRCheck
+ if ( iIDg<0 || fSigmaPar[iIDg] <= 0.) continue; // fixed parameter RRRCheck
if (fLocFitAdd) fVecBGlo[ iIDg ] += weight*resid*derGlo[ig]; //!!!
else fVecBGlo[ iIDg ] -= weight*resid*derGlo[ig]; //!!!
//
int nfill = 0;
for (int jg=ig+1;jg--;) { // matCGlo is symmetric by construction
int jIDg = indGlo[jg];
- if ( fSigmaPar[jIDg] <= 0.) continue; // fixed parameter RRRCheck
+ if ( jIDg<0 || fSigmaPar[jIDg] <= 0.) continue; // fixed parameter RRRCheck
if ( !IsZero(vl = weight*derGlo[ig]*derGlo[jg]) ) {
fFillIndex[nfill] = jIDg;
fFillValue[nfill++] = fLocFitAdd ? vl:-vl;
AliInfo(Form("Global fit %s, CPU time: %.1f",res ? "Converged":"Failed",sw.CpuTime()));
if (!res) return 0;
//
- if (par) for (int i=fNGloPar;i--;) par[i] = fInitPar[i]+fDeltaPar[i];
+ if (par) for (int i=fNGloParIni;i--;) par[i] = GetFinalParam(i);
//
if (fGloSolveStatus==kInvert) { // errors on params are available
- if (error) for (int i=fNGloPar;i--;) error[i] = fProcPnt[i]>0 ? TMath::Sqrt(TMath::Abs(fMatCGlo->QueryDiag(i))) : 0.;
- if (pull) for (int i=fNGloPar;i--;) pull[i] = fProcPnt[i]>0 && (fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QueryDiag(i))>0. && fSigmaPar[i]>0
- ? fDeltaPar[i]/TMath::Sqrt(fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QueryDiag(i)) : 0;
+ if (error) for (int i=fNGloParIni;i--;) error[i] = GetFinalError(i);
+ if (pull) for (int i=fNGloParIni;i--;) pull[i] = GetPull(i);
}
//
return 1;
double *der = fRecord->GetValue()+2;
int csize = fRecord->GetSize()-2;
//
+ if (fkReGroup) {
+ for (int jp=csize;jp--;) {
+ int idp = indV[jp];
+ if (fkReGroup[idp]<0) AliFatal(Form("Constain is requested for suppressed parameter #%d",indV[jp]));
+ indV[jp] = idp;
+ }
+ }
// check if after suppression of fixed variables there are non-0 derivatives
// and determine the max statistics of involved params
int nSuppressed = 0;
//
sws.Start();
+#ifdef _DUMP_EQ_BEFORE_
+ const char* faildumpB = Form("mp2eq_before%d.dat",fIter);
+ int defoutB = dup(1);
+ if (defoutB<0) {AliFatal("Failed on dup"); exit(1);}
+ int slvDumpB = open(faildumpB, O_RDWR|O_CREAT, 0666);
+ if (slvDumpB>=0) {
+ dup2(slvDumpB,1);
+ printf("Solving%d for %d params\n",fIter,fNGloSize);
+ matCGlo.Print("10");
+ for (int i=0;i<fNGloSize;i++) printf("b%2d : %+.10f\n",i,fVecBGlo[i]);
+ }
+ dup2(defoutB,1);
+ close(slvDumpB);
+ close(defoutB);
+
+#endif
/*
printf("Solving:\n");
- matCGlo.Print();
+ matCGlo.Print("l");
for (int i=0;i<fNGloSize;i++) printf("b%2d : %+e\n",i,fVecBGlo[i]);
*/
+#ifdef _DUMPEQ_BEFORE_
+ const char* faildumpB = Form("mp2eq_before%d.dat",fIter);
+ int defoutB = dup(1);
+ int slvDumpB = open(faildumpB, O_RDWR|O_CREAT, 0666);
+ dup2(slvDumpB,1);
+ //
+ printf("#Equation before step %d\n",fIter);
+ fMatCGlo->Print("10");
+ printf("#RHS/STAT : NGlo:%d NGloSize:%d\n",fNGloPar,fNGloSize);
+ for (int i=0;i<fNGloSize;i++) printf("%d %+.10f %d\n",i,fVecBGlo[i],fProcPnt[i]);
+ //
+ dup2(defoutB,1);
+ close(slvDumpB);
+ close(defoutB);
+#endif
+ //
fGloSolveStatus = SolveGlobalMatEq(); // obtain solution for this step
+#ifdef _DUMPEQ_AFTER_
+ const char* faildumpA = Form("mp2eq_after%d.dat",fIter);
+ int defoutA = dup(1);
+ int slvDumpA = open(faildumpA, O_RDWR|O_CREAT, 0666);
+ dup2(slvDumpA,1);
+ //
+ printf("#Matrix after step %d\n",fIter);
+ fMatCGlo->Print("10");
+ printf("#RHS/STAT : NGlo:%d NGloSize:%d\n",fNGloPar,fNGloSize);
+ for (int i=0;i<fNGloSize;i++) printf("%d %+.10f %d\n",i,fVecBGlo[i],fProcPnt[i]);
+ //
+ dup2(defoutA,1);
+ close(slvDumpA);
+ close(defoutA);
+#endif
+ //
sws.Stop();
printf("Solve %d |",fIter); sws.Print();
//
if (fGloSolveStatus==kFailed) return 0;
//
for (int i=fNGloPar;i--;) fDeltaPar[i] += fVecBGlo[i]; // Update global parameters values (for iterations)
+
+#ifdef _DUMP_EQ_AFTER_
+ const char* faildumpA = Form("mp2eq_after%d.dat",fIter);
+ int defoutA = dup(1);
+ if (defoutA<0) {AliFatal("Failed on dup"); exit(1);}
+ int slvDumpA = open(faildumpA, O_RDWR|O_CREAT, 0666);
+ if (slvDumpA>=0) {
+ dup2(slvDumpA,1);
+ printf("Solving%d for %d params\n",fIter,fNGloSize);
+ matCGlo.Print("10");
+ for (int i=0;i<fNGloSize;i++) printf("b%2d : %+.10f\n",i,fVecBGlo[i]);
+ }
+ dup2(defoutA,1);
+ close(slvDumpA);
+ close(defoutA);
+#endif
//
- // PrintGlobalParameters();
+ /*
+ printf("Solved:\n");
+ matCGlo.Print("l");
+ for (int i=0;i<fNGloSize;i++) printf("b%2d : %+e (->%+e)\n",i,fVecBGlo[i], fDeltaPar[i]);
+ */
+
+ PrintGlobalParameters();
return 1;
}
//
/*
printf("GlobalMatrix\n");
- fMatCGlo->Print();
+ fMatCGlo->Print("l");
printf("RHS\n");
for (int i=0;i<fNGloPar;i++) printf("%d %+e\n",i,fVecBGlo[i]);
*/
//
dup2(defout,1);
close(slvDump);
+ close(defout);
printf("#Dumped failed matrix and RHS to %s\n",faildump);
}
else AliInfo("Failed on file open for matrix dumping");
return kFailed;
}
for (int i=fNGloSize;i--;) fVecBGlo[i] = sol[i];
+ //
return kNoInversion;
//
}
{
// return error for parameter iPar
if (fGloSolveStatus==kInvert) {
+ if (fkReGroup) iPar = fkReGroup[iPar];
+ if (iPar<0) {
+ // AliDebug(2,Form("Parameter %d was suppressed in the regrouping",iPar));
+ return 0;
+ }
double res = fMatCGlo->QueryDiag(iPar);
if (res>=0) return TMath::Sqrt(res);
}
- return -1.;
+ return 0.;
+}
+
+//_____________________________________________________________________________________________
+Double_t AliMillePede2::GetPull(int iPar) const
+{
+ // return pull for parameter iPar
+ if (fGloSolveStatus==kInvert) {
+ if (fkReGroup) iPar = fkReGroup[iPar];
+ if (iPar<0) {
+ // AliDebug(2,Form("Parameter %d was suppressed in the regrouping",iPar));
+ return 0;
+ }
+ //
+ return fProcPnt[iPar]>0 && (fSigmaPar[iPar]*fSigmaPar[iPar]-fMatCGlo->QueryDiag(iPar))>0. && fSigmaPar[iPar]>0
+ ? fDeltaPar[iPar]/TMath::Sqrt(fSigmaPar[iPar]*fSigmaPar[iPar]-fMatCGlo->QueryDiag(iPar)) : 0;
+ }
+ return 0.;
}
double lError = 0.;
double lGlobalCor =0.;
- AliInfo("");
- AliInfo(" Result of fit for global parameters");
- AliInfo(" ===================================");
- AliInfo(" I initial final differ lastcor error gcor Npnt");
- AliInfo("----------------------------------------------------------------------------------------------");
- //
- for (int i=0; i<fNGloPar; i++) {
- lError = GetParError(i);
+ printf("\nMillePede2 output\n");
+ printf(" Result of fit for global parameters\n");
+ printf(" ===================================\n");
+ printf(" I initial final differ lastcor error gcor Npnt\n");
+ printf("----------------------------------------------------------------------------------------------\n");
+ //
+ int lastPrintedId = -1;
+ for (int i0=0; i0<fNGloParIni; i0++) {
+ int i = GetRGId(i0); if (i<0) continue;
+ if (i!=i0 && lastPrintedId>=0 && i<=lastPrintedId) continue; // grouped param
+ lastPrintedId = i;
+ lError = GetParError(i0);
lGlobalCor = 0.0;
//
double dg;
if (fGloSolveStatus==kInvert && TMath::Abs( (dg=fMatCGlo->QueryDiag(i)) *fDiagCGlo[i]) > 0) {
lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(dg*fDiagCGlo[i])));
- AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %6d",
- i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor,fProcPnt[i]));
+ printf("%4d(%4d)\t %+.6f\t %+.6f\t %+.6f\t %.6f\t %.6f\t %.6f\t %6d\n",
+ i,i0,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor,fProcPnt[i]);
}
else {
- AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t OFF\t OFF\t %6d",i,fInitPar[i],fInitPar[i]+fDeltaPar[i],
- fDeltaPar[i],fVecBGlo[i],fProcPnt[i]));
+ printf("%4d (%4d)\t %+.6f\t %+.6f\t %+.6f\t %.6f\t OFF\t OFF\t %6d\n",i,i0,fInitPar[i],fInitPar[i]+fDeltaPar[i],
+ fDeltaPar[i],fVecBGlo[i],fProcPnt[i]);
}
}
return 1;
}
//_____________________________________________________________________________________________
-Bool_t AliMillePede2::IsRecordAcceptable() const
+Bool_t AliMillePede2::IsRecordAcceptable()
{
// validate record according run lists set by the user
static Long_t prevRunID = kMaxInt;
Long_t runID = fRecord->GetRunID();
if (runID!=prevRunID) {
int n = 0;
+ fRunWgh = 1.;
prevRunID = runID;
// is run to be rejected?
if (fRejRunList && (n=fRejRunList->GetSize())) {
prevAns = kTRUE;
for (int i=n;i--;) if (runID == (*fRejRunList)[i]) {
prevAns = kFALSE;
- printf("New Run to reject: %ld -> %d\n",runID,prevAns);
+ AliInfo(Form("New Run to reject: %ld",runID));
break;
}
}
else if (fAccRunList && (n=fAccRunList->GetSize())) { // is run specifically selected
prevAns = kFALSE;
- for (int i=n;i--;) if (runID == (*fAccRunList)[i]) {prevAns = kTRUE; break;}
+ for (int i=n;i--;) if (runID == (*fAccRunList)[i]) {
+ prevAns = kTRUE;
+ if (fAccRunListWgh) fRunWgh = (*fAccRunListWgh)[i];
+ AliInfo(Form("New Run to accept explicitly: %ld, weight=%f",runID,fRunWgh));
+ break;
+ }
+ if (!prevAns) AliInfo(Form("New Run is not in the list to accept: %ld",runID));
}
}
//
}
//_____________________________________________________________________________________________
-void AliMillePede2::SetAccRunList(const UInt_t *runs, Int_t nruns)
+void AliMillePede2::SetAccRunList(const UInt_t *runs, Int_t nruns, const Float_t* wghList)
{
// set the list of runs to be selected
if (fAccRunList) delete fAccRunList;
+ if (fAccRunListWgh) delete fAccRunListWgh;
fAccRunList = 0;
if (nruns<1 || !runs) return;
fAccRunList = new TArrayL(nruns);
- for (int i=0;i<nruns;i++) (*fAccRunList)[i] = runs[i];
+ fAccRunListWgh = new TArrayF(nruns);
+ for (int i=0;i<nruns;i++) {
+ (*fAccRunList)[i] = runs[i];
+ (*fAccRunListWgh)[i] =wghList ? wghList[i] : 1.0;
+ }
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetInitPars(const Double_t* par)
+{
+ // initialize parameters, account for eventual grouping
+ for (int i=0;i<fNGloParIni;i++) {
+ int id = GetRGId(i); if (id<0) continue;
+ fInitPar[id] = par[i];
+ }
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetSigmaPars(const Double_t* par)
+{
+ // initialize sigmas, account for eventual grouping
+ for (int i=0;i<fNGloParIni;i++) {
+ int id = GetRGId(i); if (id<0) continue;
+ fSigmaPar[id] = par[i];
+ }
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetInitPar(Int_t i,Double_t par)
+{
+ // initialize param, account for eventual grouping
+ int id = GetRGId(i); if (id<0) return;
+ fInitPar[id] = par;
}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetSigmaPar(Int_t i,Double_t par)
+{
+ // initialize sigma, account for eventual grouping
+ int id = GetRGId(i); if (id<0) return;
+ fSigmaPar[id] = par;
+}
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff