#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;
fConsRecFile(0),
fCurrRecDataID(0),
fCurrRecConstrID(0),
- fLocFitAdd(kTRUE)
+ fLocFitAdd(kTRUE),
+ fSelFirst(1),
+ fSelLast(-1),
+ fRejRunList(0),
+ fAccRunList(0)
{}
//_____________________________________________________________________________________________
fMatCLoc(0),fMatCGlo(0),fMatCGloLoc(0),fFillIndex(0),fFillValue(0),
fDataRecFName(0),fRecord(0),fDataRecFile(0),
fTreeData(0),fRecFileStatus(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),
- fCurrRecConstrID(0),fLocFitAdd(kTRUE)
+ 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) {AliFatal(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) {AliFatal(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;
+ recfName.ReadLine(inpf);
+ 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;
//
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;
//
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)
{
//
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, trying Gaussian Elimination");
}
//
// 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)));
}
} // 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);
//
//
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();
fNGloFix++;
fVecBGlo[i] = 0.;
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