- Three classes by MinJung Kweon AliHFEpriVtx, AliHFEsecVtx and AliHFEmcQA for primar...
[u/mrichter/AliRoot.git] / STEER / AliMillePede2.cxx
1 /**********************************************************************************************/
2 /* General class for alignment with large number of degrees of freedom                        */
3 /* Based on the original milliped2 by Volker Blobel                                           */
4 /* http://www.desy.de/~blobel/mptalks.html                                                    */
5 /*                                                                                            */ 
6 /* Author: ruben.shahoyan@cern.ch                                                             */
7 /*                                                                                            */ 
8 /**********************************************************************************************/
9
10 #include "AliMillePede2.h"
11 #include "AliLog.h"
12 #include <TStopwatch.h>
13 #include <TFile.h>
14 #include <TMath.h>
15 #include "AliMatrixSq.h"
16 #include "AliSymMatrix.h"
17 #include "AliRectMatrix.h"
18 #include "AliMatrixSparse.h"
19
20
21 using namespace std;
22
23
24 ClassImp(AliMillePede2)
25
26 Bool_t   AliMillePede2::fgInvChol        = kTRUE;     // Invert global matrix with Cholesky solver
27 Bool_t   AliMillePede2::fgWeightSigma    = kTRUE;     // weight local constraint by module statistics
28 Bool_t   AliMillePede2::fgIsMatGloSparse = kFALSE;    // use faster dense matrix by default
29 Int_t    AliMillePede2::fgMinResCondType = 1;         // Jacoby preconditioner by default
30 Double_t AliMillePede2::fgMinResTol      = 1.e-11;    // default tolerance
31 Int_t    AliMillePede2::fgMinResMaxIter  = 10000;     // default max number of iterations 
32 Int_t    AliMillePede2::fgIterSol        = AliMinResSolve::kSolMinRes; // default iterative solver
33 Int_t    AliMillePede2::fgNKrylovV       = 240;        // default number of Krylov vectors to keep
34
35 //_____________________________________________________________________________________________
36 AliMillePede2::AliMillePede2() 
37 : fNLocPar(0),
38   fNGloPar(0),
39   fNGloSize(0),
40 //
41   fNLocEquations(0),
42   fIter(0),
43   fMaxIter(10),
44   fNStdDev(3),
45   fNGloConstraints(0),
46   fNLagrangeConstraints(0),
47   fNLocFits(0),
48   fNLocFitsRejected(0),
49   fNGloFix(0),
50   fGloSolveStatus(gkFailed),
51 //
52   fChi2CutFactor(1.),
53   fChi2CutRef(1.),
54   fResCutInit(100.),
55   fResCut(100.),
56   fMinPntValid(1),
57 //
58   fNGroupsSet(0),
59   fParamGrID(0),
60   fProcPnt(0),
61   fVecBLoc(0),
62   fDiagCGlo(0),
63   fVecBGlo(0),
64   fInitPar(0),
65   fDeltaPar(0),
66   fSigmaPar(0),
67   fIsLinear(0),
68   fConstrUsed(0),
69 //
70   fGlo2CGlo(0),
71   fCGlo2Glo(0),
72 //
73   fMatCLoc(0),
74   fMatCGlo(0),
75   fMatCGloLoc(0),
76   //
77   fFillIndex(0),
78   fFillValue(0),
79   //
80   fDataRecFName("/tmp/mp2_data_records.root"),
81   fRecord(0),
82   fDataRecFile(0),  
83   fTreeData(0),
84   //
85   fConstrRecFName("/tmp/mp2_constraints_records.root"),
86   fTreeConstr(0),
87   fConsRecFile(0),
88   fCurrRecDataID(0),
89   fCurrRecConstrID(0),
90   fLocFitAdd(kTRUE)
91 {}
92
93 //_____________________________________________________________________________________________
94 AliMillePede2::AliMillePede2(const AliMillePede2& src) : 
95   TObject(src),fNLocPar(0),fNGloPar(0),fNGloSize(0),fNLocEquations(0),fIter(0),
96   fMaxIter(10),fNStdDev(3),fNGloConstraints(0),fNLagrangeConstraints(0),
97   fNLocFits(0),fNLocFitsRejected(0),
98   fNGloFix(0),fGloSolveStatus(0),fChi2CutFactor(0),fChi2CutRef(0),fResCutInit(0),
99   fResCut(0),fMinPntValid(1),fNGroupsSet(0),fParamGrID(0),fProcPnt(0),fVecBLoc(0),fDiagCGlo(0),fVecBGlo(0),
100   fInitPar(0),fDeltaPar(0),fSigmaPar(0),fIsLinear(0),fConstrUsed(0),fGlo2CGlo(0),fCGlo2Glo(0),
101   fMatCLoc(0),fMatCGlo(0),fMatCGloLoc(0),fFillIndex(0),fFillValue(0),
102   fDataRecFName(0),fRecord(0),fDataRecFile(0),
103   fTreeData(0),fConstrRecFName(0),fTreeConstr(0),fConsRecFile(0),fCurrRecDataID(0),
104   fCurrRecConstrID(0),fLocFitAdd(kTRUE)
105 {printf("Dummy\n");}
106
107 //_____________________________________________________________________________________________
108 AliMillePede2::~AliMillePede2() 
109 {
110   CloseDataRecStorage();
111   CloseConsRecStorage();
112   //
113   delete[] fParamGrID;
114   delete[] fProcPnt;
115   delete[] fVecBLoc;
116   delete[] fDiagCGlo;
117   delete[] fVecBGlo;
118   delete[] fInitPar;
119   delete[] fDeltaPar;
120   delete[] fSigmaPar;
121   delete[] fGlo2CGlo;
122   delete[] fCGlo2Glo;
123   delete[] fIsLinear;
124   delete[] fConstrUsed;
125   delete[] fFillIndex;
126   delete[] fFillValue;
127   //
128   delete fRecord;
129   delete fMatCLoc;
130   delete fMatCGlo;
131   delete fMatCGloLoc;
132
133
134 //_____________________________________________________________________________________________
135 Int_t AliMillePede2::InitMille(int nGlo, int nLoc, int lNStdDev,double lResCut, double lResCutInit)
136 {
137   //
138   if (nLoc>0)        fNLocPar = nLoc;
139   if (nGlo>0)        fNGloPar = nGlo;
140   if (lResCutInit>0) fResCutInit = lResCutInit; 
141   if (lResCut>0)     fResCut     = lResCut;
142   if (lNStdDev>0)    fNStdDev    = lNStdDev;
143   //
144   fNGloSize = fNGloPar;
145   //
146   try {
147     //
148     if (fgIsMatGloSparse) {fMatCGlo = new AliMatrixSparse(fNGloPar); fMatCGlo->SetSymmetric(kTRUE);}
149     else                   fMatCGlo = new AliSymMatrix(fNGloPar);
150     //
151     fFillIndex    = new Int_t[fNGloPar];
152     fFillValue    = new Double_t[fNGloPar];
153     //
154     fMatCLoc      = new AliSymMatrix(fNLocPar);
155     fMatCGloLoc   = new AliRectMatrix(fNGloPar,fNLocPar);
156     //
157     fParamGrID    = new Int_t[fNGloPar];
158     fProcPnt      = new Int_t[fNGloPar];
159     fVecBLoc      = new Double_t[fNLocPar];
160     fDiagCGlo     = new Double_t[fNGloPar];
161     //
162     fInitPar      = new Double_t[fNGloPar];
163     fDeltaPar     = new Double_t[fNGloPar];
164     fSigmaPar     = new Double_t[fNGloPar];
165     fIsLinear     = new Bool_t[fNGloPar];
166     //
167     fGlo2CGlo     = new Int_t[fNGloPar];
168     fCGlo2Glo     = new Int_t[fNGloPar];
169   }
170   catch(bad_alloc&) {
171     AliInfo(Form("Failed to allocate the memory for %d global and %d local parameters",fNGloPar,fNLocPar));
172     return 0;
173   }
174   //
175   memset(fVecBLoc   ,0,fNLocPar*sizeof(Double_t));
176   memset(fDiagCGlo  ,0,fNGloPar*sizeof(Double_t));
177   memset(fInitPar   ,0,fNGloPar*sizeof(Double_t));
178   memset(fDeltaPar  ,0,fNGloPar*sizeof(Double_t));
179   memset(fSigmaPar  ,0,fNGloPar*sizeof(Double_t));
180   memset(fProcPnt   ,0,fNGloPar*sizeof(Int_t));
181   //
182   for (int i=fNGloPar;i--;) {
183     fGlo2CGlo[i] = fCGlo2Glo[i] = -1;
184     fIsLinear[i] = kTRUE;
185     fParamGrID[i] = -1;
186   }
187   //
188   return 1;
189 }
190
191 //_____________________________________________________________________________________________
192 Bool_t AliMillePede2::ImposeDataRecFile(const char* fname)
193 {
194   CloseDataRecStorage();
195   SetDataRecFName(fname);
196   return InitDataRecStorage(kTRUE); // open in read mode
197 }
198
199 //_____________________________________________________________________________________________
200 Bool_t AliMillePede2::ImposeConsRecFile(const char* fname)
201 {
202   CloseConsRecStorage();
203   SetConsRecFName(fname);
204   return InitConsRecStorage(kTRUE); // open in read mode
205 }
206
207 //_____________________________________________________________________________________________
208 Bool_t AliMillePede2::InitDataRecStorage(Bool_t read)
209 {
210   // initialize the buffer for processed measurements records
211   //
212   if (fDataRecFile) {AliInfo("Data Records File is already initialized"); return kFALSE;} 
213   //
214   if (!fRecord) fRecord = new AliMillePedeRecord();
215   //
216   fDataRecFile = TFile::Open(GetDataRecFName(),read ? "":"recreate");
217   if (!fDataRecFile) {AliInfo(Form("Failed to initialize data records file %s",GetDataRecFName())); return kFALSE;}
218   //
219   AliInfo(Form("File %s used for derivatives records",GetDataRecFName()));
220   if (read) {
221     fTreeData = (TTree*)fDataRecFile->Get("AliMillePedeRecords_Data");
222     if (!fTreeData) {AliInfo(Form("Did not find data records tree in %s",GetDataRecFName())); return kFALSE;}
223     fTreeData->SetBranchAddress("Record_Data",&fRecord);
224     AliInfo(Form("Found %d derivatives records",fTreeData->GetEntries()));
225   }
226   else {
227     fTreeData = new TTree("AliMillePedeRecords_Data","Data Records for AliMillePede2");
228     fTreeData->Branch("Record_Data","AliMillePedeRecord",&fRecord,32000,99);
229   }
230   fCurrRecDataID = -1;
231   //
232   return kTRUE;
233 }
234
235 //_____________________________________________________________________________________________
236 Bool_t AliMillePede2::InitConsRecStorage(Bool_t read)
237 {
238   // initialize the buffer for processed measurements records
239   //
240   if (fConsRecFile) {AliInfo("Constraints Records File is already initialized"); return kFALSE;} 
241   //
242   if (!fRecord) fRecord = new AliMillePedeRecord();
243   //
244   fConsRecFile = TFile::Open(GetConsRecFName(),read ? "":"recreate");
245   if (!fConsRecFile) {AliInfo(Form("Failed to initialize constraints records file %s",GetConsRecFName())); return kFALSE;}
246   //
247   AliInfo(Form("File %s used for constraints records",GetConsRecFName()));
248   if (read) {
249     fTreeConstr = (TTree*)fConsRecFile->Get("AliMillePedeRecords_Constraints");
250     if (!fTreeConstr) {AliInfo(Form("Did not find constraints records tree in %s",GetConsRecFName())); return kFALSE;}
251     fTreeConstr->SetBranchAddress("Record_Constraints",&fRecord);
252     AliInfo(Form("Found %d constraints records",fTreeConstr->GetEntries()));
253     //
254   }
255   else {
256     //
257     fTreeConstr = new TTree("AliMillePedeRecords_Constraints","Constraints Records for AliMillePede2");
258     fTreeConstr->Branch("Record_Constraints","AliMillePedeRecord",&fRecord,32000,99);
259   }
260   fCurrRecConstrID = -1;
261   //
262   return kTRUE;
263 }
264
265 //_____________________________________________________________________________________________
266 void AliMillePede2::CloseDataRecStorage()
267 {
268   if (fTreeData) {
269     if (fDataRecFile->IsWritable()) {
270       fDataRecFile->cd();
271       fTreeData->Write();
272     }
273     delete fTreeData;  
274     fTreeData = 0;
275     fDataRecFile->Close();
276     delete fDataRecFile;
277     fDataRecFile = 0;
278   }
279   //
280 }
281
282 //_____________________________________________________________________________________________
283 void AliMillePede2::CloseConsRecStorage()
284 {
285   if (fTreeConstr) {
286     if (fConsRecFile->IsWritable()) {
287       fConsRecFile->cd();
288       fTreeConstr->Write();
289     }
290     delete fTreeConstr;  
291     fTreeConstr = 0;
292     fConsRecFile->Close();
293     delete fConsRecFile;
294     fConsRecFile = 0;
295   }
296   //
297 }
298
299 //_____________________________________________________________________________________________
300 Bool_t AliMillePede2::ReadNextRecordData()
301 {
302   // read next data record (if any)
303   if (!fTreeData || ++fCurrRecDataID >= fTreeData->GetEntries()) { fCurrRecDataID--; return kFALSE;}
304   fTreeData->GetEntry(fCurrRecDataID);
305   return kTRUE;
306 }
307
308 //_____________________________________________________________________________________________
309 Bool_t AliMillePede2::ReadNextRecordConstraint()
310 {
311   // read next constraint record (if any)
312   if (!fTreeConstr || ++fCurrRecConstrID >= fTreeConstr->GetEntries()) { fCurrRecConstrID--; return kFALSE;}
313   fTreeConstr->GetEntry(fCurrRecConstrID);
314   return kTRUE;
315 }
316
317 //_____________________________________________________________________________________________
318 void AliMillePede2::SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma)
319 {
320   if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data
321   //
322   // write data of single measurement
323   if (lSigma<=0.0) { // If parameter is fixed, then no equation
324     for (int i=fNLocPar; i--;) derlc[i] = 0.0;
325     for (int i=fNGloPar; i--;) dergb[i] = 0.0;
326     return;
327   }
328   //
329   fRecord->AddResidual(lMeas);
330   //
331   // Retrieve local param interesting indices
332   for (int i=0;i<fNLocPar;i++) if (derlc[i]!=0.0) {fRecord->AddIndexValue(i,derlc[i]); derlc[i] = 0.0;}
333   //
334   fRecord->AddWeight( 1.0/lSigma/lSigma );
335   //
336   // Idem for global parameters
337   for (int i=0;i<fNGloPar;i++) if (dergb[i]!=0.0) {
338     fRecord->AddIndexValue(i,dergb[i]); dergb[i] = 0.0;
339     fRecord->MarkGroup(fParamGrID[i]);
340   }
341   //
342 }
343
344 //_____________________________________________________________________________________________
345 void AliMillePede2::SetLocalEquation(int *indgb, double *dergb, int ngb, int *indlc,
346                                      double *derlc,int nlc,double lMeas,double lSigma)
347 {       
348   // write data of single measurement
349   if (lSigma<=0.0) { // If parameter is fixed, then no equation
350     for (int i=nlc;i--;)  derlc[i] = 0.0;
351     for (int i=ngb;i--;)  dergb[i] = 0.0;
352     return;
353   }
354   //
355   if (!fDataRecFile || !fDataRecFile->IsWritable()) InitDataRecStorage(); // create a buffer to store the data
356   //
357   fRecord->AddResidual(lMeas);
358   //
359   // Retrieve local param interesting indices
360   for (int i=0;i<nlc;i++) if (derlc[i]!=0.0) {fRecord->AddIndexValue(indlc[i],derlc[i]); derlc[i]=0.; indlc[i]=0;}
361   //
362   fRecord->AddWeight( 1./lSigma/lSigma );
363   //
364   // Idem for global parameters
365   for (int i=0;i<ngb;i++) if (dergb[i]!=0.0) {fRecord->AddIndexValue(indgb[i],dergb[i]); dergb[i]=0.; indgb[i]=0;} 
366   //
367 }
368
369
370 //_____________________________________________________________________________________________
371 void AliMillePede2::SetGlobalConstraint(double *dergb, double val, double sigma)
372 {       
373   // Define a constraint equation.
374   if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data
375   //
376   fRecord->Reset();
377   fRecord->AddResidual(val);
378   fRecord->AddWeight(sigma);
379   for (int i=0; i<fNGloPar; i++) if (dergb[i]!=0)  fRecord->AddIndexValue(i,dergb[i]);
380   fNGloConstraints++;
381   if (sigma==0) fNLagrangeConstraints++;
382   //  printf("NewConstraint:\n"); fRecord->Print(); //RRR
383   SaveRecordConstraint();
384   //
385 }
386
387 //_____________________________________________________________________________________________
388 void AliMillePede2::SetGlobalConstraint(const int *indgb, double *dergb, int ngb, double val,double sigma)
389 {       
390   // Define a constraint equation.
391   if (!fConsRecFile || !fConsRecFile->IsWritable()) InitConsRecStorage(); // create a buffer to store the data
392   fRecord->Reset();
393   fRecord->AddResidual(val);
394   fRecord->AddWeight(sigma);   // dummy
395   for (int i=0; i<ngb; i++) if (dergb[i]!=0)  fRecord->AddIndexValue(indgb[i],dergb[i]);
396   fNGloConstraints++;
397   if (sigma==0) fNLagrangeConstraints++;
398   SaveRecordConstraint();
399   //
400 }
401
402 //_____________________________________________________________________________________________
403 Int_t AliMillePede2::LocalFit(double *localParams)
404 {
405   /*
406     Perform local parameters fit once all the local equations have been set
407     -----------------------------------------------------------
408     localParams = (if !=0) will contain the fitted track parameters and
409     related errors
410   */
411   static int     nrefSize = 0;
412   //  static TArrayI refLoc,refGlo,nrefLoc,nrefGlo;
413   static Int_t  *refLoc=0,*refGlo=0,*nrefLoc=0,*nrefGlo=0;
414   int nPoints = 0;
415   //
416   AliSymMatrix    &matCLoc    = *fMatCLoc;
417   AliMatrixSq     &matCGlo    = *fMatCGlo;
418   AliRectMatrix   &matCGloLoc = *fMatCGloLoc;
419   //
420   memset(fVecBLoc,0,fNLocPar*sizeof(double));
421   matCLoc.Reset();      
422   //
423   int cnt=0;
424   int recSz = fRecord->GetSize();
425   //
426   while(cnt<recSz) {  // Transfer the measurement records to matrices
427     //
428     // extract addresses of residual, weight and pointers on local and global derivatives for each point
429     if (nrefSize<=nPoints) {
430       int *tmpA = 0;
431       nrefSize = 2*(nPoints+1); 
432       tmpA = refLoc;  refLoc  = new Int_t[nrefSize]; if (tmpA) memcpy(refLoc,tmpA,nPoints*sizeof(int));
433       tmpA = refGlo;  refGlo  = new Int_t[nrefSize]; if (tmpA) memcpy(refGlo,tmpA,nPoints*sizeof(int));
434       tmpA = nrefLoc; nrefLoc = new Int_t[nrefSize]; if (tmpA) memcpy(nrefLoc,tmpA,nPoints*sizeof(int));
435       tmpA = nrefGlo; nrefGlo = new Int_t[nrefSize]; if (tmpA) memcpy(nrefGlo,tmpA,nPoints*sizeof(int));
436     }
437     //
438     refLoc[nPoints]  = ++cnt;
439     int nLoc = 0;
440     while(!fRecord->IsWeight(cnt)) {nLoc++; cnt++;}
441     nrefLoc[nPoints] = nLoc;
442     //
443     refGlo[nPoints]  = ++cnt;
444     int nGlo = 0;
445     while(!fRecord->IsResidual(cnt) && cnt<recSz) {nGlo++; cnt++;} 
446     nrefGlo[nPoints] = nGlo;
447     //
448     nPoints++;
449   }
450   double vl;
451   //
452   for (int ip=nPoints;ip--;) {  // Transfer the measurement records to matrices
453     double  resid  = fRecord->GetValue( refLoc[ip]-1 );
454     double  weight = fRecord->GetValue( refGlo[ip]-1 );    
455     double *derLoc = fRecord->GetValue()+refLoc[ip];
456     double *derGlo = fRecord->GetValue()+refGlo[ip];
457     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
458     int    *indGlo = fRecord->GetIndex()+refGlo[ip];
459     //
460     for (int i=nrefGlo[ip];i--;) {       // suppress the global part (only relevant with iterations)
461       int iID = indGlo[i];              // Global param indice
462       if (fSigmaPar[iID]<=0.) continue;                                    // fixed parameter RRRCheck
463       if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]);  // linear parameter
464       else                resid -= derGlo[i]*fDeltaPar[iID];                  // nonlinear parameter
465     }
466     //
467     // Symmetric matrix, don't bother j>i coeffs
468     for (int i=nrefLoc[ip];i--;) {         // Fill local matrix and vector
469       fVecBLoc[ indLoc[i] ] += weight*resid*derLoc[i];
470       for (int j=i+1;j--;) matCLoc(indLoc[i] ,indLoc[j]) += weight*derLoc[i]*derLoc[j];
471     }   
472     //
473   } // end of the transfer of the measurement record to matrices
474   //
475   // first try to solve by faster Cholesky decomposition, then by Gaussian elimination
476   if (!matCLoc.SolveChol(fVecBLoc,kTRUE)) {
477     AliInfo("Failed to solve locals by Cholesky, traying Gaussian Elimination");
478     if (!matCLoc.SolveSpmInv(fVecBLoc,kTRUE)) { 
479       AliInfo("Failed to solve locals by Gaussian Elimination, skip...");
480       matCLoc.Print("d");
481       return 0; // failed to solve
482     }
483   }
484   //
485   // If requested, store the track params and errors
486   //  printf("locfit: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n");
487   if (localParams) for (int i=fNLocPar; i--;) {
488       localParams[2*i]   = fVecBLoc[i];
489       localParams[2*i+1] = TMath::Sqrt(TMath::Abs(matCLoc.QueryDiag(i)));
490     }
491   //
492   float lChi2 = 0;
493   int   nEq   = 0;
494   //
495   for (int ip=nPoints;ip--;) {  // Calculate residuals
496     double  resid  = fRecord->GetValue( refLoc[ip]-1 );
497     double  weight = fRecord->GetValue( refGlo[ip]-1 );    
498     double *derLoc = fRecord->GetValue()+refLoc[ip];
499     double *derGlo = fRecord->GetValue()+refGlo[ip];
500     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
501     int    *indGlo = fRecord->GetIndex()+refGlo[ip];
502     //
503     // Suppress local and global contribution in residuals;
504     for (int i=nrefLoc[ip];i--;) resid -= derLoc[i]*fVecBLoc[ indLoc[i] ];     // local part 
505     //
506     for (int i=nrefGlo[ip];i--;) {                                             // global part
507       int iID = indGlo[i];
508       if ( fSigmaPar[iID] <= 0.) continue;                                     // fixed parameter RRRCheck
509       if (fIsLinear[iID]) resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]);   // linear parameter
510       else                resid -= derGlo[i]*fDeltaPar[iID];                   // nonlinear parameter
511     }
512     //
513     // reject the track if the residual is too large (outlier)
514     double absres = TMath::Abs(resid);
515     if ( (absres >= fResCutInit && fIter ==1 ) ||
516          (absres >= fResCut     && fIter > 1)) {
517       if (fLocFitAdd) fNLocFitsRejected++;      
518       return 0;
519     }
520     // 
521     lChi2 += weight*resid*resid ; // total chi^2
522     nEq++;                        // number of equations                        
523   } // end of Calculate residuals
524   //
525   int nDoF = nEq-fNLocPar;
526   lChi2 = (nDoF>0) ? lChi2/nDoF : 0;  // Chi^2/dof  
527   //
528   if (fNStdDev != 0 && nDoF>0 && lChi2 > Chi2DoFLim(fNStdDev,nDoF)*fChi2CutFactor) { // check final chi2
529     if (fLocFitAdd) fNLocFitsRejected++;      
530     return 0;
531   }
532   //
533   if (fLocFitAdd) {
534     fNLocFits++;
535     fNLocEquations += nEq;
536   }
537   else {
538     fNLocFits--;
539     fNLocEquations -= nEq;
540   }
541   //
542   //  local operations are finished, track is accepted 
543   //  We now update the global parameters (other matrices)
544   //
545   int nGloInFit = 0;
546   //
547   for (int ip=nPoints;ip--;) {  // Update matrices
548     double  resid  = fRecord->GetValue( refLoc[ip]-1 );
549     double  weight = fRecord->GetValue( refGlo[ip]-1 );    
550     double *derLoc = fRecord->GetValue()+refLoc[ip];
551     double *derGlo = fRecord->GetValue()+refGlo[ip];
552     int    *indLoc = fRecord->GetIndex()+refLoc[ip];
553     int    *indGlo = fRecord->GetIndex()+refGlo[ip];
554     //
555     for (int i=nrefGlo[ip];i--;) {    // suppress the global part
556       int iID = indGlo[i];           // Global param indice
557       if ( fSigmaPar[iID] <= 0.) continue;                                         // fixed parameter RRRCheck
558       if (fIsLinear[iID])  resid -= derGlo[i]*(fInitPar[iID]+fDeltaPar[iID]);      // linear parameter
559       else                 resid -= derGlo[i]*fDeltaPar[iID];                      // nonlinear parameter
560     }
561     //
562     for (int ig=nrefGlo[ip];ig--;) {
563       int iIDg = indGlo[ig];   // Global param indice (the matrix line)          
564       if ( fSigmaPar[iIDg] <= 0.) continue;                                    // fixed parameter RRRCheck
565       if (fLocFitAdd) fVecBGlo[ iIDg ] += weight*resid*derGlo[ig]; //!!!
566       else            fVecBGlo[ iIDg ] -= weight*resid*derGlo[ig]; //!!!
567       //
568       // First of all, the global/global terms (exactly like local matrix)
569       int nfill = 0;
570       for (int jg=ig+1;jg--;) { // matCGlo is symmetric by construction  
571         int jIDg = indGlo[jg];
572         if ( fSigmaPar[jIDg] <= 0.) continue;                                    // fixed parameter RRRCheck
573         if ( (vl = weight*derGlo[ig]*derGlo[jg])!=0 ) {
574           fFillIndex[nfill]   = jIDg;
575           fFillValue[nfill++] = fLocFitAdd ? vl:-vl;
576         }
577       }
578       if (nfill) matCGlo.AddToRow(iIDg,fFillValue,fFillIndex,nfill);
579       //
580       // Now we have also rectangular matrices containing global/local terms.
581       int iCIDg = fGlo2CGlo[iIDg];  // compressed Index of index          
582       if (iCIDg == -1) {
583         Double_t *rowGL = matCGloLoc(nGloInFit);
584         for (int k=fNLocPar;k--;) rowGL[k] = 0.0;  // reset the row
585         iCIDg = fGlo2CGlo[iIDg] = nGloInFit;
586         fCGlo2Glo[nGloInFit++] = iIDg;
587       }
588       //
589       Double_t *rowGLIDg = matCGloLoc(iCIDg);
590       for (int il=nrefLoc[ip];il--;) rowGLIDg[ indLoc[il] ] += weight*derGlo[ig]*derLoc[il];
591       fProcPnt[iIDg] += fLocFitAdd ? 1:-1;       // update counter
592       //
593     }
594   } // end of Update matrices
595   //
596   // calculate fMatCGlo -= fMatCGloLoc * fMatCLoc * fMatCGloLoc^T
597   // and       fVecBGlo -= fMatCGloLoc * fVecBLoc
598   //
599   //-------------------------------------------------------------- >>>
600   double vll;
601   for (int iCIDg=0; iCIDg<nGloInFit; iCIDg++) {
602     int iIDg = fCGlo2Glo[iCIDg];
603     //
604     vl = 0;
605     Double_t *rowGLIDg =  matCGloLoc(iCIDg);
606     for (int kl=0;kl<fNLocPar;kl++) if (rowGLIDg[kl]) vl += rowGLIDg[kl]*fVecBLoc[kl];
607     if  (vl!=0) fVecBGlo[iIDg] -= fLocFitAdd ? vl : -vl;
608     //
609     int nfill = 0;
610     for (int jCIDg=0;jCIDg<=iCIDg; jCIDg++) {  
611       int jIDg = fCGlo2Glo[jCIDg];
612       //
613       vl = 0;
614       Double_t *rowGLJDg =  matCGloLoc(jCIDg);
615       for (int kl=0;kl<fNLocPar;kl++) {
616         // diag terms
617         if ( (vll=rowGLIDg[kl]*rowGLJDg[kl])!=0 ) vl += matCLoc.QueryDiag(kl)*vll;
618         //
619         // off-diag terms
620         for (int ll=0;ll<kl;ll++) {
621           if ( (vll=rowGLIDg[kl]*rowGLJDg[ll])!=0 ) vl += matCLoc(kl,ll)*vll;
622           if ( (vll=rowGLIDg[ll]*rowGLJDg[kl])!=0 ) vl += matCLoc(kl,ll)*vll;
623         }
624       }
625       if (vl!=0) {
626         fFillIndex[nfill]   = jIDg;
627         fFillValue[nfill++] = fLocFitAdd ? -vl : vl;
628       }
629     }
630     if (nfill)  matCGlo.AddToRow(iIDg,fFillValue,fFillIndex,nfill);
631   }
632   //
633   // reset compressed index array
634   //
635   for (int i=nGloInFit;i--;) { 
636     fGlo2CGlo[ fCGlo2Glo[i] ] = -1;
637     fCGlo2Glo[i] = -1;
638   }
639   //
640   //---------------------------------------------------- <<<
641   return 1;
642 }
643
644 //_____________________________________________________________________________________________
645 Int_t AliMillePede2::GlobalFit(Double_t *par, Double_t *error, Double_t *pull)
646 {
647   // performs a requested number of global iterations
648   fIter = 1;
649   //
650   TStopwatch sw; sw.Start();
651   //
652   int res = 0;
653   AliInfo("Starting Global fit.");
654   while (fIter<=fMaxIter) {
655     //
656     res = GlobalFitIteration();
657     if (!res) break;
658     //
659     if (fChi2CutFactor != fChi2CutRef) {    
660       fChi2CutFactor = TMath::Sqrt(fChi2CutFactor);
661       if (fChi2CutFactor < 1.2*fChi2CutRef) {
662         fChi2CutFactor = fChi2CutRef;
663         fIter = fMaxIter - 1;     // Last iteration
664       }
665     }
666     fIter++;
667   }
668   //
669   sw.Stop();
670   AliInfo(Form("Global fit %s, CPU time: %.1f",res ? "Converged":"Failed",sw.CpuTime()));  
671   if (!res) return 0;
672   //
673   if (par) for (int i=fNGloPar;i--;) par[i] = fInitPar[i]+fDeltaPar[i]; 
674   //
675   if (fGloSolveStatus==gkInvert) { // errors on params are available
676     if (error) for (int i=fNGloPar;i--;) error[i] = fProcPnt[i]>0 ? TMath::Sqrt(TMath::Abs(fMatCGlo->QueryDiag(i))) : 0.;
677     if (pull)  for (int i=fNGloPar;i--;) pull[i]  = fProcPnt[i]>0 && (fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QueryDiag(i))>0. && fSigmaPar[i]>0 
678                                            ? fDeltaPar[i]/TMath::Sqrt(fSigmaPar[i]*fSigmaPar[i]-fMatCGlo->QueryDiag(i)) : 0;
679   }
680   //
681   return 1;
682 }
683
684 //_____________________________________________________________________________________________
685 Int_t AliMillePede2::GlobalFitIteration()
686 {
687   // perform global parameters fit once all the local equations have been fitted
688   //
689   AliInfo(Form("Global Fit Iteration#%2d (Local Fit Chi^2 cut factor: %.2f)",fIter,fChi2CutFactor));
690   //
691   if (!fNGloPar || !fTreeData) {
692     AliInfo("No data was stored, aborting iteration");
693     return 0;
694   }
695   TStopwatch sw,sws; 
696   sw.Start();
697   sws.Stop();
698   //
699   if (!fConstrUsed) {
700     fConstrUsed = new Bool_t[fNGloConstraints];
701     memset(fConstrUsed,0,fNGloConstraints*sizeof(Bool_t));
702   }
703   // Reset all info specific for this step
704   AliMatrixSq& matCGlo = *fMatCGlo;
705   matCGlo.Reset();
706   memset(fProcPnt,0,fNGloPar*sizeof(Int_t));
707   //
708   fNGloConstraints = fTreeConstr ? fTreeConstr->GetEntries() : 0;
709   //
710   // count number of Lagrange constraints: they need new row/cols to be added
711   fNLagrangeConstraints = 0;
712   for (int i=0; i<fNGloConstraints; i++) {
713     ReadRecordConstraint(i);
714     if ( fRecord->GetValue(1)==0 ) fNLagrangeConstraints++; // exact constraint (no error) -> Lagrange multiplier 
715   }
716   //
717   // if needed, readjust the size of the global vector (for matrices this is done automatically)
718   if (!fVecBGlo || fNGloSize!=fNGloPar+fNLagrangeConstraints) {
719     delete[] fVecBGlo;   // in case some constraint was added between the two manual iterations
720     fNGloSize = fNGloPar+fNLagrangeConstraints;
721     fVecBGlo = new Double_t[fNGloSize];
722   }
723   memset(fVecBGlo,0,fNGloSize*sizeof(double));
724   //
725   fNLocFits         = 0;
726   fNLocFitsRejected = 0;
727   fNLocEquations    = 0;
728   //
729   //  Process data records and build the matrices
730   Long_t ndr = fTreeData->GetEntries();
731   AliInfo(Form("Building the Global matrix from %d data records",ndr));
732   if (!ndr) return 0;
733   //
734   TStopwatch swt; swt.Start();
735   fLocFitAdd = kTRUE;  // add contributions of matching tracks
736   for (Long_t i=0;i<ndr;i++) {
737     ReadRecordData(i);
738     LocalFit();
739     if ( (i%int(0.1*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr);
740   }
741   swt.Stop();
742   printf("%ld local fits done: ", ndr);
743   swt.Print();
744   sw.Start(kFALSE);
745   //
746   //
747   // ---------------------- Reject parameters with low statistics ------------>>
748   fNGloFix = 0;
749   if (fMinPntValid>1 && fNGroupsSet) {
750     //
751     printf("Checking parameters with statistics < %d\n",fMinPntValid);
752     TStopwatch swsup;
753     swsup.Start();
754     // 1) build the list of parameters to fix
755     Int_t fixArrSize = 10;
756     Int_t nFixedGroups = 0;
757     TArrayI fixGroups(fixArrSize);
758     //
759     for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones
760       int grID = fParamGrID[i];
761       if (grID<0) continue; // not in the group
762       if (fProcPnt[i]>=fMinPntValid) continue;
763       fProcPnt[i] = 0;
764       // check if the group is already accounted
765       Bool_t fnd = kFALSE;
766       for (int j=nFixedGroups;j--;) if (fixGroups[j]==grID) {fnd=kTRUE; break;}
767       if (fnd) continue;
768       if (nFixedGroups>=fixArrSize) {fixArrSize*=2; fixGroups.Set(fixArrSize);}
769       fixGroups[nFixedGroups++] = grID; // add group to fix
770       //
771     }
772     // 2) loop over records and add contributions of fixed groups with negative sign
773     fLocFitAdd = kFALSE;
774     //
775     for (Long_t i=0;i<ndr;i++) {
776       ReadRecordData(i);
777       Bool_t suppr = kFALSE;
778       for (int ifx=nFixedGroups;ifx--;)if (fRecord->IsGroupPresent(fixGroups[ifx])) suppr = kTRUE;
779       if (suppr) LocalFit();
780     }
781     fLocFitAdd = kTRUE;
782     //
783     if (nFixedGroups) {
784       printf("Suppressed contributions of groups with NPoints<%d :\n",fMinPntValid);
785       for (int i=0;i<nFixedGroups;i++) printf("%d ",fixGroups[i]); printf("\n");
786     }
787     swsup.Stop();
788     swsup.Print();
789   }
790   // ---------------------- Reject parameters with low statistics ------------<<
791   //
792   // add large number to diagonal of fixed params  
793   //
794   for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones
795     if (fProcPnt[i]<1) {
796       fNGloFix++; 
797       fVecBGlo[i] = 0.;
798       matCGlo.DiagElem(i) = float(fNLocEquations*fNLocEquations);
799     }
800     else matCGlo.DiagElem(i) += (fgWeightSigma ? fProcPnt[i] : 1.)/(fSigmaPar[i]*fSigmaPar[i]);
801   }
802   //
803   for (int i=fNGloPar;i--;) fDiagCGlo[i] = matCGlo.QueryDiag(i); // save the diagonal elements  
804   //
805   // add constraint equations
806   int nVar = fNGloPar;                    // Current size of global matrix      
807   for (int i=0; i<fNGloConstraints; i++) {
808     ReadRecordConstraint(i);
809     double val   = fRecord->GetValue(0);  
810     double sig   = fRecord->GetValue(1);  
811     int    *indV = fRecord->GetIndex()+2;
812     double *der  = fRecord->GetValue()+2;
813     int    csize = fRecord->GetSize()-2;
814     //
815     // check if after suppression of fixed variables there are non-0 derivatives
816     // and determine the max statistics of involved params
817     int nSuppressed = 0;
818     int maxStat = 1;
819     for (int j=csize;j--;) {
820       if (fProcPnt[indV[j]]<1) nSuppressed++; 
821       else {
822         maxStat = TMath::Max(maxStat,fProcPnt[indV[j]]);
823       }
824     }
825     //
826     if (nSuppressed==csize) {
827       //      AliInfo(Form("Neglecting constraint %d of %d derivatives since no free parameters left",i,csize));
828       //
829       // was this constraint ever created ? 
830       if ( sig==0 && fConstrUsed[i] ) { // this is needed only for constraints with Lagrange multiplier
831         // to avoid empty row impose dummy constraint on "Lagrange multiplier"
832         matCGlo.DiagElem(nVar) = 1.;
833         fVecBGlo[nVar++] = 0;
834       }
835       continue;
836     }
837     //
838     // account for already accumulated corrections
839     for (int j=csize; j--;) val -= der[j]*(fInitPar[ indV[j] ]+fDeltaPar[ indV[j] ]);
840     //
841     if (sig > 0) {  // this is a gaussian constriant: no Lagrange multipliers are added
842       //
843       double sig2i = (fgWeightSigma ? TMath::Sqrt(maxStat) : 1.)/sig/sig;
844       for (int ir=0;ir<csize;ir++) {
845         int iID = indV[ir];
846         for (int ic=0;ic<=ir;ic++) { // matrix is symmetric
847           int jID = indV[ic];
848           double vl = der[ir]*der[ic]*sig2i;
849           if (vl!=0) matCGlo(iID,jID) += vl;
850         }
851         fVecBGlo[iID] += val*der[ir]*sig2i;
852       }
853     }
854     else {   // this is exact constriant:  Lagrange multipliers must be added
855       for (int j=csize; j--;) {
856         int jID = indV[j];
857         if (fProcPnt[jID]<1) continue;                      // this parameter was fixed, don't put it into constraint
858         matCGlo(nVar,jID) = float(fNLocEquations)*der[j];   // fMatCGlo is symmetric, only lower triangle is filled  
859       }
860       //
861       if (matCGlo.QueryDiag(nVar)) matCGlo.DiagElem(nVar) = 0.0;
862       fVecBGlo[nVar++] = float(fNLocEquations)*val; //RS ? should we use here fNLocFits ? 
863       fConstrUsed[i] = kTRUE;
864     }
865   }
866   //
867   AliInfo(Form("Obtained %-7ld equations from %-7ld records (%-7ld rejected). Fixed %-4d globals",
868                fNLocEquations,fNLocFits,fNLocFitsRejected,fNGloFix));
869
870   //
871   sws.Start();
872   fGloSolveStatus = SolveGlobalMatEq();                     // obtain solution for this step
873   sws.Stop();
874   printf("Solve %d |",fIter); sws.Print();
875   //
876   sw.Stop();
877   AliInfo(Form("Iteration#%2d %s. CPU time: %.1f",fIter,fGloSolveStatus==gkFailed ? "Failed":"Converged",sw.CpuTime()));
878   if (fGloSolveStatus==gkFailed) return 0;
879   //
880   for (int i=fNGloPar;i--;) fDeltaPar[i] += fVecBGlo[i];    // Update global parameters values (for iterations)
881   //
882   //  PrintGlobalParameters();
883   return 1;
884 }
885
886 //_____________________________________________________________________________________________
887 Int_t AliMillePede2::SolveGlobalMatEq()
888 {
889   //
890   // solve global matrix equation MatCGlob*X=VecBGlo and store the result in the VecBGlo
891   //
892   /*
893   printf("GlobalMatrix\n");
894   fMatCGlo->Print();
895   printf("RHS\n");
896   for (int i=0;i<fNGloPar;i++) printf("%d %+e\n",i,fVecBGlo[i]);
897   */
898   //
899   if (!fgIsMatGloSparse) {
900     //
901     if (fNLagrangeConstraints==0) { // pos-def systems are faster to solve by Cholesky
902       if ( ((AliSymMatrix*)fMatCGlo)->SolveChol(fVecBGlo, fgInvChol) ) return fgInvChol ? gkInvert:gkNoInversion;
903       else AliInfo("Solution of Global Dense System by Cholesky failed, trying Gaussian Elimiation");
904     }
905     //
906     if (((AliSymMatrix*)fMatCGlo)->SolveSpmInv(fVecBGlo, kTRUE)) return gkInvert;
907     else AliInfo("Solution of Global Dense System by Gaussian Elimiation failed, trying iterative methods");
908   }
909   // try to solve by minres
910   TVectorD sol(fNGloSize);
911   //
912   AliMinResSolve *slv = new AliMinResSolve(fMatCGlo,fVecBGlo);
913   if (!slv) return gkFailed;
914   //
915   Bool_t res = kFALSE;
916   if      (fgIterSol == AliMinResSolve::kSolMinRes) 
917     res =  slv->SolveMinRes(sol,fgMinResCondType,fgMinResMaxIter,fgMinResTol);
918   else if (fgIterSol == AliMinResSolve::kSolFGMRes) 
919     res =  slv->SolveFGMRES(sol,fgMinResCondType,fgMinResMaxIter,fgMinResTol,fgNKrylovV);
920   else 
921     AliInfo(Form("Undefined Iteritive Solver ID=%d, only %d are defined",fgIterSol,AliMinResSolve::kNSolvers));
922   //
923   if (!res) return gkFailed;
924   for (int i=fNGloSize;i--;) fVecBGlo[i] = sol[i];
925   return gkNoInversion;
926   //
927 }
928
929 //_____________________________________________________________________________________________
930 Float_t AliMillePede2::Chi2DoFLim(int nSig, int nDoF) const
931 {
932   /// return the limit in chi^2/nd for n sigmas stdev authorized
933   // Only n=1, 2, and 3 are expected in input
934   int lNSig;
935   float sn[3]        =  {0.47523, 1.690140, 2.782170};
936   float table[3][30] = {{1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630,
937                          1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321,
938                          1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136,
939                          1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040},
940                         {4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736,
941                          1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658,
942                          1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346,
943                          1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742},
944                         {9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468,
945                          2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635,
946                          2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901,
947                          1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681}};
948   
949   if (nDoF < 1) {
950     return 0.0;
951   }
952   else {  
953     lNSig = TMath::Max(1,TMath::Min(nSig,3));
954     
955     if (nDoF <= 30) {    
956       return table[lNSig-1][nDoF-1];
957     }
958     else { // approximation
959       return ((sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3)))*
960               (sn[lNSig-1]+TMath::Sqrt(float(2*nDoF-3))))/float(2*nDoF-2);
961     }
962   }
963 }
964
965 //_____________________________________________________________________________________________
966 Int_t AliMillePede2::SetIterations(double lChi2CutFac)
967 {
968   // Number of iterations is calculated from lChi2CutFac 
969   fChi2CutFactor = TMath::Max(1.0, lChi2CutFac);
970   //
971   AliInfo(Form("Initial cut factor is %f",fChi2CutFactor));
972   fIter = 1; // Initializes the iteration process
973   return 1;
974 }
975
976 //_____________________________________________________________________________________________
977 Double_t AliMillePede2::GetParError(int iPar) const
978 {
979   // return error for parameter iPar
980   if (fGloSolveStatus==gkInvert) {
981     double res = fMatCGlo->QueryDiag(iPar);
982     if (res>=0) return TMath::Sqrt(res);
983   } 
984   return -1.;
985 }
986
987
988 //_____________________________________________________________________________________________
989 Int_t AliMillePede2::PrintGlobalParameters() const
990 {
991   ///  Print the final results into the logfile
992   double lError = 0.;
993   double lGlobalCor =0.;
994         
995   AliInfo("");
996   AliInfo("   Result of fit for global parameters");
997   AliInfo("   ===================================");
998   AliInfo("    I       initial       final       differ        lastcor        error       gcor       Npnt");
999   AliInfo("----------------------------------------------------------------------------------------------");
1000   //
1001   for (int i=0; i<fNGloPar; i++) {
1002     lError = GetParError(i);
1003     lGlobalCor = 0.0;
1004     //          
1005     double dg;
1006     if (fGloSolveStatus==gkInvert && TMath::Abs( (dg=fMatCGlo->QueryDiag(i)) *fDiagCGlo[i]) > 0) {    
1007       lGlobalCor = TMath::Sqrt(TMath::Abs(1.0-1.0/(dg*fDiagCGlo[i])));
1008       AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %6d",
1009                    i,fInitPar[i],fInitPar[i]+fDeltaPar[i],fDeltaPar[i],fVecBGlo[i],lError,lGlobalCor,fProcPnt[i]));
1010     }
1011     else {    
1012       AliInfo(Form("%d\t %.6f\t %.6f\t %.6f\t %.6f\t OFF\t OFF\t %6d",i,fInitPar[i],fInitPar[i]+fDeltaPar[i],
1013                    fDeltaPar[i],fVecBGlo[i],fProcPnt[i]));
1014     }
1015   }
1016   return 1;
1017 }