84fe5fcd77eb179f0983cf5462b6fe6a2731b88d
[u/mrichter/AliRoot.git] / ITS / AliITSAlignMille2.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //-----------------------------------------------------------------------------
19 //
20 //  Interface to AliMillePede2 alignment class for the ALICE ITS detector
21 // 
22 //  ITS specific alignment class which interface to AliMillepede.   
23 //  For each track ProcessTrack calculates the local and global derivatives
24 //  at each hit and fill the corresponding local equations. Provide methods for
25 //  fixing or constraning detection elements for best results. 
26 // 
27 //  author M. Lunardon (thanks to J. Castillo), ruben.shahoyan@cern.ch
28 //-----------------------------------------------------------------------------
29
30 #include <TFile.h>
31 #include <TGrid.h>
32 #include <TClonesArray.h>
33 #include <TMath.h>
34 #include <TVirtualFitter.h>
35 #include <TGeoManager.h>
36 #include <TSystem.h>
37 #include <TRandom.h>
38 #include <TCollection.h>
39 #include <TGeoPhysicalNode.h>
40 #include <TMap.h>
41 #include <TObjString.h>
42 #include <TString.h>
43 #include "AliITSAlignMille2.h"
44 #include "AliITSgeomTGeo.h"
45 #include "AliGeomManager.h"
46 #include "AliMillePede2.h"
47 #include "AliTrackPointArray.h"
48 #include "AliAlignObjParams.h"
49 #include "AliLog.h"
50 #include "AliTrackFitterRieman.h"
51 #include "AliITSAlignMille2Constraint.h"
52 #include "AliITSAlignMille2ConstrArray.h"
53 #include "AliITSresponseSDD.h"
54 #include "AliITSTPArrayFit.h"
55 #include "AliCDBManager.h"
56 #include "AliCDBStorage.h"
57 #include "AliCDBEntry.h"
58 #include "AliITSsegmentationSDD.h"
59 #include "AliITSDriftSpeedArraySDD.h"
60 #include "AliITSCorrectSDDPoints.h"
61 #include "AliESDVertex.h"
62
63 ClassImp(AliITSAlignMille2)
64
65 const Char_t* AliITSAlignMille2::fgkRecKeys[] = {
66   "OCDB_PATH",
67   "OCDB_SPECIFIC",
68   "GEOMETRY_FILE",
69   "SUPERMODULE_FILE",
70   "CONSTRAINTS_REFERENCE_FILE",
71   "PREALIGNMENT_FILE",
72   "PRECALIBSDD_FILE",
73   "PREVDRIFTSDD_FILE",
74   "PRECORRMAPSDD_FILE",
75   "INITCORRMAPSDD_FILE",
76   "INITCALBSDD_FILE",
77   "INITVDRIFTSDD_FILE",
78   "INITDELTA_FILE",
79   "INITGEOM_FILE",
80   "SET_GLOBAL_DELTAS",
81   "CONSTRAINT_LOCAL",
82   "MODULE_VOLUID",
83   "MODULE_INDEX",
84   "SET_PSEUDO_PARENTS",
85   "SET_TRACK_FIT_METHOD",
86   "SET_MINPNT_TRA",
87   "SET_NSTDDEV",
88   "SET_RESCUT_INIT",
89   "SET_RESCUT_OTHER",
90   "SET_LOCALSIGMAFACTOR",
91   "SET_STARTFAC",
92   "SET_FINALFAC",
93   "SET_B_FIELD",
94   "SET_SPARSE_MATRIX",
95   "REQUIRE_POINT",
96   "CONSTRAINT_ORPHANS",
97   "CONSTRAINT_SUBUNITS",
98   "APPLY_CONSTRAINT",
99   "SET_EXTRA_CLUSTERS_MODE",
100   "SET_USE_TPAFITTER",
101   "SET_USE_LOCAL_YERROR",
102   "SET_MIN_POINTS_PER_MODULE",
103   "SET_USE_SDDVDCORRMULT",
104   "SET_WEIGHT_PT",
105   "SET_USE_DIAMOND",
106   "CORRECT_DIAMOND",
107   "SET_USE_VERTEX",
108   "SET_SAME_SDDT0"
109 };
110
111 const Char_t AliITSAlignMille2::fgkXYZ[] = "XYZ";
112
113 //========================================================================================================
114
115 AliITSAlignMille2* AliITSAlignMille2::fgInstance = 0;  
116 Int_t              AliITSAlignMille2::fgInstanceID = 0;
117
118 //________________________________________________________________________________________________________
119 AliITSAlignMille2::AliITSAlignMille2(const Char_t *configFilename,TList *userInfo ) 
120 : TObject(),
121   fMillepede(0),
122   fStartFac(16.), 
123   fFinalFac(1.), 
124   fResCutInitial(100.), 
125   fResCut(100.),
126   fNGlobal(0),
127   fNLocal(4),
128   fNStdDev(3),
129   fIsMilleInit(kFALSE),
130   fAllowPseudoParents(kFALSE),
131   //
132   fTPAFitter(0),
133   fCurrentModule(0),
134   fTrack(0),
135   fTrackBuff(0),
136   fCluster(),
137   fCurrentSensID(-1),
138   fClusLoc(12*3),
139   fClusGlo(12*3),
140   fClusSigLoc(12*3),
141   fGlobalDerivatives(0),
142   fMeasLoc(0),
143   fMeasGlo(0),
144   fSigmaLoc(0),
145   fConstrPT(-1),
146   fConstrPTErr(-1),
147   fConstrCharge(0),
148   fRunID(0),
149   //
150   fMinNPtsPerTrack(3),
151   fIniTrackParamsMeth(1),
152   fTotBadLocEqPoints(0),
153   fRieman(0),
154   //
155   fConstraints(0),
156   fCacheMatrixOrig(kMaxITSSensID+1),
157   fCacheMatrixCurr(kMaxITSSensID+1),
158   //
159   fUseGlobalDelta(kFALSE),
160   fTempExcludedModule(-1),
161   fUserProvided(0),
162   //
163   fIniUserInfo(userInfo),
164   fIniGeomPath(""),
165   fIniDeltaPath(""),
166   fIniSDDRespPath(""),
167   fPreCalSDDRespPath(""),
168   fIniSDDVDriftPath(""),
169   fPreSDDVDriftPath(""),
170   fIniSDDCorrMapPath(""),
171   fPreSDDCorrMapPath(""),
172   fConvertPreDeltas(kFALSE),
173   fGeometryPath(""),
174   fPreDeltaPath(""),
175   fConstrRefPath(""),
176   fDiamondPath(""),
177   fGeoManager(0),
178   fIsConfigured(kFALSE),
179   fPreAlignQF(0),
180 //
181   fIniRespSDD(0),
182   fPreRespSDD(0),
183   fIniVDriftSDD(0),
184   fPreVDriftSDD(0),
185   fIniCorrMapSDD(0),
186   fPreCorrMapSDD(0),
187   fSegmentationSDD(0),
188   fPrealignment(0),
189   fConstrRef(0),
190   fMilleModule(2),
191   fSuperModule(2),
192   fNModules(0),
193   fNSuperModules(0),
194   fUsePreAlignment(kFALSE),
195   fUseLocalYErr(kFALSE),
196   fBOn(kFALSE),
197   fBField(0.0),
198   fDataType(kCosmics),
199   fMinPntPerSens(0),
200   fBug(0),
201   fMilleVersion(2),
202   fExtraClustersMode(0),
203   fTrackWeight(1),
204   fWeightPt(0.),
205   fIsSDDVDriftMult(kFALSE),
206   fDiamond(),
207   fDiamondI(),
208   fUseDiamond(kFALSE),
209   fUseVertex(kFALSE),
210   fVertexSet(kFALSE),
211   fDiamondPointID(-1),
212   fDiamondModID(-1),
213   fCheckDiamondPoint(kDiamondCheckIfPrompt),
214   fFixCurvIfConstraned(kTRUE),
215   fCurvFitWasConstrained(kFALSE),
216   fConvAlgMatOld(100)
217 {
218   /// main constructor that takes input from configuration file
219   for (int i=3;i--;) fSigmaFactor[i] = 1.0;
220   //
221   // new RS
222   for (int i=0;i<3;i++) {
223     fCorrDiamond[i] = 0;
224   }
225   for (int itp=0;itp<kNDataType;itp++) {
226     fRequirePoints[itp] = kFALSE;
227     for (Int_t i=0; i<6; i++) {
228       fNReqLayUp[itp][i]=0;
229       fNReqLayDown[itp][i]=0;
230       fNReqLay[itp][i]=0;
231     }
232     for (Int_t i=0; i<3; i++) {
233       fNReqDetUp[itp][i]=0;
234       fNReqDetDown[itp][i]=0;
235       fNReqDet[itp][i]=0;
236     }
237   }
238   //
239   //  if (ProcessUserInfo(userInfo)) exit(1);
240   //
241   fDiamond.SetVolumeID(kVtxSensVID);
242   fDiamondI.SetVolumeID(kVtxSensVID);
243   float xyzd[3] = {0,0,0};
244   float covd[6] = {1,0,0,1,0,1e4}; 
245   fDiamond.SetXYZ(xyzd,covd); // dummy diamond
246   covd[5] = 1e-4;
247   fDiamondI.SetXYZ(xyzd,covd);
248   //
249   Int_t lc=LoadConfig(configFilename);
250   if (lc) {
251     AliError(Form("Error %d loading configuration from %s",lc,configFilename));
252     exit(1);
253   }
254   //
255   fMillepede = new AliMillePede2();  
256   fgInstance = this;
257   fgInstanceID++;
258   ResetCovIScale();
259   //
260 }
261
262 //________________________________________________________________________________________________________
263 AliITSAlignMille2::~AliITSAlignMille2()
264 {
265   /// Destructor
266   delete fMillepede;
267   delete[] fGlobalDerivatives;
268   delete fRieman;
269   delete fPrealignment;
270   delete fConstrRef;
271   delete fPreRespSDD;
272   delete fIniRespSDD;
273   delete fSegmentationSDD;
274   delete fIniVDriftSDD;
275   delete fPreVDriftSDD;
276   delete fIniCorrMapSDD;
277   delete fPreCorrMapSDD;
278   delete fTPAFitter;
279   fCacheMatrixOrig.Delete();
280   fCacheMatrixCurr.Delete();
281   fTrackBuff.Delete();
282   fConstraints.Delete();
283   fMilleModule.Delete();
284   fSuperModule.Delete();
285   if (--fgInstanceID==0) fgInstance = 0;
286 }
287
288 ///////////////////////////////////////////////////////////////////////
289
290 //________________________________________________________________________________________________________
291 TObjArray* AliITSAlignMille2::GetConfigRecord(FILE* stream, TString& recTitle, TString& recOpt, Bool_t rew)
292 {
293   // read new record from config file
294   TString record;
295   static TObjArray* recElems = 0;
296   if (recElems) {delete recElems; recElems = 0;}
297   recOpt = "";
298   //
299   TString keyws = recTitle;
300   if (!keyws.IsNull()) {
301     keyws.ToUpper();
302     //    keyws += " ";
303   }
304   while (record.Gets(stream)) {
305     int cmt=record.Index("#"); 
306     if (cmt>=0) record.Remove(cmt);  // skip comment
307     record.ReplaceAll("\t"," ");
308     record.ReplaceAll("\r"," ");
309     record.Remove(TString::kBoth,' '); 
310     if (record.IsNull()) continue;      // nothing to decode 
311     if (!keyws.IsNull() && !record.BeginsWith(keyws.Data())) continue; // specific record was requested
312     //
313     recElems = record.Tokenize(" ");
314     recTitle = recElems->At(0)->GetName();
315     recTitle.ToUpper();
316     recOpt = recElems->GetLast()>0 ? recElems->At(1)->GetName() : "";
317     break;
318   }
319   if (rew || !recElems) rewind(stream);
320   return recElems;
321 }
322
323 //________________________________________________________________________________________________________
324 Int_t AliITSAlignMille2::CheckConfigRecords(FILE* stream)
325 {  
326   TString record,recTitle;
327   int lineCnt = 0;
328   rewind(stream);
329   while (record.Gets(stream)) {
330     int cmt=record.Index("#"); 
331     lineCnt++;
332     if (cmt>=0) record.Remove(cmt);  // skip comment
333     record.ReplaceAll("\t"," ");
334     record.ReplaceAll("\r"," ");
335     record.Remove(TString::kBoth,' ');
336     if (record.IsNull()) continue;   // nothing to decode  
337     // extract keyword
338     int spc = record.Index(" ");
339     if (spc>0) recTitle = record(0,spc);
340     else     recTitle = record;
341     recTitle.ToUpper();
342     Bool_t strOK = kFALSE;
343     for (int ik=kNKeyWords;ik--;) if (recTitle == fgkRecKeys[ik]) {strOK = kTRUE; break;}
344     if (strOK) continue;
345     //
346     AliError(Form("Unknown keyword %s at line %d",
347                   recTitle.Data(),lineCnt));
348     return -1;
349     //
350   }
351   //
352   rewind(stream);
353   return 0;
354 }
355
356
357 //________________________________________________________________________________________________________
358 Int_t AliITSAlignMille2::LoadConfig(const Char_t *cfile)
359 {  
360   // return 0 if success
361   //        1 if error in module index or voluid
362   //
363   AliInfo(Form("Loading MillePede2 configuration from %s",cfile));
364   AliCDBManager::Instance()->SetCacheFlag(kFALSE);
365   FILE *pfc=fopen(cfile,"r");
366   if (!pfc) return -1;
367   //
368   TString record,recTitle,recOpt,recExt;
369   Int_t nrecElems,irec;
370   TObjArray *recArr=0;
371   //
372   fNModules = 0;
373   Bool_t stopped = kFALSE;
374   //
375   if (CheckConfigRecords(pfc)<0) return -1;
376   //
377   while(1) { 
378     //
379     // ============= 1: we read some important records in predefined order ================
380     //  
381     recTitle = fgkRecKeys[kOCDBDefaultPath];
382     if ( GetConfigRecord(pfc,recTitle,recOpt,1) && !recOpt.IsNull() ) {
383       AliInfo(Form("Configuration sets OCDB default storage to %s",recOpt.Data()));
384       AliCDBManager::Instance()->SetDefaultStorage( gSystem->ExpandPathName(recOpt.Data()) );
385       TObjString* objStr = (TObjString*)AliCDBManager::Instance()->GetStorageMap()->GetValue("default");
386       if (!objStr) {stopped = kTRUE; break;}
387       objStr->SetUniqueID(1); // mark as user set
388     }
389     //
390     if (fIniUserInfo && ProcessUserInfo(fIniUserInfo)) { AliError("Failed to process intial User Info"); stopped = kTRUE; break;}
391     //  
392     recTitle = fgkRecKeys[kGeomFile];
393     if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) fGeometryPath = gSystem->ExpandPathName(recOpt.Data()); 
394     if ( LoadGeometry(fGeometryPath) ) { AliError("Failed to find/load target ideal Geometry"); stopped = kTRUE; break;}
395     //
396     // Do we use new TrackPointArray fitter ?
397     recTitle = fgkRecKeys[kTPAFitter];
398     if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) fTPAFitter = new AliITSTPArrayFit(kNLocal);
399     //
400     recTitle = fgkRecKeys[kSuperModileFile];
401     if ( !GetConfigRecord(pfc,recTitle,recOpt,1) || 
402          recOpt.IsNull()                         || 
403          gSystem->ExpandPathName(recOpt)         ||
404          gSystem->AccessPathName(recOpt.Data())  ||
405          LoadSuperModuleFile(recOpt.Data()))
406       { AliError("Failed to find/load SuperModules"); stopped = kTRUE; break;}
407     //
408     recTitle = fgkRecKeys[kConstrRefFile];      // LOCAL_CONSTRAINTS are defined wrt these deltas
409     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
410       if (recOpt.IsNull() || recOpt=="IDEAL") SetConstraintWrtRef( "IDEAL" );
411       else {
412         for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
413         if ( SetConstraintWrtRef(recOpt.Data()) )
414           { AliError("Failed to load reference deltas for local constraints"); stopped = kTRUE; break;}
415       }
416     }
417     //
418     recTitle = fgkRecKeys[kInitGeomFile];
419     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1))  && !recOpt.IsNull() ) {
420       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
421       fIniGeomPath = recOpt;
422       gSystem->ExpandPathName(fIniGeomPath);
423       fUserProvided |= kSameInitGeomBit;
424       AliInfo(Form("Configuration sets Production Geometry to %s",fIniGeomPath.Data()));
425     }
426     if (fIniGeomPath.IsNull()) fIniGeomPath = fGeometryPath;
427     //   
428     recTitle = fgkRecKeys[kInitDeltaFile];
429     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1))  && !recOpt.IsNull() ) {
430       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
431       fIniDeltaPath = recOpt;
432       gSystem->ExpandPathName(fIniDeltaPath);
433       fUserProvided |= kSameInitDeltasBit;
434       AliInfo(Form("Configuration sets Production Deltas to %s",fIniDeltaPath.Data()));
435     }
436     //
437     recTitle = fgkRecKeys[kPreDeltaFile];
438     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
439       if (!recOpt.IsNull()) {
440         for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
441         fPreDeltaPath = recOpt;
442         gSystem->ExpandPathName(fPreDeltaPath);
443       }
444       else if (!fIniDeltaPath.IsNull()) {
445         AliInfo("PreAlignment Deltas keyword is present but empty, will set to Init Deltas of the first tree");
446         fPreDeltaPath = fIniDeltaPath;  
447         if (fIniGeomPath != fGeometryPath) fConvertPreDeltas = kTRUE; // production and target geometries are different, request conversion
448       }
449       AliInfo(Form("Configuration sets PreAlignment Deltas to %s",fPreDeltaPath.Data()));
450     }
451     //
452     // if initial deltas were provided, load them, apply to geometry and store are "original" matrices
453     if (CacheMatricesOrig()) {stopped = kTRUE; break;}
454     //   
455     // then load prealignment deltas
456     if (!fPreDeltaPath.IsNull()) {
457       if (fConvertPreDeltas) ConvertDeltas();   // Prealignment deltas are the same as production ones, but need conversion to new geometry
458       else if (LoadDeltas(fPreDeltaPath,fPrealignment)) {stopped = kTRUE; break;} // read deltas from the file
459     }
460     if (fPrealignment && ApplyToGeometry()) {stopped = kTRUE; break;}
461     //
462     recTitle = fgkRecKeys[ kInitCalSDDFile ];
463     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
464       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
465       fIniSDDRespPath = recOpt;
466       gSystem->ExpandPathName(fIniSDDRespPath);
467       fUserProvided |= kSameInitSDDRespBit;
468       AliInfo(Form("Configuration sets Production SDD Response to %s",fIniSDDRespPath.Data()));
469     }
470     if (LoadSDDResponse(fIniSDDRespPath, fIniRespSDD) ) {stopped = kTRUE; break;}
471     //
472     //
473     recTitle = fgkRecKeys[ kInitCorrMapSDDFile ];
474     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
475       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
476       fIniSDDCorrMapPath = recOpt;
477       gSystem->ExpandPathName(fIniSDDCorrMapPath);
478       fUserProvided |= kSameInitSDDCorrMapBit;
479       AliInfo(Form("Configuration sets Production SDD Correction Map to %s",fIniSDDCorrMapPath.Data()));
480     }
481     if (LoadSDDCorrMap(fIniSDDCorrMapPath, fIniCorrMapSDD) ) {stopped = kTRUE; break;}
482     //
483     recTitle = fgkRecKeys[kPreCalSDDFile];
484     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
485       if (!recOpt.IsNull()) {
486         for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
487         fPreCalSDDRespPath = recOpt;
488         gSystem->ExpandPathName(fPreCalSDDRespPath);
489       }
490       else if (!fIniSDDRespPath.IsNull()) {
491         AliInfo("PreCalibration SDD response keyword is present but empty, will set to Init SDD repsonse");
492         fPreCalSDDRespPath = fIniSDDRespPath;   
493       }
494       AliInfo(Form("Configuration sets PreCalibration SDD Response to %s",fPreCalSDDRespPath.Data()));
495     }
496     //
497     if (LoadSDDResponse(fPreCalSDDRespPath, fPreRespSDD) ) {stopped = kTRUE; break;}
498     //
499     recTitle = fgkRecKeys[kPreCorrMapSDDFile];
500     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
501       if (!recOpt.IsNull()) {
502         for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
503         fPreSDDCorrMapPath = recOpt;
504         gSystem->ExpandPathName(fPreSDDCorrMapPath);
505       }
506       else if (!fIniSDDCorrMapPath.IsNull()) {
507         AliInfo("PreCalibration SDD Correction Map keyword is present but empty, will set to Init SDD Correction Map");
508         fPreSDDCorrMapPath = fIniSDDCorrMapPath;
509       }
510       AliInfo(Form("Configuration sets PreCalibration SDD Correction Map to %s",fPreSDDCorrMapPath.Data()));
511     }
512     //
513     if (LoadSDDCorrMap(fPreSDDCorrMapPath, fPreCorrMapSDD) ) {stopped = kTRUE; break;}
514     //    //
515     recTitle = fgkRecKeys[ kInitVDriftSDDFile ];
516     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
517       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
518       fIniSDDVDriftPath = recOpt;
519       gSystem->ExpandPathName(fIniSDDVDriftPath);
520       fUserProvided |= kSameInitSDDVDriftBit;
521       AliInfo(Form("Configuration sets Production SDD VDrift to %s",fIniSDDVDriftPath.Data()));
522     }
523     if (LoadSDDVDrift(fIniSDDVDriftPath, fIniVDriftSDD) ) {stopped = kTRUE; break;}
524     //
525     recTitle = fgkRecKeys[ kPreVDriftSDDFile ];
526     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
527       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
528       fPreSDDVDriftPath = recOpt;
529       gSystem->ExpandPathName(fPreSDDVDriftPath);
530       AliInfo(Form("Configuration sets PreCalibration SDD VDrift to %s",fPreSDDVDriftPath.Data()));
531       if (LoadSDDVDrift(fPreSDDVDriftPath, fPreVDriftSDD) ) {stopped = kTRUE; break;}
532     }
533     //
534     recTitle = fgkRecKeys[ kGlobalDeltas ];
535     if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) SetUseGlobalDelta(kTRUE);
536     //
537     recTitle = fgkRecKeys[ kUseDiamond ];
538     if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
539       if (!GetUseGlobalDelta()) {
540         AliError("Diamond constraint is supported only for Global Frame mode");
541         stopped = kTRUE; 
542         break;
543       }
544       fUseDiamond = kTRUE;
545       if (!recOpt.IsNull()) {
546         fDiamondPath = recOpt;
547         gSystem->ExpandPathName(fDiamondPath);
548         fUserProvided |= kSameDiamondBit;
549         AliInfo(Form("Configuration sets Diamond constraint to %s",fDiamondPath.Data()));
550       }
551     }
552     //
553     recTitle = fgkRecKeys[ kUseVertex ];
554     if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
555       if (!GetUseGlobalDelta()) {
556         AliError("Vertex constraint is supported only for Global Frame mode");
557         stopped = kTRUE; 
558         break;
559       }
560       fUseVertex = kTRUE;
561       if (fUseDiamond) {
562         AliError("Cannot use Vertex and Diamond constraints at the same time");
563         stopped = kTRUE; 
564         break;
565       }
566       AliInfo("Will use Vertex constraint when available");
567     }
568     // =========== 2: see if there are local gaussian constraints defined =====================
569     //            Note that they should be loaded before the modules declaration
570     //
571     recTitle = fgkRecKeys[ kConstrLocal ];
572     while( (recArr=GetConfigRecord(pfc,recTitle,recOpt,0)) ) {
573       nrecElems = recArr->GetLast()+1;
574       if (recOpt.IsFloat()) {stopped = kTRUE; break;} // wrong name
575       if (GetConstraint(recOpt.Data())) {
576         AliError(Form("Existing constraint %s repeated",recOpt.Data()));
577         stopped = kTRUE; break;
578       }
579       recExt = recArr->At(2)->GetName();
580       if (!recExt.IsFloat()) {stopped = kTRUE; break;}
581       double val = recExt.Atof();      
582       recExt = recArr->At(3)->GetName();
583       if (!recExt.IsFloat()) {stopped = kTRUE; break;}
584       double err = recExt.Atof();      
585       int nwgh = nrecElems - 4;
586       double *wgh = new double[nwgh];
587       for (nwgh=0,irec=4;irec<nrecElems;irec++) {
588         recExt = recArr->At(irec)->GetName();
589         if (!recExt.IsFloat()) {stopped = kTRUE; break;}
590         wgh[nwgh++] = recExt.Atof();
591       }
592       if (stopped) {delete[] wgh; break;}
593       //
594       ConstrainLocal(recOpt.Data(),wgh,nwgh,val,err);
595       delete[] wgh;
596       //
597     } // end while for loop over local constraints
598     if (stopped) break;
599     //
600     // =========== 3: now read modules to align ===================================
601     //
602     rewind(pfc);
603     // create fixed modules
604     for (int j=0; j<fNSuperModules; j++) {
605       AliITSAlignMille2Module* proto = GetSuperModule(j);
606       if (!proto->IsAlignable()) continue;
607       AliITSAlignMille2Module* mod = new AliITSAlignMille2Module(*proto);
608       // the matrix might be updated in case some prealignment was applied, check 
609       TGeoHMatrix* mup = AliGeomManager::GetMatrix(mod->GetName());
610       if (mup) *(mod->GetMatrix()) = *mup;
611       fMilleModule.AddAtAndExpand(mod,fNModules);
612       mod->SetGeomParamsGlobal(fUseGlobalDelta);
613       mod->SetUniqueID(fNModules++);
614       mod->SetNotInConf(kTRUE);
615     }
616     CreateVertexModule();
617     //
618     while( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0)) ) {
619       if (!(recTitle==fgkRecKeys[ kModVolID ] || recTitle==fgkRecKeys[ kModIndex ])) continue;
620       // Expected format: MODULE id tolX tolY tolZ tolPsi tolTh tolPhi [[sigX sigY sigZ]  extra params]
621       // where tol* is the tolerance (sigma) for given DOF. 0 means fixed
622       // sig* is the scaling parameters for the errors of the clusters of this module
623       // extra params are defined for specific modules, e.g. t0 and vdrift corrections of SDD
624       //
625       nrecElems = recArr->GetLast()+1;
626       if (nrecElems<2 || !recOpt.IsDigit()) {stopped = kTRUE; break;}
627       int idx = recOpt.Atoi(); 
628       UShort_t voluid =  (idx<=kMaxITSSensID) ? GetModuleVolumeID(idx) : idx;
629       AliITSAlignMille2Module* mod = 0;
630       //
631       if (voluid>=kMinITSSupeModuleID) { // custom supermodule
632         mod = GetMilleModuleByVID(voluid);
633         if (!mod) { // need to create
634           for (int j=0; j<fNSuperModules; j++) {
635             if (voluid==GetSuperModule(j)->GetVolumeID()) {
636               mod = new AliITSAlignMille2Module(*GetSuperModule(j));
637               // the matrix might be updated in case some prealignment was applied, check 
638               TGeoHMatrix* mup = AliGeomManager::GetMatrix(mod->GetName());
639               if (mup) *(mod->GetMatrix()) = *mup;
640               fMilleModule.AddAtAndExpand(mod,fNModules);
641               mod->SetGeomParamsGlobal(fUseGlobalDelta);
642               mod->SetUniqueID(fNModules++);
643               break;
644             }   
645           }
646         }
647         mod->SetNotInConf(kFALSE);
648       }
649       else if (idx<=kMaxITSSensVID) {
650         mod = new AliITSAlignMille2Module(voluid);
651         fMilleModule.AddAtAndExpand(mod,fNModules);
652         mod->SetGeomParamsGlobal(fUseGlobalDelta);
653         mod->SetUniqueID(fNModules++);
654       }
655       if (!mod) {stopped = kTRUE; break;}  // bad volid
656       //
657       // geometry variation settings
658       for (int i=0;i<AliITSAlignMille2Module::kMaxParGeom;i++) {
659         irec = i+2;
660         if (irec >= nrecElems) break;
661         recExt = recArr->At(irec)->GetName();
662         if (!recExt.IsFloat()) {stopped = kTRUE; break;}
663         mod->SetFreeDOF(i, recExt.Atof() );     
664       }
665       if (stopped) break;
666       //
667       // scaling factors for cluster errors
668       // first set default ones
669       for (int i=0;i<3;i++) mod->SetSigmaFactor(i, fSigmaFactor[i]);    
670       for (int i=0;i<3;i++) {
671         irec = i+8;
672         if (irec >= nrecElems) break;
673         recExt = recArr->At(irec)->GetName();
674         if (!recExt.IsFloat()) {stopped = kTRUE; break;}
675         mod->SetSigmaFactor(i, recExt.Atof() ); 
676       }     
677       if (stopped) break;
678       //
679       // now comes special detectors treatment
680       if (mod->IsSDD()) {
681         double vl = 0;
682         if (nrecElems>11) {
683           recExt = recArr->At(11)->GetName();
684           if (recExt.IsFloat()) vl = recExt.Atof();
685           else {stopped = kTRUE; break;}
686           irec = 11;
687         }
688         mod->SetFreeDOF(AliITSAlignMille2Module::kDOFT0,vl);
689         //
690         Bool_t cstLR = kFALSE;
691         for (int lr=0;lr<2;lr++) { // left right side vdrift corrections
692           vl = 0;
693           if (nrecElems>12+lr) {
694             recExt = recArr->At(12+lr)->GetName();
695             if (recExt.IsFloat()) vl = recExt.Atof();
696             else {stopped = kTRUE; break;}
697             irec = 12+lr;
698           }
699           mod->SetFreeDOF(lr==0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR,vl);
700           if (lr==1 && vl>=10) cstLR = kTRUE;  // the right side should be constrained to left one 
701         }
702         if (cstLR) mod->SetVDriftLRSame();
703       }
704       //
705       mod->EvaluateDOF();
706       //
707       // now check if there are local constraints on this module
708       for (++irec;irec<nrecElems;irec++) {
709         recExt = recArr->At(irec)->GetName();
710         if (recExt.IsFloat()) {stopped=kTRUE;break;}
711         AliITSAlignMille2ConstrArray* cstr = (AliITSAlignMille2ConstrArray*)GetConstraint(recExt.Data());
712         if (!cstr) {
713           AliInfo(Form("No Local constraint %s was declared",recExt.Data())); 
714           stopped=kTRUE; 
715           break;
716         }
717         cstr->AddModule(mod);
718       }
719       if (stopped) break;
720     } // end while for loop over modules
721     if (stopped) break;
722     //
723     if (fNModules==0) {AliError("Failed to find any MODULE"); stopped = kTRUE; break;}  
724     BuildHierarchy();  // preprocess loaded modules
725     //
726     // =========== 4: the rest may come in arbitrary order =======================================
727     rewind(pfc);
728     while ( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0))!=0 ) {
729       //
730       nrecElems = recArr->GetLast()+1;
731       //
732       // some simple flags -----------------------------------------------------------------------
733       //
734       if      (recTitle == fgkRecKeys[ kPseudoParents ])  SetAllowPseudoParents(kTRUE);
735       //
736       // some optional parameters ----------------------------------------------------------------
737       else if (recTitle == fgkRecKeys[ kTrackFitMethod ]) {
738         if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
739         SetInitTrackParamsMeth(recOpt.Atoi());
740       }
741       //
742       else if (recTitle == fgkRecKeys[ kMinPntTrack ]) {
743         if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
744         fMinNPtsPerTrack = recOpt.Atoi();
745       }
746       //
747       else if (recTitle == fgkRecKeys[ kNStDev ]) {
748         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
749         fNStdDev = (Int_t)recOpt.Atof();
750       }
751       //
752       else if (recTitle == fgkRecKeys[ kResCutInit  ]) {
753         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
754         fResCutInitial = recOpt.Atof();
755       }
756       //
757       else if (recTitle == fgkRecKeys[ kResCutOther ]) {
758         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
759         fResCut = recOpt.Atof();
760       }
761       //
762       else if (recTitle == fgkRecKeys[ kLocalSigFactor ]) { //-------------------------
763         for (irec=0;irec<3;irec++) if (nrecElems>irec+1) {
764             fSigmaFactor[irec] = ((TObjString*)recArr->At(irec+1))->GetString().Atof();
765             if (fSigmaFactor[irec]<=0.) stopped = kTRUE;
766           }
767         if (stopped) break; 
768       }
769       //
770       else if (recTitle == fgkRecKeys[ kStartFactor ]) {        //-------------------------
771         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
772         fStartFac = recOpt.Atof();
773       }
774       //
775       else if (recTitle == fgkRecKeys[ kFinalFactor ]) {        //-------------------------
776         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
777         fFinalFac = recOpt.Atof();
778       }
779       //
780       // pepo2708909
781       else if (recTitle == fgkRecKeys[ kExtraClustersMode ]) {        //-------------------------
782         if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
783         fExtraClustersMode = recOpt.Atoi();
784       }
785       // endpepo270809
786       //
787       else if (recTitle == fgkRecKeys[ kBField ]) {         //-------------------------
788         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
789         SetBField( recOpt.Atof() );
790       }
791       //
792       else if (recTitle == fgkRecKeys[ kSDDVDCorrMult ]) {         //-------------------------
793         SetSDDVDCorrMult( recOpt.IsNull() || (recOpt.IsFloat() && (recOpt.Atof())>-0.5) ); 
794       }
795       //
796       else if (recTitle == fgkRecKeys[ kWeightPt ]) {         //-------------------------
797         double wgh = 1;
798         if (!recOpt.IsNull()) {
799           if (!recOpt.IsFloat()) {stopped = kTRUE; break;}
800           else wgh = recOpt.Atof();
801         }
802         SetWeightPt(wgh);
803       }
804       //
805       else if (recTitle == fgkRecKeys[ kSparseMatrix ]) {   // matrix solver type
806         //
807         AliMillePede2::SetGlobalMatSparse(kTRUE);
808         if (recOpt.IsNull()) continue;
809         // solver type and settings
810         if      (recOpt == "MINRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolMinRes );
811         else if (recOpt == "FGMRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolFGMRes );
812         else {stopped = kTRUE; break;}
813         //
814         if (nrecElems>=3) { // preconditioner type
815           recExt = recArr->At(2)->GetName();
816           if (!recExt.IsDigit()) {stopped = kTRUE; break;}
817           AliMillePede2::SetMinResPrecondType( recExt.Atoi() );
818         }
819         //
820         if (nrecElems>=4) { // tolerance
821           recExt = recArr->At(3)->GetName();
822           if (!recExt.IsFloat()) {stopped = kTRUE; break;}
823           AliMillePede2::SetMinResTol( recExt.Atof() );
824         }
825         //
826         if (nrecElems>=5) { // maxIter
827           recExt = recArr->At(4)->GetName();
828           if (!recExt.IsDigit()) {stopped = kTRUE; break;}
829           AliMillePede2::SetMinResMaxIter( recExt.Atoi() );
830         }       
831       }
832       //
833       else if (recTitle == fgkRecKeys[ kRequirePoint ]) {       //-------------------------
834         // syntax:   REQUIRE_POINT where ndet updw nreqpts
835         //    where = LAYER or DETECTOR
836         //    ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
837         //    updw = 1 for Y>0, -1 for Y<0, 0 if not specified
838         //    nreqpts = minimum number of points of that type
839         if (nrecElems>=5) {
840           recOpt.ToUpper();
841           int lr = ((TObjString*)recArr->At(2))->GetString().Atoi() - 1;
842           int hb = ((TObjString*)recArr->At(3))->GetString().Atoi();
843           int np = ((TObjString*)recArr->At(4))->GetString().Atoi();
844           //
845           int rtp = -1; // use for run type
846           if (nrecElems>5) {
847             TString tpstr = ((TObjString*)recArr->At(5))->GetString();
848             if ( tpstr.Contains("cosmics",TString::kIgnoreCase) ) rtp = kCosmics;
849             else if ( tpstr.Contains("collision",TString::kIgnoreCase) ) rtp = kCollision;
850             else {stopped = kTRUE; break;}
851           }
852           //
853           int tpmn= rtp<0 ? 0 : rtp;
854           int tpmx= rtp<0 ? kNDataType-1 : rtp;
855           for (int itp=tpmn;itp<=tpmx;itp++) {
856             fRequirePoints[itp]=kTRUE;
857             if (recOpt == "LAYER") {
858               if (lr<0 || lr>5) {stopped = kTRUE; break;}
859               if (hb>0) fNReqLayUp[itp][lr]=np;
860               else if (hb<0) fNReqLayDown[itp][lr]=np;
861               else fNReqLay[itp][lr]=np;
862             }
863             else if (recOpt == "DETECTOR") {
864               if (lr<0 || lr>2) {stopped = kTRUE; break;}
865               if (hb>0) fNReqDetUp[itp][lr]=np;
866               else if (hb<0) fNReqDetDown[itp][lr]=np;
867               else fNReqDet[itp][lr]=np;
868             }
869             else {stopped = kTRUE; break;}
870           }
871           if (stopped) break;
872         }
873         else {stopped = kTRUE; break;}
874       }
875       //
876       // global constraints on the subunits/orphans 
877       else if (recTitle == fgkRecKeys[ kConstrOrphans ]) {    //------------------------
878         // expect CONSTRAINT_ORPHANS MEAN/MEDIAN Value parID0 ... parID1 ...
879         if (nrecElems<4) {stopped = kTRUE; break;}
880         recExt = recArr->At(2)->GetName();
881         if (!recExt.IsFloat()) {stopped = kTRUE; break;}
882         double val = recExt.Atof();
883         UInt_t pattern = 0;
884         for (irec=3;irec<nrecElems;irec++) { // read params to constraint
885           recExt = recArr->At(irec)->GetName();
886           if (!recExt.IsDigit()) {stopped = kTRUE; break;}
887           pattern |= 0x1 << recExt.Atoi();
888         }
889         if (stopped) break;
890         if      (recOpt == "MEAN")   ConstrainOrphansMean(val,pattern);
891         else if (recOpt == "MEDIAN") ConstrainOrphansMedian(val,pattern);
892         else {stopped = kTRUE; break;}
893       }
894       //
895       else if (recTitle == fgkRecKeys[ kConstrSubunits ]) {    //------------------------
896         // expect CONSTRAINT_SUBUNITS MEAN/MEDIAN Value parID0 ... parID1 ... VolID1 ... VolIDn - VolIDm
897         if (nrecElems<5) {stopped = kTRUE; break;}
898         recExt = recArr->At(2)->GetName();
899         if (!recExt.IsFloat()) {stopped = kTRUE; break;}
900         double val = recExt.Atof();
901         UInt_t pattern = 0;
902         for (irec=3;irec<nrecElems;irec++) { // read params to constraint
903           recExt = recArr->At(irec)->GetName();
904           if (!recExt.IsDigit()) {stopped = kTRUE; break;}
905           int parid = recExt.Atoi();
906           if (parid<kMaxITSSensID) pattern |= 0x1 << recExt.Atoi();
907           else break;           // list of params is over 
908         }
909         if (stopped) break;
910         //
911         Bool_t meanC;
912         if      (recOpt == "MEAN")   meanC = kTRUE;
913         else if (recOpt == "MEDIAN") meanC = kFALSE;
914         else    {stopped = kTRUE; break;}
915         //
916         int curID = -1;
917         int rangeStart = -1;
918         for (;irec<nrecElems;irec++) { // read modules to apply this constraint
919           recExt = recArr->At(irec)->GetName();
920           if (recExt == "-") {rangeStart = curID; continue;}  // range is requested
921           else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
922           else curID = recExt.Atoi();
923           //
924           if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
925           // this was a range start or single 
926           int start;
927           if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
928           else start = curID;  // create constraint either for single module (or 1st in the range)
929           for (int id=start;id<=curID;id++) {
930             int id0 = IsVIDDefined(id);
931             if (id0<0) {AliDebug(3,Form("Undefined module %d requested in the SubUnits constraint, skipping",id)); continue;}
932             if (meanC) ConstrainModuleSubUnitsMean(id0,val,pattern);
933             else       ConstrainModuleSubUnitsMedian(id0,val,pattern);
934           }
935         }
936         if (rangeStart>=0) stopped = kTRUE; // unfinished range
937         if (stopped) break;
938       } 
939       // 
940       // association of modules with local constraints
941       else if (recTitle == fgkRecKeys[ kApplyConstr ]) {            //------------------------
942         // expect APPLY_CONSTRAINT NAME [NAME1...] [VolID1 ... VolIDn - VolIDm]
943         if (nrecElems<3) {stopped = kTRUE; break;}
944         int nmID0=-1,nmID1=-1;
945         for (irec=1;irec<nrecElems;irec++) { // find the range of constraint names
946           recExt = recArr->At(irec)->GetName();
947           if (recExt.IsFloat()) break;
948           // check if such a constraint was declared
949           if (!GetConstraint(recExt.Data())) {
950             AliInfo(Form("No Local constraint %s was declared",recExt.Data())); 
951             stopped=kTRUE; 
952             break;
953           }
954           if (nmID0<0) nmID0 = irec;
955           nmID1 = irec;
956         }
957         if (stopped) break;
958         //
959         if (irec>=nrecElems) {stopped = kTRUE; break;} // no modules provided
960         //
961         // now read the list of modules to constrain
962         int curID = -1;
963         int rangeStart = -1;
964         for (;irec<nrecElems;irec++) { // read modules to apply this constraint
965           recExt = recArr->At(irec)->GetName();
966           if (recExt == "-") {rangeStart = curID; continue;}  // range is requested
967           else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
968           else curID = recExt.Atoi();
969           //
970           if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
971           //
972           // this was a range start or single 
973           int start;
974           if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
975           else start = curID;  // create constraint either for single module (or 1st in the range)
976           for (int id=start;id<=curID;id++) {
977             AliITSAlignMille2Module *md = GetMilleModuleByVID(id);
978             if (!md) {AliDebug(3,Form("Undefined module %d requested in the Local constraint, skipping",id)); continue;}
979             for (int nmid=nmID0;nmid<=nmID1;nmid++) 
980               ((AliITSAlignMille2ConstrArray*)GetConstraint(recArr->At(nmid)->GetName()))->AddModule(md);
981           }
982         }
983         if (rangeStart>=0) stopped = kTRUE; // unfinished range
984         if (stopped) break;
985       }
986       //
987       // request of the same T0 for group of SDD modules
988       else if (recTitle == fgkRecKeys[ kSameSDDT0 ]) {            //------------------------
989         // expect SET_SAME_SDDT0 [SensID1 ... SensIDn - SensIDm]
990         if (nrecElems<3) {stopped = kTRUE; break;}
991         //
992         // now read the list of modules to constrain
993         int curID = -1;
994         int rangeStart = -1;
995         AliITSAlignMille2ConstrArray *cstrT0 = new AliITSAlignMille2ConstrArray("SDDT0",0,0,0,0);
996         int naddM = 0;
997         cstrT0->SetPattern(BIT(AliITSAlignMille2Module::kDOFT0));
998         for (irec=1;irec<nrecElems;irec++) { // read modules to apply this constraint
999           recExt = recArr->At(irec)->GetName();
1000           if (recExt == "-") {rangeStart = curID; continue;}  // range is requested
1001           else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
1002           else curID = recExt.Atoi();
1003           //
1004           if (curID<kSDDoffsID || curID>=kSDDoffsID+kNSDDmod) {stopped = kTRUE; break;}
1005           //
1006           // this was a range start or single 
1007           int start;
1008           if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
1009           else start = curID;  // create constraint either for single module (or 1st in the range)
1010           for (int id=start;id<=curID;id++) {
1011             int vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
1012             if (vid<=1) {AliDebug(3,Form("Undefined module index %d requested in the SAME_SDDT0 constraint, skipping",id)); continue;}
1013             AliITSAlignMille2Module *md = GetMilleModuleByVID(vid);
1014             if (!md) {AliDebug(3,Form("Undefined module %d requested in the Local constraint, skipping",id)); continue;}
1015             cstrT0->AddModule(md,kFALSE);
1016             naddM++;
1017           }       
1018         }
1019         if (rangeStart>=0) stopped = kTRUE; // unfinished range
1020         if (stopped) break;
1021         if (naddM<2) delete cstrT0;
1022         else {
1023           cstrT0->SetConstraintID(GetNConstraints());
1024           fConstraints.Add(cstrT0);
1025         }
1026       }
1027       //
1028       // Do we use new local Y errors?
1029       else if (recTitle == fgkRecKeys[ kUseLocalYErr ]) {
1030         // expect SET_TPAFITTER 
1031         fUseLocalYErr = kTRUE;
1032       }
1033       //
1034       else if (recTitle == fgkRecKeys[ kMinPointsSens ]) {         //-------------------------
1035         if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
1036         SetMinPointsPerSensor( recOpt.Atoi() );
1037       }
1038       //
1039       else if (recTitle == fgkRecKeys[ kOCDBSpecificPath ]) {         //-------------------------
1040         if (recOpt.IsNull() || nrecElems<3 ) {stopped = kTRUE; break;}
1041         AliCDBManager::Instance()->SetSpecificStorage(recOpt.Data(), gSystem->ExpandPathName(recArr->At(2)->GetName()));
1042         AliInfo(Form("Configuration sets OCDB specific storage %s to %s",recOpt.Data(),recArr->At(2)->GetName()));
1043         TObjString *pths = (TObjString*)AliCDBManager::Instance()->GetStorageMap()->GetValue(recOpt.Data());
1044         if (!pths) { stopped = kTRUE; break; }
1045         pths->SetUniqueID(1); // mark as set by user
1046       }
1047       //
1048       else if (recTitle == fgkRecKeys[ kCorrectDiamond ] && fUseDiamond) {         //-------------------------
1049         if (nrecElems<4) {stopped = kTRUE; break;}
1050         for (int i=0;i<3;i++) fCorrDiamond[i] = ((TObjString*)recArr->At(i+1))->GetString().Atof();
1051         AliInfo(Form("Correction %+.4f %+.4f %+.4f will be applied to diamond",fCorrDiamond[0],fCorrDiamond[1],fCorrDiamond[2]));
1052       }
1053       //
1054       else continue; // already processed record
1055       //
1056     } // end of while loop 4 over the various params 
1057     //
1058     break;
1059   } // end of while(1) loop 
1060   //
1061   fclose(pfc);
1062   if (!fDiamondPath.IsNull() && IsDiamondUsed() && LoadDiamond(fDiamondPath) ) stopped = kTRUE;
1063   if (stopped) {
1064     AliError(Form("Failed on record %s %s ...\n",recTitle.Data(),recOpt.Data()));
1065     return -1;
1066   }
1067   //
1068   if (CacheMatricesCurr()) return -1;
1069   SetUseLocalYErrors(fUseLocalYErr); // YErr used only with TPAFitter 
1070   fSegmentationSDD = new AliITSsegmentationSDD();
1071   //
1072   fIsConfigured = kTRUE;
1073   return 0;
1074 }
1075
1076 //________________________________________________________________________________________________________
1077 void AliITSAlignMille2::BuildHierarchy()
1078 {
1079   // build the hieararhy of the modules to align
1080   //
1081   if (!GetUseGlobalDelta() && PseudoParentsAllowed()) {
1082     AliInfo("PseudoParents mode is allowed only when the deltas are global\n"
1083             "Since Deltas are local, switching to NoPseudoParents");
1084     SetAllowPseudoParents(kFALSE);
1085   }
1086   // set parent/child relationship for modules to align
1087   AliInfo("Setting parent/child relationships\n");
1088   //
1089   // 1) child -> parent reference
1090   for (int ipar=0;ipar<fNModules;ipar++) {
1091     AliITSAlignMille2Module* parent = GetMilleModule(ipar);
1092     if (parent->IsSensor()) continue; // sensor cannot be a parent
1093     //
1094     for (int icld=0;icld<fNModules;icld++) {
1095       if (icld==ipar) continue;
1096       AliITSAlignMille2Module* child = GetMilleModule(icld);
1097       if (!child->BelongsTo(parent)) continue;
1098       // child cannot have more sensors than the parent
1099       if (child->GetNSensitiveVolumes() > parent->GetNSensitiveVolumes()) continue;
1100       //
1101       AliITSAlignMille2Module* parOld = child->GetParent();
1102       // is this parent candidate closer than the old parent ? 
1103       if (parOld && parOld->GetNSensitiveVolumes()<parent->GetNSensitiveVolumes()) continue; // parOld is closer
1104       child->SetParent(parent);
1105     }
1106     //
1107   }
1108   //
1109   // add parent -> children reference
1110   for (int icld=0;icld<fNModules;icld++) {
1111     AliITSAlignMille2Module* child = GetMilleModule(icld);
1112     AliITSAlignMille2Module* parent = child->GetParent();
1113     if (parent) parent->AddChild(child);
1114   }  
1115   //
1116   // reorder the modules in such a way that parents come first
1117   for (int icld=0;icld<fNModules;icld++) {
1118     AliITSAlignMille2Module* child  = GetMilleModule(icld);
1119     AliITSAlignMille2Module* parent; 
1120     while ( (parent=child->GetParent()) &&  (parent->GetUniqueID()>child->GetUniqueID()) ) {
1121       // swap
1122       fMilleModule[icld] = parent;
1123       fMilleModule[parent->GetUniqueID()] = child;
1124       child->SetUniqueID(parent->GetUniqueID());
1125       parent->SetUniqueID(icld);
1126       child = parent;
1127     }
1128     //
1129   }  
1130   //
1131   // Go over the child->parent chain and mark modules with explicitly provided sensors.
1132   // If the sensors of the unit are explicitly declared, all undeclared sensors are 
1133   // suppresed in this unit.
1134   for (int icld=fNModules;icld--;) {
1135     AliITSAlignMille2Module* child = GetMilleModule(icld);
1136     AliITSAlignMille2Module* parent = child->GetParent();
1137     if (!parent) continue;
1138     //
1139     // check if this parent was already processed
1140     if (!parent->AreSensorsProvided()) {
1141       parent->DelSensitiveVolumes();
1142       parent->SetSensorsProvided(kTRUE);
1143     }
1144     // reattach sensors to parent
1145     for (int isc=child->GetNSensitiveVolumes();isc--;) {
1146       UShort_t senVID = child->GetSensVolVolumeID(isc);
1147       if (!parent->IsIn(senVID)) parent->AddSensitiveVolume(senVID);
1148     }
1149   }
1150   //
1151 }
1152
1153 // pepo
1154 //________________________________________________________________________________________________________
1155 void AliITSAlignMille2::SetCurrentModule(Int_t id)
1156 {
1157   // set the current supermodule
1158   // new meaning
1159   if (fMilleVersion>=2) {
1160     fCurrentModule = GetMilleModule(id);
1161     return;
1162   }
1163   // old meaning
1164   if (fMilleVersion<=1) {
1165     Int_t index=id;
1166     /// set as current the SuperModule that contains the 'index' sens.vol.
1167     if (index<0 || index>2197) {
1168       AliInfo("index does not correspond to a sensitive volume!");
1169       return;
1170     }
1171     UShort_t voluid=AliITSAlignMille2Module::GetVolumeIDFromIndex(index);
1172     Int_t k=IsContained(voluid);
1173     if (k>=0){
1174       fCurrentSensID = index;
1175       fCluster.SetVolumeID(voluid);
1176       fCluster.SetXYZ(0,0,0);
1177       InitModuleParams();
1178     }
1179     else
1180       AliInfo(Form("module %d not defined\n",index));    
1181   }
1182 }
1183 // endpepo
1184 //________________________________________________________________________________________________________
1185 void AliITSAlignMille2::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts,Int_t runtype) 
1186 {
1187   // set minimum number of points in specific detector or layer
1188   // where = LAYER or DETECTOR
1189   // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
1190   // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
1191   // nreqpts = minimum number of points of that type
1192   ndet--;
1193   int tpmn= runtype<0 ? 0 : runtype;
1194   int tpmx= runtype<0 ? kNDataType-1 : runtype;
1195   //
1196   for (int itp=tpmn;itp<=tpmx;itp++) {
1197     fRequirePoints[itp]=kTRUE;
1198     if (strstr(where,"LAYER")) {
1199       if (ndet<0 || ndet>5) return;
1200       if (updw>0) fNReqLayUp[itp][ndet]=nreqpts;
1201       else if (updw<0) fNReqLayDown[itp][ndet]=nreqpts;
1202       else fNReqLay[itp][ndet]=nreqpts;
1203     }
1204     else if (strstr(where,"DETECTOR")) {
1205       if (ndet<0 || ndet>2) return;
1206       if (updw>0) fNReqDetUp[itp][ndet]=nreqpts;
1207       else if (updw<0) fNReqDetDown[itp][ndet]=nreqpts;
1208       else fNReqDet[itp][ndet]=nreqpts; 
1209     }
1210   }
1211 }
1212
1213 //________________________________________________________________________________________________________
1214 Int_t AliITSAlignMille2::GetModuleIndex(const Char_t *symname) 
1215 {
1216   /// index from symname
1217   if (!symname) return -1;
1218   for (Int_t i=0;i<=kMaxITSSensID; i++) {
1219     if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i;
1220   }
1221   return -1;
1222 }
1223
1224 //________________________________________________________________________________________________________
1225 Int_t AliITSAlignMille2::GetModuleIndex(UShort_t voluid) 
1226 {
1227   /// index from volume ID
1228   AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
1229   if (lay<1|| lay>6) return -1;
1230   Int_t idx=Int_t(voluid)-2048*lay;
1231   if (idx>=AliGeomManager::LayerSize(lay)) return -1;
1232   for (Int_t ilay=1; ilay<lay; ilay++) 
1233     idx += AliGeomManager::LayerSize(ilay);
1234   return idx;
1235 }
1236
1237 //________________________________________________________________________________________________________
1238 UShort_t AliITSAlignMille2::GetModuleVolumeID(const Char_t *symname) 
1239 {
1240   /// volume ID from symname
1241   /// works for sensitive volumes only
1242   if (!symname) return 0;
1243
1244   for (UShort_t voluid=2000; voluid<13300; voluid++) {
1245     Int_t modId;
1246     AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
1247     if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
1248       if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
1249     }
1250   }
1251
1252   return 0;
1253 }
1254
1255 //________________________________________________________________________________________________________
1256 UShort_t AliITSAlignMille2::GetModuleVolumeID(Int_t index) 
1257 {
1258   /// volume ID from index
1259   if (index<0) return 0;
1260   if (index<2198)
1261     return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index));
1262   else {
1263     for (int i=0; i<fNSuperModules; i++) {
1264       if (GetSuperModule(i)->GetIndex()==index) return GetSuperModule(i)->GetVolumeID();
1265     }
1266   }
1267   return 0;
1268 }
1269
1270 //________________________________________________________________________________________________________
1271 Int_t AliITSAlignMille2::LoadGeometry(TString& path) 
1272 {
1273   // initialize ideal geometry
1274   AliInfo(Form("Loading ideal geometry %s",path.Data()));
1275   if (path.IsNull()) {
1276     AliError("Path to geometry is not provided");
1277     return -1;
1278   }
1279   //
1280   AliCDBEntry *entry = 0;
1281   TGeoManager *gm = 0;
1282   while(1) {
1283     if (path.BeginsWith("path: ")) { // must load from OCDB
1284       entry = GetCDBEntry(path.Data());
1285       if (!entry) break;
1286       gm = (TGeoManager*) entry->GetObject();
1287       entry->SetObject(NULL);
1288       entry->SetOwner(kTRUE);
1289       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
1290       //      delete cdbId;
1291       //      delete entry;
1292       break;
1293     }
1294     //
1295     if (gSystem->AccessPathName(path.Data())) break;
1296     TFile* precf = TFile::Open(path.Data());
1297     if (precf->FindKey("ALICE")) gm = (TGeoManager*)precf->Get("ALICE");
1298     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
1299       gm = (TGeoManager*) entry->GetObject();
1300       if (gm && gm->InheritsFrom(TGeoManager::Class())) entry->SetObject(NULL);
1301       else gm = 0;
1302       entry->SetObject(NULL);
1303       entry->SetOwner(kTRUE);
1304       delete entry;
1305     }
1306     precf->Close();
1307     delete precf;
1308     break;
1309   } 
1310   //
1311   if (!gm) {AliError(Form("Failed to load geometry from %s",path.Data())); return -1;}
1312   AliGeomManager::SetGeometry(gm);
1313   fGeoManager = AliGeomManager::GetGeometry();
1314   if (!fGeoManager) {
1315     AliInfo("Couldn't initialize geometry");
1316     return -1;
1317   }
1318   return 0;
1319 }
1320
1321 //________________________________________________________________________________________________________
1322 Int_t AliITSAlignMille2::SetConstraintWrtRef(const char* reffname) 
1323 {
1324   // Load the global deltas from this file. The local gaussian constraints on some modules 
1325   // will be defined with respect to the deltas from this reference file, converted to local
1326   // delta format. Note: conversion to local format requires reloading the geometry!
1327   //
1328   AliInfo(Form("Loading reference deltas for local constraints from %s",reffname));
1329   if (!fGeoManager) return -1; 
1330   fConstrRefPath = reffname;
1331   if (fConstrRefPath == "IDEAL") { // the reference is the ideal geometry, just create dummy reference array
1332     fConstrRef = new TClonesArray("AliAlignObjParams",1);
1333     return 0;
1334   }
1335   if (LoadDeltas(fConstrRefPath,fConstrRef)) return -1;
1336   //
1337   // we need ideal geometry to convert global deltas to local ones
1338   if (fUsePreAlignment) {
1339     AliError("The call of SetConstraintWrtRef must be done before application of the prealignment");
1340     return -1;
1341   }
1342   //
1343   AliInfo("Converting global reference deltas to local ones");
1344   Int_t nprea = fConstrRef->GetEntriesFast();
1345   for (int ix=0; ix<nprea; ix++) {
1346     AliAlignObjParams *preo=(AliAlignObjParams*) fConstrRef->At(ix);
1347     if (!preo->ApplyToGeometry()) return -1;
1348   }
1349   //
1350   // now convert the global reference deltas to local ones
1351   for (int i=fConstrRef->GetEntriesFast();i--;) {
1352     AliAlignObjParams *preo = (AliAlignObjParams*)fConstrRef->At(i);
1353     TGeoHMatrix * mupd = AliGeomManager::GetMatrix(preo->GetSymName());
1354     if (!mupd) {  // this is not alignable entry, need to look in the supermodules
1355       for (int im=fNSuperModules;im--;) {
1356         AliITSAlignMille2Module* mod = GetSuperModule(im);
1357         if ( strcmp(mod->GetName(), preo->GetSymName()) ) continue;
1358         mupd = mod->GetMatrix();
1359         break;
1360       }
1361       if (!mupd) {
1362         AliError(Form("Failed to find the volume for reference %s",preo->GetSymName()));
1363         return -1;
1364       }
1365     } 
1366     TGeoHMatrix preMat;
1367     preo->GetMatrix(preMat);                     //  Delta_Glob
1368     TGeoHMatrix tmpMat    = *mupd;               //  Delta_Glob * Delta_Glob_Par * M
1369     preMat.MultiplyLeft( &tmpMat.Inverse() );    //  M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
1370     tmpMat.MultiplyLeft( &preMat );              //  (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
1371     preo->SetMatrix(tmpMat);     // local corrections 
1372   }
1373   //
1374   // we need to reload the geometry spoiled by this reference deltas...
1375   delete fGeoManager;
1376   AliInfo("Reloading target ideal geometry");
1377   return LoadGeometry(fGeometryPath);
1378   //
1379 }
1380
1381 //________________________________________________________________________________________________________
1382 void AliITSAlignMille2::Init()
1383 {
1384   // perform global initialization
1385   //
1386   if (fIsMilleInit) {
1387     AliInfo("Millepede has been already initialized!");
1388     return;
1389   }
1390   // range constraints in such a way that the childs are constrained before their parents
1391   // orphan constraints come last
1392   for (int ic=0;ic<GetNConstraints();ic++) {
1393     for (int ic1=ic+1;ic1<GetNConstraints();ic1++) {
1394       AliITSAlignMille2Constraint *cst0 = GetConstraint(ic);
1395       AliITSAlignMille2Constraint *cst1 = GetConstraint(ic1);
1396       if (cst0->GetModuleID()<cst1->GetModuleID()) {
1397         // swap
1398         fConstraints[ic] = cst1;
1399         fConstraints[ic1] = cst0;
1400       }
1401     }
1402   }
1403   //
1404   if (!GetUseGlobalDelta()) {
1405     AliInfo("ATTENTION: The parameters are defined in the local frame, no check for degeneracy will be done");
1406     for (int imd=fNModules;imd--;) {
1407       AliITSAlignMille2Module* mod = GetMilleModule(imd);
1408       int npar = mod->GetNParTot();
1409       // the parameter may have max 1 free instance, otherwise the equations are underdefined
1410       for (int ipar=0;ipar<npar;ipar++) {
1411         if (!mod->IsFreeDOF(ipar)) continue;
1412         mod->SetParOffset(ipar,fNGlobal++);
1413       }
1414     }
1415   }
1416   else {
1417     // init millepede, decide which parameters are to be fitted explicitly
1418     for (int imd=fNModules;imd--;) {
1419       AliITSAlignMille2Module* mod = GetMilleModule(imd);
1420       if (mod->IsNotInConf()) continue; // dummy module
1421       int npar = mod->GetNParTot();
1422       // the parameter may have max 1 free instance, otherwise the equations are underdefined
1423       for (int ipar=0;ipar<npar;ipar++) {
1424         if (!mod->IsFreeDOF(ipar)) continue;  // fixed
1425         //
1426         int nFreeInstances = 0;
1427         //
1428         AliITSAlignMille2Module* parent = mod;
1429         Bool_t cstMeanMed=kFALSE,cstGauss=kFALSE;
1430         //
1431         Bool_t addToFit = kFALSE;       
1432         // the parameter may be ommitted from explicit fit (if PseudoParentsAllowed is true) if
1433         // 1) it is not explicitly constrained or its does not participate in Gaussian constraint
1434         // 2) the same applies to all of its parents
1435         // 3) it has at least 1 unconstrained direct child
1436         while(parent) {
1437           if (parent->IsNotInConf()) {parent = parent->GetParent(); continue;}
1438           if (!parent->IsFreeDOF(ipar)) {parent = parent->GetParent(); continue;}
1439           nFreeInstances++;
1440           if (IsParModConstrained(parent,ipar, cstMeanMed, cstGauss)) nFreeInstances--;
1441           if (cstGauss) addToFit = kTRUE;
1442           parent = parent->GetParent();
1443         }
1444         if (nFreeInstances>1) {
1445           AliError(Form("Parameter#%d of module %s\nhas %d free instances in the "
1446                         "unconstrained parents\nSystem is undefined",ipar,mod->GetName(),nFreeInstances));
1447           exit(1);
1448         }
1449         //
1450         // i) Are PseudoParents allowed?
1451         if (!PseudoParentsAllowed()) addToFit = kTRUE;
1452         // ii) check if this module has no child with such a free parameter. Since the order of this check 
1453         // goes from child to parent, by this moment such a parameter must have been already added
1454         else if (!IsParModFamilyVaried(mod,ipar))  addToFit = kTRUE;  // no varied children at all
1455         else if (!IsParFamilyFree(mod,ipar,1))     addToFit = kTRUE;  // no unconstrained direct children
1456         // otherwise the value of this parameter can be extracted from simple contraint and the values of 
1457         // the relevant parameters of its children the fit is done. Hence it is not included
1458         if (!addToFit) continue;
1459         //
1460         // shall add this parameter to explicit fit
1461         //      printf("Adding %s %d -> %d\n",mod->GetName(), ipar, fNGlobal);
1462         mod->SetParOffset(ipar,fNGlobal++);
1463       }
1464     }
1465   }
1466   //
1467   AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, kNLocal, fNStdDev));
1468   fGlobalDerivatives = new Double_t[fNGlobal];
1469   memset(fGlobalDerivatives,0,fNGlobal*sizeof(Double_t));
1470   //
1471   fMillepede->InitMille(fNGlobal,kNLocal,fNStdDev,fResCut,fResCutInitial);
1472   fMillepede->SetMinPntValid(fMinPntPerSens);
1473   fIsMilleInit = kTRUE;
1474   //
1475   ResetLocalEquation();    
1476   AliInfo("Parameters initialized to zero");
1477   //
1478   /// Fix non free parameters
1479   for (Int_t i=0; i<fNModules; i++) {
1480     AliITSAlignMille2Module* mod = GetMilleModule(i);
1481     for (Int_t j=0; j<mod->GetNParTot(); j++) {
1482       if (mod->GetParOffset(j)<0) continue; // not varied
1483       FixParameter(mod->GetParOffset(j),mod->GetParConstraint(j));
1484       fMillepede->SetParamGrID(i, mod->GetParOffset(j));
1485     }
1486   }
1487   //
1488   ResetCovIScale();
1489   // Set iterations
1490   if (fStartFac>1) fMillepede->SetIterations(fStartFac);    
1491   if (fFinalFac>1) fMillepede->SetChi2CutRef(fFinalFac);    
1492   fTrackBuff.Expand(24);
1493   //
1494 }
1495
1496 //________________________________________________________________________________________________________
1497 void AliITSAlignMille2::AddConstraint(Double_t *par, Double_t value, Double_t sigma) 
1498 {
1499   /// Constrain equation defined by par to value
1500   if (!fIsMilleInit) Init();
1501   fMillepede->SetGlobalConstraint(par, value, sigma);
1502   AliInfo("Adding constraint");
1503 }
1504
1505 //________________________________________________________________________________________________________
1506 void AliITSAlignMille2::InitGlobalParameters(Double_t *par) 
1507 {
1508   /// Initialize global parameters with par array
1509   if (!fIsMilleInit) Init();
1510   fMillepede->SetGlobalParameters(par);
1511   AliInfo("Init Global Parameters");
1512 }
1513
1514 //________________________________________________________________________________________________________ 
1515 void AliITSAlignMille2::FixParameter(Int_t iPar, Double_t value) 
1516 {
1517   /// Parameter iPar is encourage to vary in [-value;value]. 
1518   /// If value == 0, parameter is fixed
1519   if (!fIsMilleInit) {
1520     AliInfo("Millepede has not been initialized!");
1521     return;
1522   }
1523   fMillepede->SetParSigma(iPar, value);
1524   if (IsZero(value)) AliInfo(Form("Parameter %i Fixed", iPar));
1525 }
1526
1527 //________________________________________________________________________________________________________
1528 void AliITSAlignMille2::ResetLocalEquation()
1529 {
1530   /// Reset the derivative vectors
1531   for(int i=kNLocal;i--;)  fLocalDerivatives[i] = 0.0;
1532   memset(fGlobalDerivatives, 0, fNGlobal*sizeof(double) );
1533 }
1534
1535 //________________________________________________________________________________________________________
1536 Int_t AliITSAlignMille2::ApplyToGeometry() 
1537 {
1538   // apply prealignment to ideal geometry
1539   Int_t nprea = fPrealignment->GetEntriesFast();
1540   AliInfo(Form("Array of prealignment deltas: %d entries",nprea));
1541   //
1542   for (int ix=0; ix<nprea; ix++) {
1543     AliAlignObjParams *preo=(AliAlignObjParams*) fPrealignment->At(ix);
1544     Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(preo->GetVolUID());
1545     if (index>=0) {
1546       if (index>=fPreAlignQF.GetSize()) fPreAlignQF.Set(index+10);
1547       fPreAlignQF[index] = (int) preo->GetUniqueID()+1;
1548     }
1549     if (!preo->ApplyToGeometry()) {
1550       AliError(Form("Failed on ApplyToGeometry at %s",preo->GetSymName()));
1551       return -1;
1552     }
1553   }
1554   //
1555   fUsePreAlignment = kTRUE;
1556   return 0;
1557 }
1558
1559 //________________________________________________________________________________________________________
1560 Int_t AliITSAlignMille2::GetPreAlignmentQualityFactor(Int_t index) const
1561 {
1562   // quality factors from prealignment
1563   if (!fUsePreAlignment || index<0 || index>=fPreAlignQF.GetSize()) return -1;
1564   return fPreAlignQF[index]-1;
1565 }
1566
1567 //________________________________________________________________________________________________________
1568 AliTrackPointArray *AliITSAlignMille2::PrepareTrack(const AliTrackPointArray *atp) 
1569 {
1570   /// create a new AliTrackPointArray keeping only defined modules
1571   /// move points according to a given prealignment, if any
1572   /// sort alitrackpoints w.r.t. global Y direction, if cosmics
1573   const Double_t kRad2L[6] = {5*5,10*10,18*18,30*30,40*40,60*60};
1574   const Float_t kSensSigY2[6] = {200e-4*200e-4/12, 200e-4*200e-4/12, 
1575                                  300e-4*300e-4/12, 300e-4*300e-4/12, 
1576                                  300e-4*300e-4/12, 300e-4*300e-4/12}; // thickness^2/12
1577   //
1578   fTrack = NULL;
1579   Int_t   idx[20] = {0};
1580   Short_t lrID[20] = {0};
1581   Int_t npts=atp->GetNPoints();
1582   if (npts<fMinNPtsPerTrack) return NULL;
1583   TGeoHMatrix hcov;
1584   //
1585   /// checks if AliTrackPoints belong to defined modules
1586   Int_t ngoodpts=0;
1587   Int_t intidx[20];
1588   for (int j=0; j<npts; j++) {
1589     intidx[j] = GetRequestedModID(atp->GetVolumeID()[j]);
1590     if (intidx[j]<0) continue;
1591     ngoodpts++;
1592     Float_t xx=atp->GetX()[j];
1593     Float_t yy=atp->GetY()[j];
1594     Float_t r=xx*xx + yy*yy;
1595     int lay;
1596     for (lay=0;lay<6;lay++) if (r<kRad2L[lay]) break;
1597     if (lay>5) continue;
1598     lrID[j] = lay;
1599   }
1600   //
1601   AliDebug(3,Form("Number of points in defined modules: %d out of %d",ngoodpts,npts));
1602
1603   // pepo270809
1604   Int_t nextra=0;
1605   // extra clusters selection mode  
1606   if (fExtraClustersMode) {
1607     // 1 = keep one cluster, remove randomly the extra
1608     // 2 = keep one cluster, remove the internal one
1609     // 10 = keep tracks only if at least one extra is present
1610     
1611     int iextra1[20],iextra2[20],layovl[20];
1612     // extra clusters mapping
1613     for (Int_t ipt=0; ipt<npts; ipt++) {
1614       if (intidx[ipt]<0) continue; // looks only defined modules...
1615       float p1x=atp->GetX()[ipt];
1616       float p1y=atp->GetY()[ipt];
1617       float p1z=atp->GetZ()[ipt];
1618       int lay1=int(AliGeomManager::VolUIDToLayer(atp->GetVolumeID()[ipt]));
1619       float r1 = p1x*p1x + p1y*p1y;
1620       UShort_t volid1=atp->GetVolumeID()[ipt];
1621       
1622       for (int ik=ipt+1; ik<npts; ik++) {
1623         if (intidx[ik]<0) continue;
1624         // compare point ipt with next ones
1625         int lay2=int(AliGeomManager::VolUIDToLayer(atp->GetVolumeID()[ik]));
1626         // check if same layer
1627         if (lay2 != lay1) continue;
1628         UShort_t volid2=atp->GetVolumeID()[ik];
1629         // check if different module
1630         if (volid1 == volid2) continue;
1631
1632         float p2x=atp->GetX()[ik];
1633         float p2y=atp->GetY()[ik];
1634         float p2z=atp->GetZ()[ik];
1635         float r2 = p2x*p2x + p2y*p2y;   
1636         float dr= (p1x-p2x)*(p1x-p2x) + (p1y-p2y)*(p1y-p2y) + (p1z-p2z)*(p1z-p2z);
1637         
1638         // looks for pairs with dr<1 cm, same layer but different module
1639         if (dr<1.0) {
1640           // extra1 is the one with smaller radius in rphi plane
1641           if (r1<r2) {
1642             iextra1[nextra]=ipt;
1643             iextra2[nextra]=ik;
1644           }
1645           else {
1646             iextra1[nextra]=ik;
1647             iextra2[nextra]=ipt;
1648           }
1649           layovl[nextra]=lay1;    
1650           nextra++;
1651         }
1652       }
1653     } // end overlaps mapping
1654     
1655     // mode=1: keep only one clusters and remove the other randomly
1656     if (fExtraClustersMode==1 && nextra) {
1657       for (int ie=0; ie<nextra; ie++) {
1658         if (gRandom->Rndm()<0.5) 
1659           intidx[iextra1[ie]]=-1;
1660         else
1661           intidx[iextra2[ie]]=-1;         
1662       }
1663     }
1664
1665     // mode=2: keep only one clusters and remove the other...
1666     if (fExtraClustersMode==2 && nextra) {
1667       for (int ie=0; ie<nextra; ie++) {
1668         if (layovl[ie]==1) intidx[iextra2[ie]]=-1;
1669         else if (layovl[ie]==2) intidx[iextra1[ie]]=-1;
1670         else intidx[iextra1[ie]]=-1;      
1671       }
1672     }
1673
1674     // mode=10: reject track if no overlaps are present
1675     if (fExtraClustersMode==10 && nextra==0) {
1676       AliInfo("Track with no extra clusters: rejected!");
1677       return NULL;
1678     }
1679     
1680     // recalculate ngoodpts
1681     ngoodpts=0;
1682     for (int i=0; i<npts; i++) {
1683       if (intidx[i]>=0) ngoodpts++;
1684     }
1685   }
1686   // endpepo270809
1687
1688   // reject track if not enough points are left
1689   if (ngoodpts<fMinNPtsPerTrack) {
1690     AliDebug(2,"Track with not enough points!");
1691     return NULL;
1692   }
1693   // >> RS
1694   AliTrackPoint p;
1695   // check points in specific places
1696   if (fRequirePoints[fDataType]) {
1697     Int_t nlayup[6],nlaydown[6],nlay[6];
1698     Int_t ndetup[3],ndetdown[3],ndet[3];
1699     for (Int_t j=0; j<6; j++) {nlayup[j]=0; nlaydown[j]=0; nlay[j]=0;}
1700     for (Int_t j=0; j<3; j++) {ndetup[j]=0; ndetdown[j]=0; ndet[j]=0;}
1701     
1702     for (int i=0; i<npts; i++) {
1703       // skip not defined points
1704       if (intidx[i]<0) continue;
1705       //      
1706       Float_t yy=atp->GetY()[i];
1707       int lay = lrID[i];
1708       int det=lay/2;
1709       //printf("Point %d - x=%f  y=%f  R=%f  lay=%d  det=%d\n",i,xx,yy,r,lay,det);
1710
1711       if (yy>=0.0) { // UP point
1712         nlayup[lay]++;
1713         nlay[lay]++;
1714         ndetup[det]++;
1715         ndet[det]++;
1716       }
1717       else {
1718         nlaydown[lay]++;
1719         nlay[lay]++;
1720         ndetdown[det]++;
1721         ndet[det]++;
1722       }
1723     }
1724     //
1725     // checks minimum values
1726     Bool_t isok=kTRUE;
1727     for (Int_t j=0; j<6; j++) {
1728       if (nlayup[j]<fNReqLayUp[fDataType][j]) isok=kFALSE; 
1729       if (nlaydown[j]<fNReqLayDown[fDataType][j]) isok=kFALSE; 
1730       if (nlay[j]<fNReqLay[fDataType][j]) isok=kFALSE; 
1731     }
1732     for (Int_t j=0; j<3; j++) {
1733       if (ndetup[j]<fNReqDetUp[fDataType][j]) isok=kFALSE; 
1734       if (ndetdown[j]<fNReqDetDown[fDataType][j]) isok=kFALSE; 
1735       if (ndet[j]<fNReqDet[fDataType][j]) isok=kFALSE; 
1736     }
1737     if (!isok) {
1738       AliDebug(2,Form("Track does not meet all location point requirements!"));
1739       return NULL;
1740     }
1741   }
1742   // build a new track with (sorted) (prealigned) good points
1743   // pepo200709
1744   //fTrack = (AliTrackPointArray*)fTrackBuff[ngoodpts-fMinNPtsPerTrack];
1745   Int_t addVertex = IsTypeCollision()&&((fUseDiamond&&(fCheckDiamondPoint!=kDiamondIgnore))||(fUseVertex&&fVertexSet)) ? 1 : 0;
1746   fTrack = (AliTrackPointArray*)fTrackBuff[ngoodpts + addVertex ];
1747   if (!fTrack) {
1748     fTrack = new AliTrackPointArray(ngoodpts + addVertex);
1749     //    fTrackBuff.AddAtAndExpand(fTrack,ngoodpts-fMinNPtsPerTrack);
1750     fTrackBuff.AddAtAndExpand(fTrack,ngoodpts + addVertex);
1751   }  
1752   //  fTrack = new AliTrackPointArray(ngoodpts);
1753   // endpepo200709
1754   //
1755   //
1756   for (int i=0; i<npts; i++) idx[i]=i;
1757   // sort track if required
1758   if (IsTypeCosmics()) TMath::Sort(npts,atp->GetY(),idx); // sort descending...
1759   //
1760   Int_t npto=0;
1761   if (fClusLoc.GetSize()<3*npts)    fClusLoc.Set(3*npts);
1762   if (fClusGlo.GetSize()<3*npts)    fClusGlo.Set(3*npts);
1763   if (fClusSigLoc.GetSize()<3*npts) fClusSigLoc.Set(3*npts);
1764   //
1765   for (int i=0; i<npts; i++) {
1766     // skip not defined points
1767     if (intidx[idx[i]]<0) continue;
1768     atp->GetPoint(p,idx[i]);
1769     int sid = AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());
1770     //
1771     // prealign point if required
1772     // get matrix used to produce the digits
1773     AliITSAlignMille2Module *mod = GetMilleModule(intidx[idx[i]]);
1774     TGeoHMatrix *svOrigMatrix = GetSensorOrigMatrixSID(sid); //mod->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
1775     // get back real local coordinate
1776     fMeasLoc  = fClusLoc.GetArray() + npto*3;
1777     fMeasGlo  = fClusGlo.GetArray() + npto*3;
1778     fSigmaLoc = fClusSigLoc.GetArray() + npto*3;
1779     fMeasGlo[0]=p.GetX();
1780     fMeasGlo[1]=p.GetY();
1781     fMeasGlo[2]=p.GetZ();
1782     AliDebug(3,Form("Global coordinates of measured point : X=%+f  Y=%+f  Z=%+f \n",fMeasGlo[0],fMeasGlo[1],fMeasGlo[2]));
1783     svOrigMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
1784     AliDebug(3,Form("Local coordinates of measured point : X=%+f  Y=%+f  Z=%+f \n",fMeasLoc[0],fMeasLoc[1],fMeasLoc[2]));
1785     //
1786     if (p.GetDriftTime()>0) ProcessSDDPointInfo(&p,sid, npto);     // for SDD points extract vdrift
1787     //
1788     // update covariance matrix
1789     Double_t hcovel[9];
1790     hcovel[0]=double(p.GetCov()[0]);
1791     hcovel[1]=double(p.GetCov()[1]);
1792     hcovel[2]=double(p.GetCov()[2]);
1793     hcovel[3]=double(p.GetCov()[1]);
1794     hcovel[4]=double(p.GetCov()[3]);
1795     hcovel[5]=double(p.GetCov()[4]);
1796     hcovel[6]=double(p.GetCov()[2]);
1797     hcovel[7]=double(p.GetCov()[4]);
1798     hcovel[8]=double(p.GetCov()[5]);
1799     hcov.SetRotation(hcovel);
1800     //
1801     if (AliLog::GetGlobalDebugLevel()>=2) {
1802       AliInfo("Original Global Cov Matrix");
1803       printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovel[0],hcovel[1],hcovel[2],hcovel[4],hcovel[5],hcovel[8]);
1804     } 
1805     //
1806     // now rotate in local system
1807     hcov.Multiply(svOrigMatrix);
1808     hcov.MultiplyLeft(&svOrigMatrix->Inverse());
1809     // now hcov is LOCAL COVARIANCE MATRIX
1810     // apply sigma scaling
1811     Double_t *hcovscl = hcov.GetRotationMatrix();
1812     /*
1813     const float *cv = p.GetCov();
1814     printf("## %d %d  %+.3e %+.3e %+.3e   %+.3e %+.3e %+.3e   %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e   %+.3e %+.3e %+.3e\n",
1815            sid,p.GetClusterType(),
1816            fMeasGlo[0],fMeasGlo[1],fMeasGlo[2],
1817            fMeasLoc[0],fMeasLoc[1],fMeasLoc[2],
1818            cv[0],cv[1],cv[2],cv[3],cv[4],cv[5],
1819            hcovscl[0],hcovscl[4],hcovscl[8]);
1820
1821     */
1822     if (AliLog::GetGlobalDebugLevel()>=2) {
1823       AliInfo("Original Local Cov Matrix");
1824       printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1825     } 
1826     hcovscl[4] = fUseLocalYErr ? kSensSigY2[lrID[idx[i]]] : 1E-8; // error due to the sensor thickness
1827     //
1828     for (int ir=3;ir--;) for (int ic=3;ic--;) {
1829         if (ir==ic) {     
1830           if ( IsZero(hcovscl[ir*3+ic],1e-8) ) hcovscl[ir*3+ic] = 1E-8;
1831           else hcovscl[ir*3+ic] *= mod->GetSigmaFactor(ir)*mod->GetSigmaFactor(ic); //RRR
1832           fSigmaLoc[ir] = TMath::Sqrt(hcovscl[ir*3+ic]);
1833         }
1834         else hcovscl[ir*3+ic]  = 0;
1835       }
1836     //
1837     if (AliLog::GetGlobalDebugLevel()>=2) {
1838       AliInfo("Modified Local Cov Matrix");
1839       printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1840     } 
1841     //
1842     if (fBug==1) {
1843       // correzione bug LAYER 5  SSD temporanea..
1844       int ssdidx=AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());
1845       if (ssdidx>=500 && ssdidx<1248) {
1846         int ladder=(ssdidx-500)%22;
1847         if (ladder==18) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx+1));
1848         if (ladder==19) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx-1));
1849       }
1850     }
1851     /// get (evenctually prealigned) matrix of sens. vol.
1852     TGeoHMatrix *svMatrix = GetSensorCurrMatrixSID(sid);    //mod->GetSensitiveVolumeMatrix(p.GetVolumeID());
1853     // modify global coordinates according with pre-aligment
1854     svMatrix->LocalToMaster(fMeasLoc,fMeasGlo);
1855     // now rotate in local system
1856     hcov.Multiply(&svMatrix->Inverse());
1857     hcov.MultiplyLeft(svMatrix);         // hcov is back in GLOBAL RF
1858     // cure once more
1859     for (int ir=3;ir--;) for (int ic=3;ic--;) if (IsZero(hcovscl[ir*3+ic])) hcovscl[ir*3+ic] = 0.;
1860     //    printf("\nErrMatGlob: after\n"); hcov.Print(""); //RRR
1861     //
1862     if (AliLog::GetGlobalDebugLevel()>=2) {
1863       AliInfo("Modified Global Cov Matrix");
1864       printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1865     } 
1866     //
1867     Float_t pcov[6];
1868     pcov[0]=hcovscl[0];
1869     pcov[1]=hcovscl[1];
1870     pcov[2]=hcovscl[2];
1871     pcov[3]=hcovscl[4];
1872     pcov[4]=hcovscl[5];
1873     pcov[5]=hcovscl[8];
1874     // 
1875     // make sure the matrix is positive definite
1876     {
1877       enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5};
1878       if (pcov[kXX]*pcov[kYY]*0.999<pcov[kXY]*pcov[kXY]) pcov[kXY] = 0.999*TMath::Sign((float)TMath::Sqrt(pcov[kXX]*pcov[kYY]),pcov[kXY]);
1879       if (pcov[kXX]*pcov[kZZ]*0.999<pcov[kXZ]*pcov[kXZ]) pcov[kXZ] = 0.999*TMath::Sign((float)TMath::Sqrt(pcov[kXX]*pcov[kZZ]),pcov[kXZ]);
1880       if (pcov[kYY]*pcov[kZZ]*0.999<pcov[kYZ]*pcov[kYZ]) pcov[kYZ] = 0.999*TMath::Sign((float)TMath::Sqrt(pcov[kYY]*pcov[kZZ]),pcov[kYZ]);
1881     }
1882     //
1883     p.SetXYZ(fMeasGlo[0],fMeasGlo[1],fMeasGlo[2],pcov);
1884     //    printf("New Gl coordinates of measured point : X=%f  Y=%f  Z=%f \n",fMeasGlo[0],fMeasGlo[1],fMeasGlo[2]);
1885     AliDebug(3,Form("New global coordinates of measured point : X=%+f  Y=%+f  Z=%+f \n",fMeasGlo[0],fMeasGlo[1],fMeasGlo[2]));
1886     fTrack->AddPoint(npto,&p);
1887     AliDebug(2,Form("Adding point[%d] = ( %+f , %+f , %+f )     volid = %d",npto,fTrack->GetX()[npto],
1888                     fTrack->GetY()[npto],fTrack->GetZ()[npto],fTrack->GetVolumeID()[npto] ));
1889     //    printf("Adding %d %d %f\n",npto, p.GetVolumeID(), p.GetY()); 
1890     npto++;
1891   }
1892   //
1893   fDiamondPointID = -1;
1894   if (addVertex) {
1895     fTrack->AddPoint(npto,&fDiamond);
1896     fMeasLoc  = fClusLoc.GetArray() + npto*3;
1897     fMeasGlo  = fClusGlo.GetArray() + npto*3;
1898     fSigmaLoc = fClusSigLoc.GetArray() + npto*3;
1899     fMeasLoc[0] = fMeasGlo[0] = fDiamond.GetX();
1900     fMeasLoc[1] = fMeasGlo[1] = fDiamond.GetY();
1901     fMeasLoc[2] = fMeasGlo[2] = fDiamond.GetZ();
1902     fSigmaLoc[0] = TMath::Sqrt(fDiamond.GetCov()[0]);
1903     fSigmaLoc[1] = TMath::Sqrt(fDiamond.GetCov()[3]);
1904     fSigmaLoc[2] = TMath::Sqrt(fDiamond.GetCov()[5]);
1905     fDiamondPointID = npto++;
1906   }
1907   //
1908   return fTrack;
1909 }
1910
1911 //________________________________________________________________________________________________________
1912 AliTrackPointArray *AliITSAlignMille2::SortTrack(const AliTrackPointArray *atp) 
1913 {
1914   /// sort alitrackpoints w.r.t. global Y direction
1915   AliTrackPointArray *atps=NULL;
1916   Int_t idx[20];
1917   Int_t npts=atp->GetNPoints();
1918   AliTrackPoint p;
1919   atps=new AliTrackPointArray(npts);
1920
1921   TMath::Sort(npts,atp->GetY(),idx);
1922
1923   for (int i=0; i<npts; i++) {
1924     atp->GetPoint(p,idx[i]);
1925     atps->AddPoint(i,&p);
1926     AliDebug(2,Form("Point[%d] = ( %+f , %+f , %+f )     volid = %d",i,atps->GetX()[i],atps->GetY()[i],atps->GetZ()[i],atps->GetVolumeID()[i] ));
1927   }
1928   return atps;
1929 }
1930
1931 //________________________________________________________________________________________________________
1932 Int_t AliITSAlignMille2::GetCurrentLayer() const 
1933 {
1934   // get current layer id
1935   if (!fGeoManager) {
1936     AliInfo("ITS geometry not initialized!");
1937     return -1;
1938   }
1939   return (Int_t)AliGeomManager::VolUIDToLayer(fCluster.GetVolumeID());
1940 }
1941
1942 //________________________________________________________________________________________________________
1943 Int_t AliITSAlignMille2::InitModuleParams() 
1944 {
1945   /// initialize geometry parameters for a given detector
1946   /// for current cluster (fCluster)
1947   /// fGlobalInitParam[] is set as:
1948   ///    [tx,ty,tz,psi,theta,phi]
1949   ///    (old was [tx,ty,tz,theta,psi,phi] ROOT's angles...)
1950   /// *** At the moment: using Raffalele's angles definition ***
1951   ///
1952   /// return 0 if success
1953   /// If module is found but has no parameters to vary, return 1
1954
1955   if (!fGeoManager) {
1956     AliInfo("ITS geometry not initialized!");
1957     return -1;
1958   }
1959
1960   // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster)
1961
1962   // set the internal index (index in module list)
1963   UShort_t voluid=fCluster.GetVolumeID();
1964   fCurrentSensID = AliITSAlignMille2Module::GetIndexFromVolumeID(voluid);
1965   //
1966   if (fCurrentSensID==-1) { // this is a special "vertex" module
1967     fCurrentModule = GetMilleModuleByVID(voluid);
1968     fCurrentSensID = fCurrentModule->GetIndex();
1969
1970   }
1971   else {
1972     //
1973     // IT IS VERY IMPORTANT: start from the end of the list, where the childs are located !!!
1974     Int_t k=fNModules-1;
1975     fCurrentModule = 0;
1976     // VERY IMPORTANT: if the sensors were explicitly provided, don't look in the supermodules  
1977     while (k>=0 && ! (fCurrentModule=GetMilleModule(k))->IsIn(voluid)) k--;
1978     if (k<0) return -3;
1979   }
1980   //
1981   for (int i=AliITSAlignMille2Module::kMaxParTot;i--;) fModuleInitParam[i] = 0.0;
1982   //
1983   int clID = fCluster.GetUniqueID()-1;
1984   if (clID<0) { // external cluster
1985     fMeasGlo  = &fExtClusterPar[0];
1986     fMeasLoc  = &fExtClusterPar[3];
1987     fSigmaLoc = &fExtClusterPar[6];
1988     fExtClusterPar[0] = fCluster.GetX();
1989     fExtClusterPar[1] = fCluster.GetY();
1990     fExtClusterPar[2] = fCluster.GetZ();
1991     //
1992     TGeoHMatrix *svMatrix = fCurrentModule->GetSensitiveVolumeMatrix(voluid);
1993     svMatrix->MasterToLocal(fMeasGlo,fMeasLoc);  
1994     TGeoHMatrix hcov;
1995     Double_t hcovel[9];
1996     hcovel[0]=double(fCluster.GetCov()[0]);
1997     hcovel[1]=double(fCluster.GetCov()[1]);
1998     hcovel[2]=double(fCluster.GetCov()[2]);
1999     hcovel[3]=double(fCluster.GetCov()[1]);
2000     hcovel[4]=double(fCluster.GetCov()[3]);
2001     hcovel[5]=double(fCluster.GetCov()[4]);
2002     hcovel[6]=double(fCluster.GetCov()[2]);
2003     hcovel[7]=double(fCluster.GetCov()[4]);
2004     hcovel[8]=double(fCluster.GetCov()[5]);
2005     hcov.SetRotation(hcovel);
2006     // now rotate in local system
2007     hcov.Multiply(svMatrix);
2008     hcov.MultiplyLeft(&svMatrix->Inverse());
2009     if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
2010     if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
2011     //
2012   }
2013   else {
2014     int offs = 3*clID;
2015     fMeasGlo  = fClusGlo.GetArray()  + offs;
2016     fMeasLoc  = fClusLoc.GetArray()  + offs;
2017     fSigmaLoc = fClusSigLoc.GetArray() + offs;
2018   }
2019   //
2020   // set minimum value for SigmaLoc to 10 micron 
2021   if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
2022   if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
2023   if (fCurrentSensID==kVtxSensID || fUseLocalYErr) if (fSigmaLoc[1]<0.0010) fSigmaLoc[1]=0.0010;
2024   //
2025   AliDebug(2,Form("Local coordinates of measured point : X=%+f  Y=%+f  Z=%+f \n",fMeasLoc[0] ,fMeasLoc[1] ,fMeasLoc[2] ));
2026   AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%g  fSigmaLocY=%g fSigmaLocZ=%g \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
2027   //   
2028   return 0;
2029 }
2030
2031 //________________________________________________________________________________________________________
2032 void AliITSAlignMille2::Print(Option_t*) const 
2033 {
2034   // print current status 
2035   printf("*** AliMillepede for ITS ***\n");
2036   printf("    Number of defined super modules: %d\n",fNModules);
2037   printf("    Obtained parameters refer to %s Deltas\n",fUseGlobalDelta ? "GLOBAL":"LOCAL");
2038   //
2039   if (fGeoManager)
2040     printf("    geometry loaded from %s\n",fGeometryPath.Data());
2041   else
2042     printf("    geometry not loaded\n");
2043   //  
2044   if (fUsePreAlignment) 
2045     printf("    using prealignment from %s \n",fPreDeltaPath.Data());
2046   else
2047     printf("    prealignment not used\n");    
2048   //
2049   //
2050   if (fBOn) 
2051     printf("    B Field set to %+f T - using helices\n",fBField);
2052   else
2053     printf("    B Field OFF - using straight lines \n");
2054   //
2055   if (fTPAFitter)
2056     printf("    Using AliITSTPArrayFit class for track fitting\n");
2057   else 
2058     printf("    Using StraightLine/Riemann fitter for track fitting\n");
2059   //
2060   printf("Using local Y error due to the sensor thickness: %s\n",(fUseLocalYErr && fTPAFitter) ? "ON":"OFF");
2061   //
2062   for (int itp=0;itp<kNDataType;itp++) {
2063     if (fRequirePoints[itp]) printf("    Required points in %s tracks:\n",itp==kCosmics? "cosmics" : "collisions");
2064     for (Int_t i=0; i<6; i++) {
2065       if (fNReqLayUp[itp][i]>0) printf("        Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[itp][i]);
2066       if (fNReqLayDown[itp][i]>0) printf("        Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[itp][i]);
2067       if (fNReqLay[itp][i]>0) printf("        Layer %d : %d points \n",i+1,fNReqLay[itp][i]);
2068     }
2069     for (Int_t i=0; i<3; i++) {
2070       if (fNReqDetUp[itp][i]>0) printf("        Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[itp][i]);
2071       if (fNReqDetDown[itp][i]>0) printf("        Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[itp][i]);
2072       if (fNReqDet[itp][i]>0) printf("        Detector %d : %d points \n",i+1,fNReqDet[itp][i]);
2073     }
2074   }
2075   printf("        SDD VDrift correction         : %s",fIsSDDVDriftMult ? "Mult":"Add");
2076   printf("        Weight acc. to pT in power    : %f",fWeightPt);
2077   //  
2078   printf("\n    Millepede configuration parameters:\n");
2079   printf("        init factor for chi2 cut      : %.4f\n",fStartFac);
2080   printf("        final factor for chi2 cut     : %.4f\n",fFinalFac);
2081   printf("        first iteration cut value     : %.4f\n",fResCutInitial);
2082   printf("        other iterations cut value    : %.4f\n",fResCut);
2083   printf("        number of stddev for chi2 cut : %d\n",fNStdDev);
2084   printf("        def.scaling for local sigmas  : %.4f %.4f %.4f\n",fSigmaFactor[0],fSigmaFactor[1],fSigmaFactor[2]);
2085   printf("        min.tracks per module         : %d\n",fMinPntPerSens);
2086   //
2087   printf("List of defined modules:\n");
2088   printf("  intidx\tindex\tvoluid\tname\n");
2089   for (int i=0; i<fNModules; i++) {
2090     AliITSAlignMille2Module* md = GetMilleModule(i); 
2091     printf("  %d\t%d\t%d\t%s\n",i,md->GetIndex(),md->GetVolumeID(),md->GetName());
2092   }
2093 }
2094
2095 //________________________________________________________________________________________________________
2096 AliITSAlignMille2Module  *AliITSAlignMille2::GetMilleModuleByVID(UShort_t voluid) const
2097 {
2098   // return pointer to a defined supermodule
2099   // return NULL if error
2100   Int_t i=IsVIDDefined(voluid);
2101   if (i<0) return NULL;
2102   return GetMilleModule(i);
2103 }
2104
2105 //________________________________________________________________________________________________________
2106 AliITSAlignMille2Module  *AliITSAlignMille2::GetMilleModuleBySymName(const Char_t* symname) const
2107 {
2108   // return pointer to a defined supermodule
2109   // return NULL if error
2110   UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
2111   if (vid>0) return GetMilleModuleByVID(vid);
2112   else {    // this is not alignable module, need to look within defined supermodules
2113     int i = IsSymDefined(symname);
2114     if (i>=0) return  GetMilleModule(i);
2115   }
2116   return 0;
2117 }
2118
2119 //________________________________________________________________________________________________________
2120 AliITSAlignMille2Module  *AliITSAlignMille2::GetMilleModuleIfContained(const Char_t* symname) const
2121 {
2122   // return pointer to a defined/contained supermodule
2123   // return NULL otherwise
2124   int i = IsSymContained(symname);
2125   return i<0 ? 0 : GetMilleModule(i);
2126 }
2127
2128 //________________________________________________________________________________________________________
2129 AliAlignObjParams* AliITSAlignMille2::GetPrealignedObject(const Char_t* symname) const
2130 {
2131   // get delta from prealignment for given volume
2132   if (!fPrealignment) return 0;
2133   for (int ipre=fPrealignment->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
2134     AliAlignObjParams* preob = (AliAlignObjParams*)fPrealignment->At(ipre);
2135     if (!strcmp(preob->GetSymName(),symname)) return preob;
2136   }
2137   return 0;
2138 }
2139
2140 //________________________________________________________________________________________________________
2141 AliAlignObjParams* AliITSAlignMille2::GetConstrRefObject(const Char_t* symname) const
2142 {
2143   // get delta with respect to which the constraint is declared
2144   if (!fConstrRef) return 0;
2145   for (int ipre=fConstrRef->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
2146     AliAlignObjParams* preob = (AliAlignObjParams*)fConstrRef->At(ipre);
2147     if (!strcmp(preob->GetSymName(),symname)) return preob;
2148   }
2149   return 0;
2150 }
2151
2152 //________________________________________________________________________________________________________
2153 Bool_t AliITSAlignMille2::InitRiemanFit() 
2154 {
2155   // Initialize Riemann Fitter for current track
2156   // return kFALSE if error
2157
2158   if (!fBOn) return kFALSE;
2159
2160   Int_t npts=0;
2161   AliTrackPoint ap;
2162   npts = fTrack->GetNPoints();
2163   AliDebug(3,Form("Fitting track with %d points",npts));
2164   if (!fRieman) fRieman = new AliTrackFitterRieman();
2165   fRieman->Reset();
2166   fRieman->SetTrackPointArray(fTrack);
2167
2168   TArrayI ai(npts);
2169   for (Int_t ipt=0; ipt<npts; ipt++) ai[ipt]=fTrack->GetVolumeID()[ipt];
2170   
2171   // fit track with 5 params in his own tracking-rotated reference system
2172   // xc = -p[1]/p[0];
2173   // yc = 1/p[0];
2174   // R  = sqrt( x0*x0 + y0*y0 - y0*p[2]);
2175   if (!fRieman->Fit(&ai,NULL,(AliGeomManager::ELayerID)1,(AliGeomManager::ELayerID)6)) {
2176     return kFALSE;
2177   }
2178
2179   for (int i=0; i<5; i++)
2180     fLocalInitParam[i] = fRieman->GetParam()[i];
2181   
2182   return kTRUE;
2183 }
2184
2185 //________________________________________________________________________________________________________
2186 void trackFit2D(Int_t &, Double_t *, double &chi2, double *par, int flag)
2187 {
2188   // local function for minuit
2189   const double kTiny = 1.e-14;
2190   chi2 = 0;
2191   static AliTrackPoint pnt;
2192   static Bool_t fullErr2D;
2193   //
2194   if (flag==1) fullErr2D = kFALSE;//kTRUE;
2195   //  fullErr2D = kTRUE;
2196   enum {kAX,kAZ,kBX,kBZ};
2197   enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5};
2198   //
2199   AliITSAlignMille2* alig = AliITSAlignMille2::GetInstance();
2200   AliTrackPointArray* track = alig->GetCurrentTrack();
2201   //
2202   int npts = track->GetNPoints();
2203   for (int ip=0;ip<npts;ip++) {
2204     track->GetPoint(pnt,ip);
2205     const float *cov = pnt.GetCov();
2206     double y  = pnt.GetY();
2207     double dx = pnt.GetX() - (par[kAX]+y*par[kBX]);
2208     double dz = pnt.GetZ() - (par[kAZ]+y*par[kBZ]);
2209     double xxe = cov[kXX];
2210     double zze = cov[kZZ];
2211     double xze = cov[kXZ];
2212     //
2213     if (fullErr2D) {
2214       xxe += par[kBX]*par[kBX]*cov[kYY]-2.*par[kBX]*cov[kXY];
2215       zze += par[kBZ]*par[kBZ]*cov[kYY]-2.*par[kBZ]*cov[kZY];
2216       xze += par[kBX]*par[kBZ]*cov[kYY]-cov[kYZ]*par[kBZ]-cov[kXY]*par[kBX];
2217     }
2218     //
2219     double det = xxe*zze - xze*xze;
2220     if (det<kTiny) {
2221       printf("Negative diag. error (det=%+e) |sxx:%+e szz:%+e sxz:%+e| bx:%+e bz:%+e|\n"
2222              "Discarding correlation term\n",det,xxe,zze,xze,par[kBX],par[kBZ]);
2223       xxe = cov[kXX];
2224       zze = cov[kZZ];
2225       xze = cov[kXZ];
2226       fullErr2D = kFALSE;
2227     }
2228     double xxeI = zze/det;
2229     double zzeI = xxe/det;
2230     double xzeI =-xze/det;
2231     //
2232     chi2 += dx*dx*xxeI + dz*dz*zzeI + 2.*dx*dz*xzeI;
2233     // 
2234     //    printf("%d | %+e %+e %+e %+e %+e -> %+e\n",ip,dx,dz,xxeI,zzeI,xzeI,  chi2);
2235   }
2236   //
2237 }
2238
2239 //________________________________________________________________________________________________________
2240 void AliITSAlignMille2::InitTrackParams(int meth) 
2241 {
2242   /// initialize local parameters with different methods
2243   /// for current track (fTrack)
2244   Int_t npts=0;
2245   AliTrackPoint ap;
2246   double sX=0,sXY=0,sZ=0,sZY=0,sY=0,sYY=0,det=0;
2247   // simple linear interpolation
2248   // get local starting parameters (to be substituted by ESD track parms)
2249   // local parms (fLocalInitParam[]) are:
2250   //      [0] = global x coord. of straight line intersection at y=0 plane
2251   //      [1] = global z coord. of straight line intersection at y=0 plane
2252   //      [2] = px/py  
2253   //      [3] = pz/py
2254   // test #1: linear fit in x(y) and z(y)
2255   npts = fTrack->GetNPoints();
2256   AliDebug(3,Form("*** initializing track with %d points ***",npts));
2257   for (int i=npts;i--;) {
2258     sY  += fTrack->GetY()[i];
2259     sYY += fTrack->GetY()[i]*fTrack->GetY()[i];
2260     sX  += fTrack->GetX()[i];
2261     sXY += fTrack->GetX()[i]*fTrack->GetY()[i];
2262     sZ  += fTrack->GetZ()[i];
2263     sZY += fTrack->GetZ()[i]*fTrack->GetY()[i];
2264   }
2265   det = sYY*npts-sY*sY;
2266   if (IsZero(det)) det = 1E-16;
2267   fLocalInitParam[0] = (sX*sYY-sY*sXY)/det;
2268   fLocalInitParam[2] = (sXY*npts-sY*sX)/det;
2269   //
2270   fLocalInitParam[1] = (sZ*sYY-sY*sZY)/det;
2271   fLocalInitParam[3] = (sZY*npts-sY*sZ)/det;
2272   // pepo200709
2273   fLocalInitParam[4] = 0.0;
2274   // endpepo200709
2275
2276   AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %+f    ugx = %+f\n",fLocalInitParam[0],fLocalInitParam[2]));
2277   //
2278   if (meth==1) return;
2279   //
2280   // perform full fit accounting for cov.matrix
2281   static TVirtualFitter *minuit = 0;
2282   static Double_t step[5]   = {1E-3,1E-3,1E-4,1E-4,1E-5};
2283   static Double_t arglist[10];
2284   //
2285   if (!minuit) {
2286     minuit = TVirtualFitter::Fitter(0,4);
2287     minuit->SetFCN(trackFit2D);
2288     arglist[0] = 1;
2289     minuit->ExecuteCommand("SET ERR",arglist, 1);
2290     //
2291     arglist[0] = -1;
2292     minuit->ExecuteCommand("SET PRINT",arglist,1);
2293     //
2294   }
2295   //
2296   minuit->SetParameter(0, "ax",   fLocalInitParam[0], step[0], 0,0);
2297   minuit->SetParameter(1, "az",   fLocalInitParam[1], step[1], 0,0);
2298   minuit->SetParameter(2, "bx",   fLocalInitParam[2], step[2], 0,0);
2299   minuit->SetParameter(3, "bz",   fLocalInitParam[3], step[3], 0,0);
2300   //
2301   arglist[0] = 1000; // number of function calls 
2302   arglist[1] = 0.001; // tolerance 
2303   minuit->ExecuteCommand("MIGRAD",arglist,2);
2304   //
2305   for (int i=0;i<4;i++) fLocalInitParam[i] = minuit->GetParameter(i);
2306   for (int i=0;i<4;i++) for (int j=0;j<4;j++) fLocalInitParEr[i][j] = minuit->GetCovarianceMatrixElement(i,j);
2307   /*
2308   double amin,edm,errdef;
2309   int nvpar,nparx;
2310   minuit->GetStats(amin,edm,errdef,nvpar,nparx);
2311   amin /= (2*npts - 4);
2312   printf("Mchi2: %+e\n",amin);
2313   */
2314   //
2315 }
2316
2317 //________________________________________________________________________________________________________
2318 Int_t AliITSAlignMille2::IsSymDefined(const Char_t* symname) const
2319 {
2320   // checks if supermodule with this symname is defined and return the internal index
2321   // return -1 if not.
2322   for (int k=fNModules;k--;) if (!strcmp(symname,GetMilleModule(k)->GetName())) return k;
2323   return -1; 
2324 }
2325
2326 //________________________________________________________________________________________________________
2327 Int_t AliITSAlignMille2::IsSymContained(const Char_t* symname) const
2328 {
2329   // checks if module with this symname is defined and return the internal index
2330   // return -1 if not.
2331   UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
2332   if (vid>0) return IsVIDContained(vid);
2333   // only sensors have real vid, but maybe we have a supermodule with fake vid? 
2334   // IMPORTANT: always start from the end to start from the sensors
2335   return IsSymDefined(symname);
2336 }
2337
2338 //________________________________________________________________________________________________________
2339 Int_t AliITSAlignMille2::IsVIDDefined(UShort_t voluid) const
2340 {
2341   // checks if supermodule 'voluid' is defined and return the internal index
2342   // return -1 if not.
2343   for (int k=fNModules;k--;) if (voluid==GetMilleModule(k)->GetVolumeID()) return k;
2344   return -1; 
2345 }
2346
2347 //________________________________________________________________________________________________________
2348 Int_t AliITSAlignMille2::IsVIDContained(UShort_t voluid) const
2349 {
2350   // checks if the sensitive module 'voluid' is contained inside a supermodule 
2351   // and return the internal index of the last identified supermodule
2352   // return -1 if error
2353   // IMPORTANT: always start from the end to start from the sensors
2354   if (AliITSAlignMille2Module::GetIndexFromVolumeID(voluid)<0) return -1;
2355   for (int k=fNModules;k--;) if (GetMilleModule(k)->IsIn(voluid)) return k;
2356   return -1; 
2357 }
2358
2359 //________________________________________________________________________________________________________
2360 Int_t AliITSAlignMille2::GetRequestedModID(UShort_t voluid) const
2361 {
2362   // checks if the sensitive module 'voluid' is contained inside a supermodule 
2363   // and return the internal index of the last identified supermodule
2364   // return -1 if error
2365   // IMPORTANT: always start from the end to start from the sensors
2366   if (AliITSAlignMille2Module::GetIndexFromVolumeID(voluid)<0) return -1;
2367   int k;
2368   for (k=fNModules;k--;) if (GetMilleModule(k)->IsIn(voluid)) break;
2369   if (k<0) return -1;
2370   AliITSAlignMille2Module* md = GetMilleModule(k);
2371   while (md && md->IsNotInConf()) md = md->GetParent();
2372   if (md) return int(md->GetUniqueID());
2373   else return -1; 
2374 }
2375
2376 //________________________________________________________________________________________________________
2377 Int_t AliITSAlignMille2::CheckCurrentTrack() 
2378 {
2379   /// checks if AliTrackPoints belongs to defined modules
2380   /// return number of good poins
2381   /// return 0 if not enough points
2382
2383   Int_t npts = fTrack->GetNPoints();
2384   Int_t ngoodpts=0;
2385   // debug points
2386   for (int j=0; j<npts; j++) if (IsVIDContained(fTrack->GetVolumeID()[j])>=0) ngoodpts++;
2387   //
2388   if (ngoodpts<fMinNPtsPerTrack) return 0;
2389
2390   return ngoodpts;
2391 }
2392
2393 //________________________________________________________________________________________________________
2394 Int_t AliITSAlignMille2::ProcessTrack(const AliTrackPointArray *track, Double_t wgh) 
2395 {
2396   /// Process track; Loop over hits and set local equations
2397   /// here 'track' is a AliTrackPointArray
2398   /// return 0 if success;
2399   //
2400   if (!fIsMilleInit) Init();
2401   //
2402   Int_t npts = track->GetNPoints();
2403   AliDebug(2,Form("*** Input track with %d points ***",npts));
2404
2405   // preprocessing of the input track: keep only points in defined volumes,
2406   // move points if prealignment is set, sort by Yglo if required
2407   fTrackWeight = wgh;
2408   fTrack=PrepareTrack(track);
2409   if (!fTrack) {
2410     RemoveHelixFitConstraint();
2411     RemoveVertexConstraint();
2412     return -1;
2413   }
2414   npts = fTrack->GetNPoints();
2415   if (npts>kMaxPoints) {
2416     AliError(Form("Compiled with kMaxPoints=%d, current track has %d points",kMaxPoints,npts));
2417   }
2418   AliDebug(2,Form("*** Processing prepared track with %d points ***",npts));
2419   //
2420   npts = FitTrack();
2421   if (npts<0) return npts;
2422   //
2423   //  printf("Params: "); for (int i=0;i<fNLocal;i++) printf("%+.2e ",fLocalInitParam[i]); printf("\n");//RRR
2424   Int_t nloceq=0;
2425   Int_t ngloeq=0;
2426   static Mille2Data md[kMaxPoints];
2427   //
2428   for (Int_t ipt=0; ipt<npts; ipt++) {
2429     fTrack->GetPoint(fCluster,ipt);
2430     fCluster.SetUniqueID(ipt+1);
2431     AliDebug(2,Form("\n--- processing point %d --- \n",ipt));    
2432
2433     // set geometry parameters for the the current module
2434     if (InitModuleParams()) continue;
2435     AliDebug(2,Form("    VolID=%d  Index=%d  InternalIdx=%d  symname=%s\n", 
2436                     track->GetVolumeID()[ipt], fCurrentModule->GetIndex(),
2437                     fCurrentModule->GetUniqueID(), AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
2438     AliDebug(2,Form("    Preprocessed Point = ( %+f , %+f , %+f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
2439     int res = fTPAFitter ? AddLocalEquationTPA(md[nloceq]) : AddLocalEquation(md[nloceq]);
2440     if (res<0) {fTotBadLocEqPoints++; nloceq = 0; break;}
2441     else if (res==0) nloceq++;
2442     else {nloceq++; ngloeq++;}
2443   } // end loop over points
2444   //
2445   fTrack=NULL;
2446   // not enough good points?
2447   if (nloceq<fMinNPtsPerTrack || ngloeq<1) return -1;
2448   //
2449   // finally send local equations to millepede
2450   SetLocalEquations(md,nloceq);
2451   fMillepede->SaveRecordData(); // RRR
2452   fCurvFitWasConstrained = kFALSE; // restore default
2453   //
2454   return 0;
2455 }
2456
2457 //________________________________________________________________________________________________________
2458 Int_t AliITSAlignMille2::FitTrack() 
2459 {
2460   // Fit the track with selected constraints
2461   //
2462   const Double_t kfDiamondTolerance = 0.1;  //diamond tolerance on top of the MS error
2463   if (!fTrack) return -1;
2464   int npts = fTrack->GetNPoints();
2465   //
2466   if (fTPAFitter) {  // use dediacted fitter
2467     //
2468     // if the diamond point is attached, for the moment don't include it in the fit
2469     fTPAFitter->AttachPoints(fTrack,0, npts-1); 
2470     fTPAFitter->SetBz(fBField);
2471     fTPAFitter->SetTypeCosmics(IsTypeCosmics());
2472     if (fIniTrackParamsMeth==1) fTPAFitter->SetIgnoreCov();
2473     //
2474     double chi2;
2475     double chi2f = 0;
2476     double dca2err;
2477     double dca2 = 0.;
2478     Bool_t fitIsDone = kFALSE;
2479     if (fUseDiamond && fDiamondPointID>0 && fCheckDiamondPoint==kDiamondCheckIfPrompt) { // diamond constraint was added, check if the track looks like prompt
2480       fTPAFitter->SetFirstLast(0,fDiamondPointID-1);
2481       if (IsCovIScaleTouched()) for (int i=npts;i--;) fTPAFitter->SetCovIScale(i,GetCovIScale(i));
2482       //
2483       chi2f = fTPAFitter->Fit(fConstrCharge,fConstrPT,fConstrPTErr);
2484       if ( chi2f<0 || (chi2f>fNStdDev*fStartFac && fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations()) ) { //RRR
2485         AliInfo(Form("Track fit failed on checking if it is prompt! skipping this track... Chi2:%+e",chi2f));
2486         fTPAFitter->Reset();
2487         //      fTrack = NULL;
2488         return -5;
2489       }
2490       double xyzRes[3];
2491       fTPAFitter->GetResiduals(xyzRes,&fDiamondI,kTRUE);
2492       dca2 = xyzRes[0]*xyzRes[0] + xyzRes[1]*xyzRes[1];
2493       double pT = IsFieldON() ? fTPAFitter->GetPt() : 0.45;
2494       if (pT<0.1) pT = 0.1;
2495       dca2err = kfDiamondTolerance + 0.02/pT;
2496       if (dca2>dca2err*dca2err) { // this is secondary
2497         int* clst = (int*) fTrack->GetClusterType();
2498         clst[fDiamondPointID] = -1;;
2499         fDiamondPointID = -1; 
2500         fitIsDone = kTRUE;
2501         npts--;
2502       }
2503       else fTPAFitter->SetFirstLast(0,fDiamondPointID); // fit with diamond
2504     }
2505     //    fTPAFitter->SetParAxis(1);
2506     if (!fitIsDone) {
2507       if (IsCovIScaleTouched()) for (int i=npts;i--;) fTPAFitter->SetCovIScale(i,GetCovIScale(i));
2508       chi2 = fTPAFitter->Fit(fConstrCharge,fConstrPT,fConstrPTErr);
2509     }
2510     //
2511     RemoveHelixFitConstraint();  // suppress eventual constraints to not affect fit of the next track
2512     RemoveVertexConstraint(); // same ...
2513     //
2514     if ( !fitIsDone && (chi2<0 || (chi2>fNStdDev*fStartFac && fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations())) ) { //RRR
2515       AliInfo(Form("Track fit failed! skipping this track... Chi2:%+e",chi2));
2516       if (fUseDiamond && fDiamondPointID>0 && fCheckDiamondPoint==kDiamondCheckIfPrompt) AliInfo(Form("VertexFree fit gave Chi2:%+e with residual %+e",chi2f,TMath::Sqrt(dca2)));
2517       /*
2518         fTrack->Print("");
2519         fTPAFitter->FitHelixCrude();
2520         fTPAFitter->SetFitDone();
2521         fTPAFitter->Print();
2522       */
2523       fTPAFitter->Reset();
2524       //      fTrack = NULL;
2525       return -5;
2526     }
2527     fNLocal = fTPAFitter->IsFieldON() ? 5:4; // Attention: the fitter might have decided to work in line mode
2528     npts  = fTPAFitter->GetLast() - fTPAFitter->GetFirst() + 1; // actual number of points
2529     /*
2530       double *pr = fTPAFitter->GetParams();
2531       printf("FtPar: %+.5e  %+.5e  %+.5e  %+.5e | chi2:%.3e\n",pr[2],pr[0],pr[3],pr[1],chi2); // RRR
2532     */
2533   }
2534   else {
2535     //
2536     if (!fBOn) { // straight lines  
2537       // set local starting parameters (to be substituted by ESD track parms)
2538       // local parms (fLocalInitParam[]) are:
2539       //      [0] = global x coord. of straight line intersection at y=0 plane
2540       //      [1] = global z coord. of straight line intersection at y=0 plane
2541       //      [2] = px/py  
2542       //      [3] = pz/py
2543       InitTrackParams(fIniTrackParamsMeth); 
2544       /*
2545       double *pr = fLocalInitParam;
2546       printf("FtPar: %+.5e  %+.5e  %+.5e  %+.5e |\n",pr[0],pr[1],pr[2],pr[3]); // RRR
2547       */
2548     } 
2549     else {
2550       // local parms (fLocalInitParam[]) are the Riemann Fitter params
2551       if (!InitRiemanFit()) {
2552         AliInfo("Riemann fit failed! skipping this track...");
2553         fTrack=NULL;
2554         return -5;
2555       }
2556     }
2557   }
2558   return npts;
2559   //
2560 }
2561
2562 //________________________________________________________________________________________________________
2563 Int_t AliITSAlignMille2::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar) 
2564 {
2565   /// calculate track intersection point in local coordinates
2566   /// according with a given set of parameters (local(4) and global(6))
2567   /// and fill fPintLoc/Glo
2568   ///    local are:   pgx0, pgz0, ugx, ugz   OR   riemann fitters pars
2569   ///    global are:  tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
2570   /// return 0 if success
2571   
2572   AliDebug(3,Form("lpar = %g %g %g %g %g\ngpar= %g %g %g %g %g %g\n",lpar[0],lpar[1],lpar[2],lpar[3],lpar[4],gpar[0],gpar[1],gpar[2],gpar[3],gpar[4],gpar[5]));
2573   AliDebug(3,Form("deltalpar = %g %g %g %g %g\n",lpar[0]-fLocalInitParam[0],lpar[1]-fLocalInitParam[1],lpar[2]-fLocalInitParam[2],lpar[3]-fLocalInitParam[3],lpar[4]-fLocalInitParam[4]));
2574
2575   
2576   // prepare the TGeoHMatrix
2577   TGeoHMatrix *tempHMat = fCurrentModule->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar,
2578                                                                            !fUseGlobalDelta);
2579   if (!tempHMat) return -1;
2580   
2581   Double_t v0g[3]; // vector with straight line direction in global coord.
2582   Double_t p0g[3]; // point of the straight line (glo)
2583   
2584   if (fBOn) { // B FIELD!
2585     AliTrackPoint prf; 
2586     for (int ip=0; ip<5; ip++)
2587       fRieman->SetParam(ip,lpar[ip]);
2588
2589     if (!fRieman->GetPCA(fCluster,prf))  {
2590       AliInfo(Form("error in GetPCA for point %d",fCluster.GetVolumeID()));
2591       return -3;
2592     }
2593     // now determine straight line passing tangent to fit curve at prf
2594     // ugx = dX/dY_glo = DeltaX/DeltaY_glo
2595     // mo' P1=(X,Y,Z)_glo_prf
2596     //       => (x,y,Z)_trk_prf ruotando di alpha...
2597     Double_t alpha=fRieman->GetAlpha();
2598     Double_t x1g = prf.GetX();
2599     Double_t y1g = prf.GetY();
2600     Double_t z1g = prf.GetZ();
2601     Double_t x1t =  x1g*TMath::Cos(alpha) + y1g*TMath::Sin(alpha);
2602     Double_t y1t = -x1g*TMath::Sin(alpha) + y1g*TMath::Cos(alpha);
2603     Double_t z1t =  z1g;    
2604
2605     Double_t x2t = x1t+1.0;
2606     Double_t y2t = y1t+fRieman->GetDYat(x1t);
2607     Double_t z2t = z1t+fRieman->GetDZat(x1t);
2608     Double_t x2g =  x2t*TMath::Cos(alpha) - y2t*TMath::Sin(alpha);
2609     Double_t y2g =  x2t*TMath::Sin(alpha) + y2t*TMath::Cos(alpha);
2610     Double_t z2g =  z2t;  
2611
2612     AliDebug(3,Form("Riemann frame:  fAlpha = %+f  =  %+f  ",alpha,alpha*180./TMath::Pi()));
2613     AliDebug(3,Form("   prf_glo=( %+f , %+f , %+f )  prf_rf=( %+f , %+f , %+f )\n", x1g,y1g,z1g, x1t,y1t,z1t));
2614     AliDebug(3,Form("   mov_glo=( %+f , %+f , %+f )      rf=( %+f , %+f , %+f )\n",x2g,y2g,z2g, x2t,y2t,z2t));
2615         
2616     if (TMath::Abs(y2g-y1g)<1e-15) {
2617       AliInfo("DeltaY=0! Cannot proceed...");
2618       return -1;
2619     }
2620     // ugx,1,ugz
2621     v0g[0] = (x2g-x1g)/(y2g-y1g);
2622     v0g[1]=1.0;
2623     v0g[2] = (z2g-z1g)/(y2g-y1g);
2624     
2625     // point: just keep prf
2626     p0g[0]=x1g;
2627     p0g[1]=y1g;
2628     p0g[2]=z1g;
2629   }  
2630   else { // staight line
2631     // vector of initial straight line direction in glob. coord
2632     v0g[0]=lpar[2];
2633     v0g[1]=1.0;
2634     v0g[2]=lpar[3];
2635     
2636     // intercept in yg=0 plane in glob coord
2637     p0g[0]=lpar[0];
2638     p0g[1]=0.0;
2639     p0g[2]=lpar[1];
2640   }
2641   AliDebug(3,Form("Line vector: ( %+f , %+f , %+f )  point:( %+f , %+f , %+f )\n",v0g[0],v0g[1],v0g[2],p0g[0],p0g[1],p0g[2]));
2642   
2643   // same in local coord.
2644   Double_t p0l[3],v0l[3];
2645   tempHMat->MasterToLocalVect(v0g,v0l);
2646   tempHMat->MasterToLocal(p0g,p0l);
2647   
2648   if (TMath::Abs(v0l[1])<1e-15) {
2649     AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
2650     return -1;
2651   }
2652   
2653   // local intersection point
2654   fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
2655   fPintLoc[1] = 0;
2656   fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
2657   
2658   // global intersection point
2659   tempHMat->LocalToMaster(fPintLoc,fPintGlo);
2660   AliDebug(3,Form("Intesect. point: L( %+f , %+f , %+f )  G( %+f , %+f , %+f )\n",fPintLoc[0],fPintLoc[1],fPintLoc[2],fPintGlo[0],fPintGlo[1],fPintGlo[2]));
2661   
2662   return 0;
2663 }
2664
2665 //________________________________________________________________________________________________________
2666 Int_t AliITSAlignMille2::CalcDerivatives(Int_t paridx, Bool_t islpar) 
2667 {
2668   /// calculate numerically (ROOT's style) the derivatives for
2669   /// local X intersection  and local Z intersection
2670   /// parlist: local  (islpar=kTRUE)  pgx0, pgz0, ugx0, ugz0  OR riemann's params
2671   ///          global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
2672   /// return 0 if success
2673   
2674   // copy initial parameters
2675   Double_t lpar[kNLocal];
2676   Double_t gpar[kNParCh];
2677   Double_t *derivative;
2678   for (Int_t i=0; i<kNLocal; i++) lpar[i]=fLocalInitParam[i];
2679   for (Int_t i=0; i<kNParCh; i++) gpar[i]=fModuleInitParam[i];
2680
2681   // trial with fixed dpar...
2682   Double_t dpar = 0.;
2683
2684   if (islpar) { // track parameters
2685     //dpar=fLocalInitParam[paridx]*0.001;
2686     // set minimum dpar
2687     derivative = fDerivativeLoc[paridx];
2688     if (!fBOn) {
2689       if (paridx<3) dpar=1.0e-4; // translations
2690       else dpar=1.0e-6; // direction
2691     }
2692     else { // B Field
2693       // pepo: proviamo con 1/1000, poi evenctually 1/100...
2694       Double_t dfrac=0.01;
2695       switch(paridx) {
2696       case 0:
2697         // RMS cosmics: 1e-4
2698         dpar = TMath::Max(1.0e-6,TMath::Abs(fLocalInitParam[paridx]*dfrac)); 
2699         break;
2700       case 1: 
2701         // RMS cosmics: 0.2
2702         dpar = TMath::Max(0.002,TMath::Abs(fLocalInitParam[paridx]*dfrac)); 
2703         break;
2704       case 2: 
2705         // RMS cosmics: 9
2706         dpar = TMath::Max(0.09,TMath::Abs(fLocalInitParam[paridx]*dfrac)); 
2707         break;
2708       case 3: 
2709         // RMS cosmics: 7
2710         dpar = TMath::Max(0.07,TMath::Abs(fLocalInitParam[paridx]*dfrac)); 
2711         break;
2712       case 4: 
2713         // RMS cosmics: 0.3
2714         dpar = TMath::Max(0.003,TMath::Abs(fLocalInitParam[paridx]*dfrac)); 
2715         break;
2716       }
2717     }
2718   }
2719   else { // alignment global parameters
2720     derivative = fDerivativeGlo[paridx];
2721     //dpar=fModuleInitParam[paridx]*0.001;
2722     if (paridx<3) dpar=1.0e-4; // translations
2723     else dpar=1.0e-2; // angles    
2724   }
2725
2726   AliDebug(3,Form("+++ using dpar=%g",dpar));
2727   
2728   // calculate derivative ROOT's like:
2729   //  using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
2730   Double_t pintl1[3]; // f(x-h)
2731   Double_t pintl2[3]; // f(x-h/2)
2732   Double_t pintl3[3]; // f(x+h/2)
2733   Double_t pintl4[3]; // f(x+h)
2734     
2735   // first values
2736   if (islpar) lpar[paridx] -= dpar;
2737   else gpar[paridx] -= dpar;
2738   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2739   for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
2740
2741   // second values
2742   if (islpar) lpar[paridx] += dpar/2;
2743   else gpar[paridx] += dpar/2;
2744   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2745   for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
2746
2747   // third values
2748   if (islpar) lpar[paridx] += dpar;
2749   else gpar[paridx] += dpar;
2750   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2751   for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
2752
2753   // fourth values
2754   if (islpar) lpar[paridx] += dpar/2;
2755   else gpar[paridx] += dpar/2;
2756   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2757   for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
2758
2759   Double_t h2 = 1./(2.*dpar);
2760   Double_t d0 = pintl4[0]-pintl1[0];
2761   Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
2762   derivative[0] = h2*(4*d2 - d0)/3.;
2763   if (TMath::Abs(derivative[0]) < 1.0e-9) derivative[0] = 0.0;
2764
2765   d0 = pintl4[2]-pintl1[2];
2766   d2 = 2.*(pintl3[2]-pintl2[2]);
2767   derivative[2] = h2*(4*d2 - d0)/3.;
2768   if (TMath::Abs(derivative[2]) < 1.0e-9) derivative[2]=0.0;
2769
2770   AliDebug(3,Form("\n+++ derivatives +++ \n"));
2771   AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",derivative[0]));
2772   AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",derivative[2]));
2773   
2774   return 0;
2775 }
2776
2777 //________________________________________________________________________________________________________
2778 Int_t AliITSAlignMille2::AddLocalEquation(Mille2Data &m) 
2779 {
2780   /// Define local equation for current cluster in X and Z coor.
2781   /// and store them to memory
2782   /// return -1 in case of failure to build some equation
2783   ///         0 if no free global parameters were found but local eq is built
2784   ///         1 if both local and global eqs are built
2785   //
2786   // store first intersection point
2787   if (CalcIntersectionPoint(fLocalInitParam, fModuleInitParam)) return -1;  
2788   for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
2789
2790   AliDebug(2,Form("Intersect. point: L( %+f , %+f , %+f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
2791   
2792   // calculate local derivatives numerically
2793   Bool_t zeroX = kTRUE;
2794   Bool_t zeroZ = kTRUE;
2795   //
2796   for (Int_t i=0; i<fNLocal; i++) {
2797     if (CalcDerivatives(i,kTRUE)) return -1;
2798     m.fDerLoc[i][kX] = fDerivativeLoc[i][0];
2799     m.fDerLoc[i][kZ] = fDerivativeLoc[i][2];
2800     if (zeroX) zeroX = IsZero(fDerivativeLoc[i][0]);
2801     if (zeroZ) zeroZ = IsZero(fDerivativeLoc[i][2]);
2802   }
2803   //  for (Int_t i=0; i<fNLocal; i++) AliDebug(2,Form("Local parameter %d - dXdpar = %g  - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
2804   //
2805   if (zeroX) {AliInfo("Skipping: zero local X derivatives!"); return -1;}
2806   if (zeroZ) {AliInfo("Skipping: zero local Z derivatives!"); return -1;}
2807   //
2808   int status = 0;
2809   int ifill = 0;
2810   //
2811   AliITSAlignMille2Module* endModule = fCurrentModule;
2812   //
2813   zeroX = zeroZ = kTRUE;
2814   Bool_t dfDone[kNParCh];
2815   for (int i=kNParCh;i--;) dfDone[i] = kFALSE;
2816   m.fNModFilled = 0;
2817   // 
2818   // special block for SDD derivatives
2819   Double_t jacobian[kNParChGeom];
2820   Int_t nmodTested = 0;
2821   //
2822   do {
2823     if (fCurrentModule->GetNParFree()==0) continue;
2824     nmodTested++;
2825     for (Int_t i=0; i<kNParChGeom; i++) {   // common for all sensors: derivatives over geom params 
2826       //
2827       if (!fUseGlobalDelta) dfDone[i] = kFALSE; // for global deltas the derivatives at diff. levels are different
2828       if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2829       if (!dfDone[i]) { 
2830         if (CalcDerivatives(i,kFALSE)) return -1; 
2831         else {
2832           dfDone[i] = kTRUE;
2833           if (zeroX) zeroX = IsZero(fDerivativeGlo[i][0]);
2834           if (zeroZ) zeroZ = IsZero(fDerivativeGlo[i][2]);
2835         }
2836       }
2837       //
2838       m.fDerGlo[ifill][kX] = fDerivativeGlo[i][0];
2839       m.fDerGlo[ifill][kZ] = fDerivativeGlo[i][2];
2840       m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
2841     }
2842     //
2843     // specific for special sensors
2844     Int_t sddLR = -1;
2845     if ( fCurrentModule->IsSDD() && 
2846          (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0  ||
2847           //      fCurrentModule->GetParOffset(sddLR = fMeasLoc[kX]>0 ?
2848           fCurrentModule->GetParOffset(sddLR = GetVDriftSDD()>0 ? 
2849                                        AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR)>=0)
2850          ) {
2851       //
2852       // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
2853       // where V0 and T are the nominal drift velocity, time and time0
2854       // and the dT0 and dV are the corrections:
2855       // dX/dT0 = dX/dxloc * dxloc/dT0 = dX/dxloc * V0
2856       // dX/dV  = dX/dxloc * dxloc/dV =  dX/dxloc * (T-T0)
2857       // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
2858       //
2859       if (!dfDone[AliITSAlignMille2Module::kDOFT0] ||  !dfDone[sddLR]) {
2860         //
2861         double dXdxlocsens=0., dZdxlocsens=0.;
2862         //
2863         // if the current module is the sensor itself and we work with local params, then 
2864         // we can directly take dX/dxloc_sens dZ/dxloc_sens
2865         if (!fUseGlobalDelta && fCurrentModule->GetVolumeID()==fCluster.GetVolumeID()) {
2866           if (!dfDone[AliITSAlignMille2Module::kDOFTX]) {
2867             CalcDerivatives(AliITSAlignMille2Module::kDOFTX,kFALSE); 
2868             dfDone[AliITSAlignMille2Module::kDOFTX] = kTRUE;
2869           }
2870           dXdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][0];
2871           dZdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][2];
2872         }
2873         //
2874         else { // need to perform some transformations
2875           // fetch the jacobian of the transformation from the sensors local frame to the frame
2876           // where the parameters are defined:
2877           // Global: dX/dxloc_sens = dX/dxgl*dxgl/dxloc_sens + ...dX/dphigl*dphigl/dxloc_sens
2878           if (fUseGlobalDelta) fCurrentModule->CalcDerivGloLoc(fCluster.GetVolumeID(),
2879                                                                AliITSAlignMille2Module::kDOFTX, jacobian);
2880           // Local:  dX/dxloc_sens = dX/dxcurr*dxcurr/dxloc_sens +..+dX/dphicurr * dphicurr/dxloc_sens 
2881           else                 fCurrentModule->CalcDerivCurLoc(fCluster.GetVolumeID(),
2882                                                                AliITSAlignMille2Module::kDOFTX, jacobian);
2883           //
2884           for (int j=0;j<kNParChGeom;j++) {
2885             // need global derivative even if the j-th param is locked
2886             if (!dfDone[j]) {CalcDerivatives(j,kFALSE); dfDone[j] = kTRUE;}
2887             dXdxlocsens += fDerivativeGlo[j][0] * jacobian[j];
2888             dZdxlocsens += fDerivativeGlo[j][2] * jacobian[j];
2889           }
2890         }
2891         //
2892         if (zeroX) zeroX = IsZero(dXdxlocsens);
2893         if (zeroZ) zeroZ = IsZero(dZdxlocsens);
2894         //
2895         double vdrift = GetVDriftSDD();
2896         double tdrift = GetTDriftSDD();
2897         //
2898         fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0] = dXdxlocsens*vdrift;
2899         fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2] = dZdxlocsens*vdrift;
2900         dfDone[AliITSAlignMille2Module::kDOFT0] = kTRUE;
2901         //
2902         double mltCorr = fIsSDDVDriftMult ? TMath::Abs(vdrift) : 1;
2903         fDerivativeGlo[sddLR][0] = -dXdxlocsens*mltCorr*TMath::Sign(tdrift,vdrift);
2904         fDerivativeGlo[sddLR][2] = -dZdxlocsens*mltCorr*TMath::Sign(tdrift,vdrift);
2905         dfDone[sddLR] = kTRUE;
2906         //
2907       }
2908       //
2909       if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
2910         m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0];
2911         m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2];
2912         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);      
2913       }
2914       //
2915       if (fCurrentModule->GetParOffset(sddLR)>=0) {
2916         m.fDerGlo[ifill][kX] = fDerivativeGlo[sddLR][0];
2917         m.fDerGlo[ifill][kZ] = fDerivativeGlo[sddLR][2];
2918         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(sddLR);      
2919       }
2920     }
2921     //
2922     m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
2923   } while( (fCurrentModule=fCurrentModule->GetParent()) );
2924   //
2925   if (nmodTested>0 && zeroX) {AliInfo("Skipping: zero global X derivatives!");return -1;}
2926   if (nmodTested>0 && zeroZ) {AliInfo("Skipping: zero global Z derivatives!");return -1;}
2927   //
2928   // ok, can copy to m
2929   AliDebug(2,Form("Adding local equation X with fMeas=%.6f  and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
2930   m.fMeas[kX] = fMeasLoc[0]-fPintLoc0[0];
2931   m.fSigma[kX] = fSigmaLoc[0];
2932   //
2933   AliDebug(2,Form("Adding local equation Z with fMeas=%.6f  and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
2934   m.fMeas[kZ] = fMeasLoc[2]-fPintLoc0[2];
2935   m.fSigma[kZ] = fSigmaLoc[2];
2936   //
2937   m.fNGlobFilled = ifill;
2938   fCurrentModule = endModule;
2939   //
2940   status += Int_t(!zeroX && !zeroZ); // 0 - only locals, 1 locals + globals
2941   return status;
2942 }
2943
2944 //________________________________________________________________________________________________________
2945 Int_t AliITSAlignMille2::AddLocalEquationTPA(Mille2Data &m) 
2946 {
2947   /// Define local equation for current cluster in X Y and Z coor.
2948   /// and store them to memory
2949   /// return -1 in case of failure to build some equation
2950   ///         0 if no free global parameters were found but local eq is built
2951   ///         1 if both local and global eqs are built
2952   //
2953   static int cnt = 0;
2954   Bool_t dbg = kFALSE;//kTRUE;
2955   if (++cnt>100000) dbg = kFALSE;
2956
2957   int curpoint = fCluster.GetUniqueID()-1;
2958   TGeoHMatrix *tempHMat = GetSensorCurrMatrixSID(fCurrentSensID);// fCurrentModule->GetSensitiveVolumeMatrix(fCluster.GetVolumeID());
2959   //
2960   fTPAFitter->GetDResDParams(&fDerivativeLoc[0][0], curpoint);    // resid. derivatives over the track parameters 
2961   if (fCurvFitWasConstrained && fFixCurvIfConstraned && !IsZero(fBField)) 
2962     for (int i=3;i--;) fDerivativeLoc[AliITSTPArrayFit::kR0][i] = 0; //Fix curvarute
2963   //
2964   for (Int_t i=fNLocal; i--;) tempHMat->MasterToLocalVect(fDerivativeLoc[i],m.fDerLoc[i]); 
2965   //
2966   int status = 0;
2967   // derivatives over the global parameters ---------------------------------------->>>
2968   Double_t dGL[3];     // derivative of global position vs local X (for SDD)
2969   Double_t dRdP[3][3]; // derivative of local residuals vs local position
2970   Double_t dPdG[AliITSAlignMille2Module::kMaxParGeom][3]; // derivatives of local position vs global params
2971   fTPAFitter->GetDResDPos(&fDerivativeGlo[0][0], curpoint);
2972   if (fCurrentSensID!=kVtxSensID) for (int i=3;i--;) tempHMat->MasterToLocalVect(fDerivativeGlo[i],dRdP[i]);
2973   else for (int i=3;i--;) for (int j=3;j--;) dRdP[i][j] = fDerivativeGlo[i][j];  // vertex constraint is in lab
2974   //
2975   if (dbg) {
2976     printf("\nCurrentMod: %s Sens:%d\n",fCurrentModule->GetName(),fCurrentSensID); //RRR
2977     printf("Module Matrix: ");
2978     fCurrentModule->GetMatrix()->Print(); //RRR
2979     for (int i=0;i<3;i++) {
2980       printf("dRdP[M%d][resI] ",i); for (int j=0;j<3;j++) printf(":[%d] %+.3e ",j,dRdP[i][j]); printf("\n");
2981     }//RRR
2982     printf("Sensor Matrix: "); tempHMat->Print();
2983   }
2984   UInt_t ifill=0, dfDone = 0;
2985   m.fNModFilled = 0;
2986   // 
2987   AliITSAlignMille2Module* endModule = fCurrentModule;
2988   //
2989   m.fModuleID[0] = fCurrentModule->GetUniqueID(); // always register id of the base module, even if it has no DOF
2990   //
2991   do {
2992     if (fCurrentModule->GetNParFree()==0) continue;
2993     status = 1;
2994     if (!fUseGlobalDelta) dfDone = 0; // for local deltas the derivatives at diff. levels are different
2995     Bool_t jacobOK = kFALSE;
2996     //
2997     for (Int_t i=0; i<kNParChGeom; i++) {              // common for all sensors: derivatives over geom params
2998       if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2999       //
3000       if (!TestWordBit(dfDone,i)) {                    // need to calculate new derivative
3001         if (!jacobOK) {
3002           if (fCurrentSensID!=kVtxSensID) {
3003             fCurrentModule->CalcDerivDPosDPar(fCluster.GetVolumeID(),fMeasLoc,&dPdG[0][0]); 
3004             if (dbg) {
3005               for (int i1=0;i1<3;i1++) {
3006                 printf("Jacob:dPdG[gpar%d][Mj]",i1); for (int j1=0;j1<3;j1++) printf(":[%d] %+.3e ",j1,dPdG[i1][j1]); printf("\n");//RRR
3007               }
3008             }
3009           }
3010           else {
3011             // this is a vertex constraint: only lateral shifts are allowed, no rotations
3012             for (int ip=AliITSAlignMille2Module::kMaxParGeom;ip--;) for (int jp=3;jp--;) dPdG[ip][jp] = (ip==jp) ? 1:0;   
3013           }
3014           jacobOK = kTRUE;
3015         }       
3016         // dRes_j/dGlo_i = \sum_{k=1:3}  dRes_j/dPos_k * dPos_k/dGlo_i
3017         fDerivativeGlo[i][kX] = dRdP[kX][kX]*dPdG[i][kX] + dRdP[kY][kX]*dPdG[i][kY] + dRdP[kZ][kX]*dPdG[i][kZ];
3018         fDerivativeGlo[i][kY] = dRdP[kX][kY]*dPdG[i][kX] + dRdP[kY][kY]*dPdG[i][kY] + dRdP[kZ][kY]*dPdG[i][kZ];
3019         fDerivativeGlo[i][kZ] = dRdP[kX][kZ]*dPdG[i][kX] + dRdP[kY][kZ]*dPdG[i][kY] + dRdP[kZ][kZ]*dPdG[i][kZ];
3020         SetWordBit(dfDone,i);
3021       }
3022       //
3023       if (dbg) {
3024         printf("Level %s DGlob[par%d][resJ] ",fCurrentModule->GetName(),i); //RRR
3025         for (int k=0;k<3;k++) printf(":[%d] %+.3e ",k, fDerivativeGlo[i][k]); printf("\n");//RRR
3026       }
3027       m.fDerGlo[ifill][kX] = fDerivativeGlo[i][kX];
3028       m.fDerGlo[ifill][kY] = fDerivativeGlo[i][kY];
3029       m.fDerGlo[ifill][kZ] = fDerivativeGlo[i][kZ];
3030       m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
3031       //
3032     }
3033     //
3034     if ( fCurrentModule->IsSDD() ) {     // specific for SDD
3035       //
3036       // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
3037       // where V0 and T are the nominal drift velocity, time and time0
3038       // and the dT0 and dV are the corrections:
3039       // drloc_i/dT0 = sum_j drloc_i/dMeasGlo_j * dMeasGlo_j/dT0 = 
3040       //             = sum_j drloc_i/dMeasGlo_j sum_k dMeasGlo_j/dMeasLoc_k * dMeasLoc_k/dT0
3041       //             = sum_j drloc_i/dMeasGlo_j dMeasGlo_j/dMeasLoc_X * V0
3042       //
3043       // drloc_i/dV0 = sum_j drloc_i/dMeasGlo_j * dMeasGlo_j/dV0 = 
3044       //             = sum_j drloc_i/dMeasGlo_j sum_k dMeasGlo_j/dMeasLoc_k * dMeasLoc_k/dV0
3045       //             = sum_j drloc_i/dMeasGlo_j dMeasGlo_j/dMeasLoc_X * T0
3046
3047       // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
3048       //
3049       Bool_t jacOK = kFALSE;
3050       //Int_t sddLR = fMeasLoc[kX]>0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR;
3051       Int_t sddLR = GetVDriftSDD()>0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR;
3052       if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
3053         if (!TestWordBit(dfDone, AliITSAlignMille2Module::kDOFT0)) {
3054           double vdrift = GetVDriftSDD();
3055           JacobianPosGloLoc(kX,dGL);
3056           jacOK = kTRUE;
3057           fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kX] = 
3058             vdrift*(dRdP[kX][kX]*dGL[kX] + dRdP[kY][kX]*dGL[kY] + dRdP[kZ][kX]*dGL[kZ]);
3059           fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kY] = 
3060             vdrift*(dRdP[kX][kY]*dGL[kX] + dRdP[kY][kY]*dGL[kY] + dRdP[kZ][kY]*dGL[kZ]);
3061           fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kZ] = 
3062             vdrift*(dRdP[kX][kZ]*dGL[kX] + dRdP[kY][kZ]*dGL[kY] + dRdP[kZ][kZ]*dGL[kZ]);
3063           //
3064           SetWordBit(dfDone, AliITSAlignMille2Module::kDOFT0);
3065         }
3066         m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kX];
3067         m.fDerGlo[ifill][kY] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kY];
3068         m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kZ];
3069         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);      
3070       }
3071       //
3072       if (fCurrentModule->GetParOffset(sddLR)>=0) {
3073         if (!TestWordBit(dfDone, sddLR)) {
3074           double tdrift = TMath::Sign(GetTDriftSDD(), GetVDriftSDD());
3075           double vdrift = fIsSDDVDriftMult ? TMath::Abs(GetVDriftSDD()) : 1;
3076           if (!jacOK) JacobianPosGloLoc(kX,dGL);
3077           fDerivativeGlo[sddLR][kX] = 
3078             -tdrift*vdrift*(dRdP[kX][kX]*dGL[kX] + dRdP[kY][kX]*dGL[kY] + dRdP[kZ][kX]*dGL[kZ]);
3079           fDerivativeGlo[sddLR][kY] = 
3080             -tdrift*vdrift*(dRdP[kX][kY]*dGL[kX] + dRdP[kY][kY]*dGL[kY] + dRdP[kZ][kY]*dGL[kZ]);
3081           fDerivativeGlo[sddLR][kZ] = 
3082             -tdrift*vdrift*(dRdP[kX][kZ]*dGL[kX] + dRdP[kY][kZ]*dGL[kY] + dRdP[kZ][kZ]*dGL[kZ]);
3083           SetWordBit(dfDone, sddLR);
3084         }
3085         m.fDerGlo[ifill][kX] = fDerivativeGlo[sddLR][kX];
3086         m.fDerGlo[ifill][kY] = fDerivativeGlo[sddLR][kY];
3087         m.fDerGlo[ifill][kZ] = fDerivativeGlo[sddLR][kZ];
3088         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(sddLR);      
3089       }
3090     }
3091     //
3092     m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
3093   } while( (fCurrentModule=fCurrentModule->GetParent()) );
3094   //
3095   // store first local residuals
3096   fTPAFitter->GetResiduals(fPintLoc , curpoint);       // lab residuals
3097   for (int i=3;i--;) fPintLoc[i] = -fPintLoc[i];
3098   if   (fCurrentSensID!=kVtxSensID) tempHMat->MasterToLocalVect(fPintLoc,m.fMeas);       // local residuals 
3099   else for (int i=3;i--;) m.fMeas[i] = fPintLoc[i];
3100   if (dbg) {
3101     printf("res(meas-loc) "); for (int k=0;k<3;k++) printf(":[%d] %+.3e ",k,m.fMeas[k]); printf("\n");
3102     printf("Fin:%s %+e %+e\n",endModule->GetName(), fDerivativeGlo[kZ][kZ], fPintLoc[kZ]);
3103   }//RRR
3104   m.fSigma[kX] = fSigmaLoc[kX];
3105   m.fSigma[kY] = fSigmaLoc[kY];
3106   m.fSigma[kZ] = fSigmaLoc[kZ];
3107   //
3108   m.fNGlobFilled = ifill;
3109   fCurrentModule = endModule;
3110   //
3111   return status;
3112 }
3113
3114 //________________________________________________________________________________________________________
3115 void AliITSAlignMille2::SetLocalEquations(const Mille2Data *marr, Int_t neq) 
3116 {
3117   /// Set local equations with data stored in m
3118   /// return 0 if success
3119   //
3120   Bool_t locPatt[kNLocal] = {0}; // pattern of lacal eq's to account
3121   for (int i=fNLocal; i--;) {
3122     if (locPatt[i]) continue; // already set
3123     for (Int_t j=0; j<neq; j++) for (int ic=3;ic--;) if (!IsZero(marr[j].fDerLoc[i][ic])) locPatt[i]=kTRUE;
3124   }
3125   //
3126   for (Int_t j=0; j<neq; j++) {
3127     //
3128     const Mille2Data &m = marr[j];
3129     //
3130     Bool_t filled = kFALSE;
3131     for (int ic=3;ic--;) {
3132       // for the diamond point (if any) the Y residual is accounted
3133       if (ic==kY && !fUseLocalYErr && !(m.fModuleID[0]==fDiamondModID)) continue;
3134       AliDebug(2,Form("setting local equation %c with fMeas=%.6f  and fSigma=%.6f",fgkXYZ[ic],m.fMeas[ic], m.fSigma[ic]));      
3135       Int_t nzero = 0, naddl = 0;
3136       for (int i=0;i<=fNLocal;i++) if (locPatt[i]) nzero += SetLocalDerivative(naddl++,m.fDerLoc[i][ic] );
3137       if (nzero==fNLocal) { 
3138         AliInfo(Form("Skipping %c residual due to the zero derivatives!",fgkXYZ[ic])); 
3139         continue; 
3140       }
3141       for (int i=m.fNGlobFilled;i--;) SetGlobalDerivative( m.fParMilleID[i] , m.fDerGlo[i][ic] );
3142       fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.fMeas[ic], m.fSigma[ic]);  
3143       filled = kTRUE;
3144       //
3145     }
3146     //
3147     if (filled) for (int i=m.fNModFilled;i--;) GetMilleModule(m.fModuleID[i])->IncNProcessedPoints();
3148   }
3149   //
3150   double wgh = 1;
3151   if (GetWeightPt() && fTPAFitter) {
3152     wgh = fTPAFitter->GetPt();
3153     if (wgh>10) wgh = 10.;
3154     if (wgh<0) wgh = fTPAFitter->IsTypeCosmics() ? 7 : 0.5;
3155     if (GetWeightPt()>0) wgh = TMath::Power(wgh,GetWeightPt());
3156   }
3157   fMillepede->SetRecordWeight(wgh*fTrackWeight);
3158   fMillepede->SetRecordRun(fRunID);
3159   //
3160 }
3161
3162 //________________________________________________________________________________________________________
3163 Int_t AliITSAlignMille2::GlobalFit()
3164 {
3165   /// Call global fit; Global parameters are stored in parameters
3166   if (!fIsMilleInit) Init();
3167   //
3168   ApplyPreConstraints();
3169   int res = fMillepede->GlobalFit();
3170   AliInfo(Form("%s fitting global parameters!",res ? "Done":"Failed"));
3171   if (res) {
3172     // fetch the parameters
3173     for (int imd=fNModules;imd--;) {
3174       AliITSAlignMille2Module* mod = GetMilleModule(imd);
3175       int nprocp = 0;
3176       for (int ip=mod->GetNParTot();ip--;) {
3177         int idp = mod->GetParOffset(ip);
3178         if (idp<0) continue;    // was not in the explicit fit
3179         mod->SetParVal(ip,fMillepede->GetFinalParam(idp));
3180         mod->SetParErr(ip,fMillepede->GetFinalError(idp));
3181         int np = fMillepede->GetProcessedPoints(idp);
3182         if (TMath::Abs(np)>TMath::Abs(nprocp)) nprocp = np;
3183       }
3184       if (!mod->GetNProcessedPoints()) mod->SetNProcessedPoints(nprocp);
3185     }
3186
3187   }
3188   ApplyPostConstraints();
3189   return res;
3190 }
3191
3192 //________________________________________________________________________________________________________
3193 void AliITSAlignMille2::PrintGlobalParameters() 
3194 {
3195   /// Print global parameters
3196   if (!fIsMilleInit) {
3197     AliInfo("Millepede has not been initialized!");
3198     return;
3199   }
3200   fMillepede->PrintGlobalParameters();
3201 }
3202
3203 //________________________________________________________________________________________________________
3204 Int_t AliITSAlignMille2::LoadSuperModuleFile(const Char_t *sfile)
3205
3206   // load definitions of supermodules from a root file
3207   // return 0 if success
3208   AliInfo(Form("Loading SuperModule definitions from %s",sfile));
3209   TFile *smf=TFile::Open(sfile);
3210   if (!smf->IsOpen()) {
3211     AliInfo(Form("Cannot open supermodule file %s",sfile));
3212     return -1;
3213   }
3214
3215   TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules");
3216   if (!sma) {
3217     AliInfo(Form("Cannot find ITSMilleSuperModules array in file"));
3218     return -2;  
3219   }  
3220   Int_t nsma=sma->GetEntriesFast();
3221   AliInfo(Form("Array of SuperModules with %d entries\n",nsma));
3222   //
3223   // pepo200709
3224   Char_t st[2048];
3225   char symname[250];
3226   // end pepo200709
3227
3228   UShort_t volid;
3229   TGeoHMatrix m;
3230   //
3231   for (Int_t i=0; i<nsma; i++) {
3232     AliAlignObjParams *a = (AliAlignObjParams*)sma->UncheckedAt(i);
3233     volid=a->GetVolUID();
3234     strncpy(st,a->GetSymName(),TMath::Min(sizeof(st),strlen(a->GetSymName())+1));
3235     a->GetMatrix(m);
3236     //
3237     symname[0] = '\0';
3238     sscanf(st,"%249s",symname);
3239     //
3240     // decode module list
3241     char *stp=strstr(st,"ModuleList:");
3242     if (!stp) return -3;
3243     stp += 11;
3244     int idx[2200];
3245     char spp[200]; int jp=0;
3246     char cl[20];
3247     strncpy(st,stp,TMath::Min(sizeof(st),strlen(stp)+1));
3248     int l=strlen(st);
3249     int j=0;
3250     int n=0;
3251     //
3252     while (j<=l) {
3253       if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) {
3254         spp[jp]=0;
3255         jp=0;
3256         if (strlen(spp)) {
3257           int k=strcspn(spp,"-");
3258           if (k<int(strlen(spp))) { // c'e' il -
3259             strncpy(cl,&(spp[k+1]), TMath::Min(sizeof(cl),strlen(&spp[k+1])+1));
3260             spp[k]=0;
3261             int ifrom=atoi(spp); int ito=atoi(cl);
3262             for (int b=ifrom; b<=ito; b++) {
3263               idx[n]=b;
3264               n++;
3265             }
3266           }
3267           else { // numerillo singolo
3268             idx[n]=atoi(spp);
3269             n++;
3270           }
3271         }
3272       }
3273       else {
3274         spp[jp]=st[j];
3275         jp++;
3276       }
3277       j++;
3278     }
3279     UShort_t volidsv[2198];
3280     for (j=0;j<n;j++) {
3281       volidsv[j]=AliITSAlignMille2Module::GetVolumeIDFromIndex(idx[j]);
3282       if (!volidsv[j]) {
3283         AliInfo(Form("Index %d not valid (range 0->%d)",idx[j],kMaxITSSensID));
3284         return -5;
3285       }
3286     }
3287     Int_t smindex=int(2198+volid-14336); // virtual index
3288     //
3289     fSuperModule.AddAtAndExpand(new AliITSAlignMille2Module(smindex,volid,symname,&m,n,volidsv),fNSuperModules);
3290     //
3291     fNSuperModules++;
3292   }
3293   //
3294   smf->Close();
3295   //
3296   return 0;
3297 }
3298
3299 //________________________________________________________________________________________________________
3300 void AliITSAlignMille2::ConstrainModuleSubUnitsMean(Int_t idm, Double_t val, UInt_t pattern)
3301 {
3302   // require that sum of modifications for the childs of this module is = val, i.e.
3303   // the internal corrections moves the module as a whole by fixed value (0 by default).
3304   // pattern is the bit pattern for the parameters to constrain
3305   //
3306   if (fIsMilleInit) {
3307     AliInfo("Millepede has been already initialized: no constrain may be added!");
3308     return;
3309   }
3310   if (!GetMilleModule(idm)->GetNChildren()) return;
3311   TString nm = "cstrSUMean";
3312   nm += GetNConstraints();
3313   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
3314                                                                       idm,val,pattern);
3315   cstr->SetConstraintID(GetNConstraints());
3316   fConstraints.Add(cstr);
3317 }
3318
3319 //________________________________________________________________________________________________________
3320 void AliITSAlignMille2::ConstrainModuleSubUnitsMedian(Int_t idm, Double_t val, UInt_t pattern)
3321 {
3322   // require that median of the modifications for the childs of this module is = val, i.e.
3323   // the internal corrections moves the module as a whole by fixed value (0 by default) 
3324   // module the outliers.
3325   // pattern is the bit pattern for the parameters to constrain
3326   // The difference between the mean and the median will be transfered to the parent
3327   if (fIsMilleInit) {
3328     AliInfo("Millepede has been already initialized: no constrain may be added!");
3329     return;
3330   }
3331   if (!GetMilleModule(idm)->GetNChildren()) return;
3332   TString nm = "cstrSUMed";
3333   nm += GetNConstraints();
3334   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
3335                                                                       idm,val,pattern);
3336   cstr->SetConstraintID(GetNConstraints());
3337   fConstraints.Add(cstr);
3338 }
3339
3340 //________________________________________________________________________________________________________
3341 void AliITSAlignMille2::ConstrainOrphansMean(Double_t val, UInt_t pattern)
3342 {
3343   // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3344   // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3345   // pattern is the bit pattern for the parameters to constrain
3346   //
3347   if (fIsMilleInit) {
3348     AliInfo("Millepede has been already initialized: no constrain may be added!");
3349     return;
3350   }
3351   TString nm = "cstrOMean";
3352   nm += GetNConstraints();
3353   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
3354                                                                       -1,val,pattern);
3355   cstr->SetConstraintID(GetNConstraints());
3356   fConstraints.Add(cstr);
3357 }
3358
3359 //________________________________________________________________________________________________________
3360 void AliITSAlignMille2::ConstrainOrphansMedian(Double_t val, UInt_t pattern)
3361 {
3362   // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3363   // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3364   // pattern is the bit pattern for the parameters to constrain
3365   //
3366   if (fIsMilleInit) {
3367     AliInfo("Millepede has been already initialized: no constrain may be added!");
3368     return;
3369   }
3370   TString nm = "cstrOMed";
3371   nm += GetNConstraints();
3372   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
3373                                                                       -1,val,pattern);
3374   cstr->SetConstraintID(GetNConstraints());
3375   fConstraints.Add(cstr);
3376 }
3377
3378 //________________________________________________________________________________________________________
3379 void AliITSAlignMille2::ConstrainLocal(const Char_t* name,Double_t *parcf,Int_t npar,Double_t val,Double_t err)
3380 {
3381   // apply constraint on parameters in the local frame
3382   if (fIsMilleInit) {
3383     AliInfo("Millepede has been already initialized: no constrain may be added!");
3384     return;
3385   }
3386   AliITSAlignMille2ConstrArray *cstr = new AliITSAlignMille2ConstrArray(name,parcf,npar,val,err);
3387   cstr->SetConstraintID(GetNConstraints());
3388   fConstraints.Add(cstr);
3389 }
3390
3391 //________________________________________________________________________________________________________
3392 void AliITSAlignMille2::ApplyGaussianConstraint(const AliITSAlignMille2ConstrArray* cstr)
3393 {
3394   // apply the constraint on the local corrections of a list of modules
3395   int nmod = cstr->GetNModules();
3396   double jacobian[AliITSAlignMille2Module::kMaxParGeom][AliITSAlignMille2Module::kMaxParGeom];
3397   //
3398   // check if this not special SDDT0 constraint
3399   if (cstr->GetPattern()==BIT(AliITSAlignMille2Module::kDOFT0)) {
3400     for (int i=0;i<cstr->GetNModules()-1;i++) {
3401       AliITSAlignMille2Module *mdI = GetMilleModule(cstr->GetModuleID(i));
3402       if (!mdI->IsFreeDOF(AliITSAlignMille2Module::kDOFT0)) continue;
3403       for (int j=i+1;j<cstr->GetNModules();j++) {
3404         AliITSAlignMille2Module *mdJ = GetMilleModule(cstr->GetModuleID(j));
3405         if (!mdJ->IsFreeDOF(AliITSAlignMille2Module::kDOFT0)) continue;
3406         //
3407         ResetLocalEquation();
3408         fGlobalDerivatives[mdI->GetParOffset(AliITSAlignMille2Module::kDOFT0)] = 1;
3409         fGlobalDerivatives[mdJ->GetParOffset(AliITSAlignMille2Module::kDOFT0)] =-1;
3410         AddConstraint(fGlobalDerivatives, 0, 1.E-6);
3411       }
3412     }
3413     return;
3414   }
3415
3416   for (int imd=nmod;imd--;) {
3417     int modID = cstr->GetModuleID(imd);
3418     AliITSAlignMille2Module* mod = GetMilleModule(modID);
3419     ResetLocalEquation();
3420     int nadded = 0;
3421     double value = cstr->GetValue();
3422     double sigma = cstr->GetError();
3423     //
3424     // in case the reference (survey) deltas were imposed for Gaussian constraints