using namespace std;
+//#define _DUMPEQ_BEFORE_
+//#define _DUMPEQ_AFTER_
ClassImp(AliMillePede2)
fCurrRecDataID(0),
fCurrRecConstrID(0),
fLocFitAdd(kTRUE),
+ fUseRecordWeight(kTRUE),
+ fMinRecordLength(1),
fSelFirst(1),
fSelLast(-1),
fRejRunList(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),
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()*fRunWgh; // 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];
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);
+#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);
+#endif
+ //
sws.Stop();
printf("Solve %d |",fIter); sws.Print();
//
for (int i=0;i<fNGloSize;i++) printf("b%2d : %+e (->%+e)\n",i,fVecBGlo[i], fDeltaPar[i]);
*/
- // PrintGlobalParameters();
+ PrintGlobalParameters();
return 1;
}
return kFailed;
}
for (int i=fNGloSize;i--;) fVecBGlo[i] = sol[i];
+ //
return kNoInversion;
//
}
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;}
Long_t fCurrRecDataID; // ID of the current data record
Long_t fCurrRecConstrID; // ID of the current constraint record
Bool_t fLocFitAdd; // Add contribution of carrent track (and not eliminate it)
+ Bool_t fUseRecordWeight; // force or ignore the record weight
+ Int_t fMinRecordLength; // ignore shorter records
Int_t fSelFirst; // event selection start
Int_t fSelLast; // event selection end
TArrayL* fRejRunList; // list of runs to reject (if any)
TArrayL* fAccRunList; // list of runs to select (if any)
TArrayF* fAccRunListWgh; // optional weights for data of accepted runs (if any)
Double_t fRunWgh; // run weight
+ Double_t fWghScl[2]; // optional rescaling for odd/even residual weights (see its usage in LocalFit)
const Int_t* fkReGroup; // optional regrouping of parameters wrt ID's from the records
//
static Bool_t fgInvChol; // Invert global matrix in Cholesky solver