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