#include "AliLog.h"
#include <TStopwatch.h>
#include <TFile.h>
+#include <TChain.h>
#include <TMath.h>
#include <TVectorD.h>
+#include <TArrayL.h>
+#include <TSystem.h>
#include "AliMatrixSq.h"
#include "AliSymMatrix.h"
#include "AliRectMatrix.h"
#include "AliMatrixSparse.h"
-
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <fstream>
using namespace std;
fRecord(0),
fDataRecFile(0),
fTreeData(0),
+ fRecFileStatus(0),
//
fConstrRecFName("/tmp/mp2_constraints_records.root"),
fTreeConstr(0),
fConsRecFile(0),
fCurrRecDataID(0),
fCurrRecConstrID(0),
- fLocFitAdd(kTRUE)
+ fLocFitAdd(kTRUE),
+ fSelFirst(1),
+ fSelLast(-1),
+ fRejRunList(0),
+ fAccRunList(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),
fDataRecFName(0),fRecord(0),fDataRecFile(0),
- fTreeData(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),
- fCurrRecConstrID(0),fLocFitAdd(kTRUE)
+ fTreeData(0),fRecFileStatus(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),
+ fCurrRecConstrID(0),fLocFitAdd(kTRUE),
+ fSelFirst(1),
+ fSelLast(-1),
+ fRejRunList(0),
+ fAccRunList(0)
{printf("Dummy\n");}
//_____________________________________________________________________________________________
delete fMatCLoc;
delete fMatCGlo;
delete fMatCGloLoc;
+ delete fRejRunList;
+ delete fAccRunList;
}
//_____________________________________________________________________________________________
{
// initialize the buffer for processed measurements records
//
- if (fDataRecFile) {AliInfo("Data Records File is already initialized"); return kFALSE;}
+ if (fTreeData) {AliInfo("Data Records File is already initialized"); return kFALSE;}
//
if (!fRecord) fRecord = new AliMillePedeRecord();
//
- fDataRecFile = TFile::Open(GetDataRecFName(),read ? "":"recreate");
- if (!fDataRecFile) {AliInfo(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;}
- //
- AliInfo(Form("File %s used for derivatives records",GetDataRecFName()));
- if (read) {
- fTreeData = (TTree*)fDataRecFile->Get("AliMillePedeRecords_Data");
- if (!fTreeData) {AliInfo(Form("Did not find data records tree in %s",GetDataRecFName())); return kFALSE;}
- fTreeData->SetBranchAddress("Record_Data",&fRecord);
- AliInfo(Form("Found %d derivatives records",fTreeData->GetEntries()));
- }
- else {
+ if (!read) { // write mode: cannot use chain
+ 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);
}
+ else { // use chain
+ TChain* ch = new TChain("AliMillePedeRecords_Data");
+ //
+ if (fDataRecFName.EndsWith(".root")) ch->AddFile(fDataRecFName);
+ else { // assume text file with list of filenames
+ //
+ ifstream inpf(fDataRecFName.Data());
+ if (!inpf.good()) {AliInfo(Form("Failed on input records list %s\n",fDataRecFName.Data())); return kFALSE;}
+ //
+ TString recfName;
+ while ( !(recfName.ReadLine(inpf)).eof() ) {
+ recfName = recfName.Strip(TString::kBoth,' ');
+ if (recfName.BeginsWith("//") || recfName.BeginsWith("#") || !recfName.EndsWith(".root")) { // comment
+ AliInfo(Form("Skip %s\n",recfName.Data()));
+ continue;
+ }
+ //
+ recfName = recfName.Strip(TString::kBoth,',');
+ recfName = recfName.Strip(TString::kBoth,'"');
+ gSystem->ExpandPathName(recfName);
+ printf("Adding %s\n",recfName.Data());
+ ch->AddFile(recfName.Data());
+ }
+ }
+ //
+ Long64_t nent = ch->GetEntries();
+ if (nent<1) { AliInfo("Obtained chain is empty"); return kFALSE;}
+ fTreeData = ch;
+ fTreeData->SetBranchAddress("Record_Data",&fRecord);
+ AliInfo(Form("Found %lld derivatives records",nent));
+ }
fCurrRecDataID = -1;
+ fRecFileStatus = read ? 1:2;
//
return kTRUE;
}
fTreeConstr = (TTree*)fConsRecFile->Get("AliMillePedeRecords_Constraints");
if (!fTreeConstr) {AliInfo(Form("Did not find constraints records tree in %s",GetConsRecFName())); return kFALSE;}
fTreeConstr->SetBranchAddress("Record_Constraints",&fRecord);
- AliInfo(Form("Found %d constraints records",fTreeConstr->GetEntries()));
+ AliInfo(Form("Found %lld constraints records",fTreeConstr->GetEntries()));
//
}
else {
void AliMillePede2::CloseDataRecStorage()
{
if (fTreeData) {
- if (fDataRecFile->IsWritable()) {
+ if (fDataRecFile && fDataRecFile->IsWritable()) {
fDataRecFile->cd();
fTreeData->Write();
}
delete fTreeData;
fTreeData = 0;
- fDataRecFile->Close();
- delete fDataRecFile;
- fDataRecFile = 0;
+ if (fDataRecFile) {
+ fDataRecFile->Close();
+ delete fDataRecFile;
+ fDataRecFile = 0;
+ }
}
+ fRecFileStatus = 0;
//
}
return kTRUE;
}
+//_____________________________________________________________________________________________
+void AliMillePede2::SetRecordWeight(double wgh)
+{
+ if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data
+ fRecord->SetWeight(wgh);
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetRecordRun(Int_t run)
+{
+ if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data
+ fRecord->SetRunID(run);
+}
+
//_____________________________________________________________________________________________
void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma)
{
- if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data
+ if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data
//
// write data of single measurement
if (lSigma<=0.0) { // If parameter is fixed, then no equation
return;
}
//
- if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data
+ if (fRecFileStatus<2) InitDataRecStorage(); // create a buffer to store the data
//
fRecord->AddResidual(lMeas);
//
}
double vl;
//
+ double 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 );
+ double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
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)
+ for (int i=nrefGlo[ip];i--;) { // suppress the global part (only relevant with iterations)
int iID = indGlo[i]; // Global param indice
if (fSigmaPar[iID]<=0.) continue; // fixed parameter RRRCheck
if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]); // linear parameter
//
matCLoc.SetSizeUsed(++maxLocUsed); // data with B=0 may use less than declared nLocals
//
+ /* //RRR
+ fRecord->Print("l");
+ printf("\nBefore\nLocalMatrix: "); matCLoc.Print("l");
+ printf("RHSLoc: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n");
+ */
// first try to solve by faster Cholesky decomposition, then by Gaussian elimination
if (!matCLoc.SolveChol(fVecBLoc,kTRUE)) {
- AliInfo("Failed to solve locals by Cholesky, traying Gaussian Elimination");
+ AliInfo("Failed to solve locals by Cholesky, trying Gaussian Elimination");
if (!matCLoc.SolveSpmInv(fVecBLoc,kTRUE)) {
AliInfo("Failed to solve locals by Gaussian Elimination, skip...");
matCLoc.Print("d");
}
//
// If requested, store the track params and errors
- // printf("locfit: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n");
+ //RRR printf("locfit: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n");
+
if (localParams) for (int i=maxLocUsed; i--;) {
localParams[2*i] = fVecBLoc[i];
localParams[2*i+1] = TMath::Sqrt(TMath::Abs(matCLoc.QueryDiag(i)));
//
for (int ip=nPoints;ip--;) { // Calculate residuals
double resid = fRecord->GetValue( refLoc[ip]-1 );
- double weight = fRecord->GetValue( refGlo[ip]-1 );
+ double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
double *derLoc = fRecord->GetValue()+refLoc[ip];
double *derGlo = fRecord->GetValue()+refGlo[ip];
int *indLoc = fRecord->GetIndex()+refLoc[ip];
nEq++; // number of equations
} // end of Calculate residuals
//
+ lChi2 /= gloWgh;
int nDoF = nEq-maxLocUsed;
lChi2 = (nDoF>0) ? lChi2/nDoF : 0; // Chi^2/dof
//
//
for (int ip=nPoints;ip--;) { // Update matrices
double resid = fRecord->GetValue( refLoc[ip]-1 );
- double weight = fRecord->GetValue( refGlo[ip]-1 );
+ double weight = fRecord->GetValue( refGlo[ip]-1 )*gloWgh;
double *derLoc = fRecord->GetValue()+refLoc[ip];
double *derGlo = fRecord->GetValue()+refGlo[ip];
int *indLoc = fRecord->GetIndex()+refLoc[ip];
}
} // end of Update matrices
//
+ /*//RRR
+ printf("After GLO\n");
+ printf("MatCLoc: "); fMatCLoc->Print("l");
+ printf("MatCGlo: "); fMatCGlo->Print("l");
+ printf("MatCGlLc:"); fMatCGloLoc->Print("l");
+ printf("BGlo: "); for (int i=0; i<fNGloPar; i++) printf("%+e |",fVecBGlo[i]); printf("\n");
+ */
// calculate fMatCGlo -= fMatCGloLoc * fMatCLoc * fMatCGloLoc^T
// and fVecBGlo -= fMatCGloLoc * fVecBLoc
//
//
// reset compressed index array
//
+ /*//RRR
+ printf("After GLOLoc\n");
+ printf("MatCGlo: "); fMatCGlo->Print("");
+ printf("BGlo: "); for (int i=0; i<fNGloPar; i++) printf("%+e |",fVecBGlo[i]); printf("\n");
+ */
for (int i=nGloInFit;i--;) {
fGlo2CGlo[ fCGlo2Glo[i] ] = -1;
fCGlo2Glo[i] = -1;
fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
if (fChi2CutFactor < 1.2*fChi2CutRef) {
fChi2CutFactor = fChi2CutRef;
- fIter = fMaxIter - 1; // Last iteration
+ //RRR fIter = fMaxIter - 1; // Last iteration
}
}
fIter++;
//
// Process data records and build the matrices
Long_t ndr = fTreeData->GetEntries();
- AliInfo(Form("Building the Global matrix from %d data records",ndr));
- if (!ndr) return 0;
+ Long_t first = fSelFirst>0 ? fSelFirst : 0;
+ Long_t last = fSelLast<1 ? ndr : (fSelLast>=ndr ? ndr : fSelLast+Long_t(1));
+ ndr = last - first;
+ //
+ AliInfo(Form("Building the Global matrix from data records %ld : %ld",first,last));
+ if (ndr<1) return 0;
//
TStopwatch swt; swt.Start();
fLocFitAdd = kTRUE; // add contributions of matching tracks
for (Long_t i=0;i<ndr;i++) {
- ReadRecordData(i);
+ Long_t iev = i+first;
+ ReadRecordData(iev);
+ if (!IsRecordAcceptable()) continue;
LocalFit();
if ( (i%int(0.2*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr);
}
swt.Stop();
printf("%ld local fits done: ", ndr);
+ /*
+ printf("MatCGlo: "); fMatCGlo->Print("l");
+ printf("BGlo: "); for (int i=0; i<fNGloPar; i++) printf("%+e |",fVecBGlo[i]); printf("\n");
swt.Print();
+ */
sw.Start(kFALSE);
//
//
Int_t nFixedGroups = 0;
TArrayI fixGroups(fixArrSize);
//
+ int grIDold = -2;
+ int oldStart = -1;
+ double oldMin = 1.e20;
+ double oldMax =-1.e20;
+ //
for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones
int grID = fParamGrID[i];
if (grID<0) continue; // not in the group
- if (fProcPnt[i]>=fMinPntValid) continue;
- fProcPnt[i] = 0;
- // check if the group is already accounted
- Bool_t fnd = kFALSE;
- for (int j=nFixedGroups;j--;) if (fixGroups[j]==grID) {fnd=kTRUE; break;}
- if (fnd) continue;
- if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);}
- fixGroups[nFixedGroups++] = grID; // add group to fix
//
+ if (grID!=grIDold) { // starting new group
+ if (grIDold>=0) { // decide if the group has enough statistics
+ if (oldMin<fMinPntValid && oldMax<2*fMinPntValid) { // suppress group
+ for (int iold=oldStart;iold>i;iold--) fProcPnt[iold] = 0;
+ Bool_t fnd = kFALSE; // check if the group is already accounted
+ for (int j=nFixedGroups;j--;) if (fixGroups[j]==grIDold) {fnd=kTRUE; break;}
+ if (!fnd) {
+ if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);}
+ fixGroups[nFixedGroups++] = grIDold; // add group to fix
+ }
+ }
+ }
+ grIDold = grID; // mark the start of the new group
+ oldStart = i;
+ oldMin = 1.e20;
+ oldMax = -1.e20;
+ }
+ if (oldMin>fProcPnt[i]) oldMin = fProcPnt[i];
+ if (oldMax<fProcPnt[i]) oldMax = fProcPnt[i];
+ //
+ }
+ // extra check for the last group
+ if (grIDold>=0 && oldMin<fMinPntValid && oldMax<2*fMinPntValid) { // suppress group
+ for (int iold=oldStart;iold--;) fProcPnt[iold] = 0;
+ Bool_t fnd = kFALSE; // check if the group is already accounted
+ for (int j=nFixedGroups;j--;) if (fixGroups[j]==grIDold) {fnd=kTRUE; break;}
+ if (!fnd) {
+ if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);}
+ fixGroups[nFixedGroups++] = grIDold; // add group to fix
+ }
}
+ //
// 2) loop over records and add contributions of fixed groups with negative sign
fLocFitAdd = kFALSE;
//
for (Long_t i=0;i<ndr;i++) {
- ReadRecordData(i);
+ Long_t iev = i+first;
+ ReadRecordData(iev);
+ if (!IsRecordAcceptable()) continue;
Bool_t suppr = kFALSE;
for (int ifx=nFixedGroups;ifx--;)if (fRecord->IsGroupPresent(fixGroups[ifx])) suppr = kTRUE;
if (suppr) LocalFit();
if (fProcPnt[i]<1) {
fNGloFix++;
fVecBGlo[i] = 0.;
- matCGlo.DiagElem(i) = float(fNLocEquations*fNLocEquations);
+ matCGlo.DiagElem(i) = 1.;//float(fNLocEquations*fNLocEquations);
+ // matCGlo.DiagElem(i) = float(fNLocEquations*fNLocEquations);
}
else matCGlo.DiagElem(i) += (fgWeightSigma ? fProcPnt[i] : 1.)/(fSigmaPar[i]*fSigmaPar[i]);
}
//
sws.Start();
+
+ /*
+ printf("Solving:\n");
+ matCGlo.Print();
+ for (int i=0;i<fNGloSize;i++) printf("b%2d : %+e\n",i,fVecBGlo[i]);
+ */
fGloSolveStatus = SolveGlobalMatEq(); // obtain solution for this step
sws.Stop();
printf("Solve %d |",fIter); sws.Print();
}
//
if (((AliSymMatrix*)fMatCGlo)->SolveSpmInv(fVecBGlo, kTRUE)) return gkInvert;
- else AliInfo("Solution of Global Dense System by Gaussian Elimiation failed, trying iterative methods");
+ else AliInfo("Solution of Global Dense System by Gaussian Elimination failed, trying iterative methods");
}
// try to solve by minres
TVectorD sol(fNGloSize);
else
AliInfo(Form("Undefined Iteritive Solver ID=%d, only %d are defined",fgIterSol,AliMinResSolve::kNSolvers));
//
- if (!res) return gkFailed;
+ if (!res) {
+ const char* faildump = "fgmr_failed.dat";
+ int defout = dup(1);
+ if (defout<0) {
+ AliInfo("Failed on dup");
+ return gkFailed;
+ }
+ int slvDump = open(faildump, O_RDWR|O_CREAT, 0666);
+ if (slvDump>=0) {
+ dup2(slvDump,1);
+ //
+ printf("#Failed to solve using solver %d with PreCond: %d MaxIter: %d Tol: %e NKrylov: %d\n",
+ fgIterSol,fgMinResCondType,fgMinResMaxIter,fgMinResTol,fgNKrylovV);
+ printf("#Dump of matrix:\n");
+ fMatCGlo->Print("10");
+ printf("#Dump of RHS:\n");
+ for (int i=0;i<fNGloSize;i++) printf("%d %+.10f\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 gkFailed;
+ }
for (int i=fNGloSize;i--;) fVecBGlo[i] = sol[i];
return gkNoInversion;
//
}
return 1;
}
+
+//_____________________________________________________________________________________________
+Bool_t AliMillePede2::IsRecordAcceptable() const
+{
+ // validate record according run lists set by the user
+ static Long_t prevRunID = kMaxInt;
+ static Bool_t prevAns = kTRUE;
+ Long_t runID = fRecord->GetRunID();
+ if (runID!=prevRunID) {
+ int n = 0;
+ 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);
+ 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;}
+ }
+ }
+ //
+ return prevAns;
+ //
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetRejRunList(const UInt_t *runs, Int_t nruns)
+{
+ // set the list of runs to be rejected
+ if (fRejRunList) delete fRejRunList;
+ fRejRunList = 0;
+ if (nruns<1 || !runs) return;
+ fRejRunList = new TArrayL(nruns);
+ for (int i=0;i<nruns;i++) (*fRejRunList)[i] = runs[i];
+}
+
+//_____________________________________________________________________________________________
+void AliMillePede2::SetAccRunList(const UInt_t *runs, Int_t nruns)
+{
+ // set the list of runs to be selected
+ if (fAccRunList) delete fAccRunList;
+ fAccRunList = 0;
+ if (nruns<1 || !runs) return;
+ fAccRunList = new TArrayL(nruns);
+ for (int i=0;i<nruns;i++) (*fAccRunList)[i] = runs[i];
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff