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