#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
fCurrRecDataID(0),
fCurrRecConstrID(0),
fLocFitAdd(kTRUE),
+ fUseRecordWeight(kTRUE),
+ fMinRecordLength(1),
fSelFirst(1),
fSelLast(-1),
fRejRunList(0),
fAccRunList(0),
+ fAccRunListWgh(0),
+ fRunWgh(1),
fkReGroup(0)
-{}
+{
+ fWghScl[0] = fWghScl[1] = -1;
+}
//_____________________________________________________________________________________________
AliMillePede2::AliMillePede2(const AliMillePede2& src) :
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),
+ fAccRunListWgh(0),
+ fRunWgh(1),
fkReGroup(0)
-{printf("Dummy\n");}
+{
+ fWghScl[0] = src.fWghScl[0];
+ fWghScl[1] = src.fWghScl[1];
+ printf("Dummy\n");
+}
//_____________________________________________________________________________________________
AliMillePede2::~AliMillePede2()
delete fMatCGloLoc;
delete fRejRunList;
delete fAccRunList;
+ delete fAccRunListWgh;
}
//_____________________________________________________________________________________________
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];}
- AliInfo(Form("Regrouping is requested: from %d raw to %d free globals, max global ID=%d",nGlo,ng,maxPID));
- if (ng != maxPID+1) AliFatal("Wrong regrouping requested");
- nGlo = ng;
+ 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;
fParamGrID[i] = -1;
}
//
+ fWghScl[0] = -1;
+ fWghScl[1] = -1;
return 1;
}
//
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];
//
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 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];
//
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("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
//
/*
printf("Solved:\n");
for (int i=0;i<fNGloSize;i++) printf("b%2d : %+e (->%+e)\n",i,fVecBGlo[i], fDeltaPar[i]);
*/
- // PrintGlobalParameters();
+ PrintGlobalParameters();
return 1;
}
//
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;}
+ 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 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;}
+ 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;
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("----------------------------------------------------------------------------------------------");
+ 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;
- lError = GetParError(i);
+ 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;
+ }
}
//_____________________________________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff