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