ITS/AliITSAlignMille2.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-----------------------------------------------------------------------------
  19 //
  20 //  Interface to AliMillePede2 alignment class for the ALICE ITS detector
  21 //
  22 //  ITS specific alignment class which interface to AliMillepede.
  23 //  For each track ProcessTrack calculates the local and global derivatives
  24 //  at each hit and fill the corresponding local equations. Provide methods for
  25 //  fixing or constraning detection elements for best results.
  26 //
  27 //  author M. Lunardon (thanks to J. Castillo), ruben.shahoyan@cern.ch
  28 //-----------------------------------------------------------------------------
  29
  30 #include <TFile.h>
  31 #include <TGrid.h>
  32 #include <TClonesArray.h>
  33 #include <TMath.h>
  34 #include <TVirtualFitter.h>
  35 #include <TGeoManager.h>
  36 #include <TSystem.h>
  37 #include <TRandom.h>
  38 #include <TCollection.h>
  39 #include <TGeoPhysicalNode.h>
  40 #include <TMap.h>
  41 #include <TObjString.h>
  42 #include <TString.h>
  43 #include "AliITSAlignMille2.h"
  44 #include "AliITSgeomTGeo.h"
  45 #include "AliGeomManager.h"
  46 #include "AliMillePede2.h"
  47 #include "AliTrackPointArray.h"
  48 #include "AliAlignObjParams.h"
  49 #include "AliLog.h"
  50 #include "AliTrackFitterRieman.h"
  51 #include "AliITSAlignMille2Constraint.h"
  52 #include "AliITSAlignMille2ConstrArray.h"
  53 #include "AliITSresponseSDD.h"
  54 #include "AliITSTPArrayFit.h"
  55 #include "AliCDBManager.h"
  56 #include "AliCDBStorage.h"
  57 #include "AliCDBEntry.h"
  58 #include "AliITSsegmentationSDD.h"
  59 #include "AliITSDriftSpeedArraySDD.h"
  60 #include "AliITSCorrectSDDPoints.h"
  61 #include "AliESDVertex.h"
  62
  63 ClassImp(AliITSAlignMille2)
  64
  65 const Char_t* AliITSAlignMille2::fgkRecKeys[] = {
  66   "OCDB_PATH",
  67   "OCDB_SPECIFIC",
  68   "GEOMETRY_FILE",
  69   "SUPERMODULE_FILE",
  70   "CONSTRAINTS_REFERENCE_FILE",
  71   "PREALIGNMENT_FILE",
  72   "PRECALIBSDD_FILE",
  73   "PREVDRIFTSDD_FILE",
  74   "PRECORRMAPSDD_FILE",
  75   "INITCORRMAPSDD_FILE",
  76   "INITCALBSDD_FILE",
  77   "INITVDRIFTSDD_FILE",
  78   "INITDELTA_FILE",
  79   "INITGEOM_FILE",
  80   "SET_GLOBAL_DELTAS",
  81   "CONSTRAINT_LOCAL",
  82   "MODULE_VOLUID",
  83   "MODULE_INDEX",
  84   "SET_PSEUDO_PARENTS",
  85   "SET_TRACK_FIT_METHOD",
  86   "SET_MINPNT_TRA",
  87   "SET_NSTDDEV",
  88   "SET_RESCUT_INIT",
  89   "SET_RESCUT_OTHER",
  90   "SET_LOCALSIGMAFACTOR",
  91   "SET_STARTFAC",
  92   "SET_FINALFAC",
  93   "SET_B_FIELD",
  94   "SET_SPARSE_MATRIX",
  95   "REQUIRE_POINT",
  96   "CONSTRAINT_ORPHANS",
  97   "CONSTRAINT_SUBUNITS",
  98   "APPLY_CONSTRAINT",
  99   "SET_EXTRA_CLUSTERS_MODE",
 100   "SET_USE_TPAFITTER",
 101   "SET_USE_LOCAL_YERROR",
 102   "SET_MIN_POINTS_PER_MODULE",
 103   "SET_USE_SDDVDCORRMULT",
 104   "SET_WEIGHT_PT",
 105   "SET_USE_DIAMOND",
 106   "CORRECT_DIAMOND",
 107   "SET_USE_VERTEX",
 108   "SET_SAME_SDDT0"
 109 };
 110
 111 const Char_t AliITSAlignMille2::fgkXYZ[] = "XYZ";
 112
 113 //========================================================================================================
 114
 115 AliITSAlignMille2* AliITSAlignMille2::fgInstance = 0;
 116 Int_t              AliITSAlignMille2::fgInstanceID = 0;
 117
 118 //________________________________________________________________________________________________________
 119 AliITSAlignMille2::AliITSAlignMille2(const Char_t *configFilename,TList *userInfo )
 120 : TObject(),
 121   fMillepede(0),
 122   fStartFac(16.),
 123   fFinalFac(1.),
 124   fResCutInitial(100.),
 125   fResCut(100.),
 126   fNGlobal(0),
 127   fNLocal(4),
 128   fNStdDev(3),
 129   fIsMilleInit(kFALSE),
 130   fAllowPseudoParents(kFALSE),
 131   //
 132   fTPAFitter(0),
 133   fCurrentModule(0),
 134   fTrack(0),
 135   fTrackBuff(0),
 136   fCluster(),
 137   fCurrentSensID(-1),
 138   fClusLoc(12*3),
 139   fClusGlo(12*3),
 140   fClusSigLoc(12*3),
 141   fGlobalDerivatives(0),
 142   fMeasLoc(0),
 143   fMeasGlo(0),
 144   fSigmaLoc(0),
 145   fConstrPT(-1),
 146   fConstrPTErr(-1),
 147   fConstrCharge(0),
 148   fRunID(0),
 149   //
 150   fMinNPtsPerTrack(3),
 151   fIniTrackParamsMeth(1),
 152   fTotBadLocEqPoints(0),
 153   fRieman(0),
 154   //
 155   fConstraints(0),
 156   fCacheMatrixOrig(kMaxITSSensID+1),
 157   fCacheMatrixCurr(kMaxITSSensID+1),
 158   //
 159   fUseGlobalDelta(kFALSE),
 160   fTempExcludedModule(-1),
 161   fUserProvided(0),
 162   //
 163   fIniUserInfo(userInfo),
 164   fIniGeomPath(""),
 165   fIniDeltaPath(""),
 166   fIniSDDRespPath(""),
 167   fPreCalSDDRespPath(""),
 168   fIniSDDVDriftPath(""),
 169   fPreSDDVDriftPath(""),
 170   fIniSDDCorrMapPath(""),
 171   fPreSDDCorrMapPath(""),
 172   fConvertPreDeltas(kFALSE),
 173   fGeometryPath(""),
 174   fPreDeltaPath(""),
 175   fConstrRefPath(""),
 176   fDiamondPath(""),
 177   fGeoManager(0),
 178   fIsConfigured(kFALSE),
 179   fPreAlignQF(0),
 180 //
 181   fIniRespSDD(0),
 182   fPreRespSDD(0),
 183   fIniVDriftSDD(0),
 184   fPreVDriftSDD(0),
 185   fIniCorrMapSDD(0),
 186   fPreCorrMapSDD(0),
 187   fSegmentationSDD(0),
 188   fPrealignment(0),
 189   fConstrRef(0),
 190   fMilleModule(2),
 191   fSuperModule(2),
 192   fNModules(0),
 193   fNSuperModules(0),
 194   fUsePreAlignment(kFALSE),
 195   fUseLocalYErr(kFALSE),
 196   fBOn(kFALSE),
 197   fBField(0.0),
 198   fDataType(kCosmics),
 199   fMinPntPerSens(0),
 200   fBug(0),
 201   fMilleVersion(2),
 202   fExtraClustersMode(0),
 203   fTrackWeight(1),
 204   fWeightPt(0.),
 205   fIsSDDVDriftMult(kFALSE),
 206   fDiamond(),
 207   fDiamondI(),
 208   fUseDiamond(kFALSE),
 209   fUseVertex(kFALSE),
 210   fVertexSet(kFALSE),
 211   fDiamondPointID(-1),
 212   fDiamondModID(-1),
 213   fCheckDiamondPoint(kDiamondCheckIfPrompt),
 214   fConvAlgMatOld(100)
 215 {
 216   /// main constructor that takes input from configuration file
 217   for (int i=3;i--;) fSigmaFactor[i] = 1.0;
 218   //
 219   // new RS
 220   for (int i=0;i<3;i++) {
 221     fCorrDiamond[i] = 0;
 222   }
 223   for (int itp=0;itp<kNDataType;itp++) {
 224     fRequirePoints[itp] = kFALSE;
 225     for (Int_t i=0; i<6; i++) {
 226       fNReqLayUp[itp][i]=0;
 227       fNReqLayDown[itp][i]=0;
 228       fNReqLay[itp][i]=0;
 229     }
 230     for (Int_t i=0; i<3; i++) {
 231       fNReqDetUp[itp][i]=0;
 232       fNReqDetDown[itp][i]=0;
 233       fNReqDet[itp][i]=0;
 234     }
 235   }
 236   //
 237   //  if (ProcessUserInfo(userInfo)) exit(1);
 238   //
 239   fDiamond.SetVolumeID(kVtxSensVID);
 240   fDiamondI.SetVolumeID(kVtxSensVID);
 241   float xyzd[3] = {0,0,0};
 242   float covd[6] = {1,0,0,1,0,1e4};
 243   fDiamond.SetXYZ(xyzd,covd); // dummy diamond
 244   covd[5] = 1e-4;
 245   fDiamondI.SetXYZ(xyzd,covd);
 246   //
 247   Int_t lc=LoadConfig(configFilename);
 248   if (lc) {
 249     AliError(Form("Error %d loading configuration from %s",lc,configFilename));
 250     exit(1);
 251   }
 252   //
 253   fMillepede = new AliMillePede2();
 254   fgInstance = this;
 255   fgInstanceID++;
 256   ResetCovIScale();
 257   //
 258 }
 259
 260 //________________________________________________________________________________________________________
 261 AliITSAlignMille2::~AliITSAlignMille2()
 262 {
 263   /// Destructor
 264   delete fMillepede;
 265   delete[] fGlobalDerivatives;
 266   delete fRieman;
 267   delete fPrealignment;
 268   delete fConstrRef;
 269   delete fPreRespSDD;
 270   delete fIniRespSDD;
 271   delete fSegmentationSDD;
 272   delete fIniVDriftSDD;
 273   delete fPreVDriftSDD;
 274   delete fIniCorrMapSDD;
 275   delete fPreCorrMapSDD;
 276   delete fTPAFitter;
 277   fCacheMatrixOrig.Delete();
 278   fCacheMatrixCurr.Delete();
 279   fTrackBuff.Delete();
 280   fConstraints.Delete();
 281   fMilleModule.Delete();
 282   fSuperModule.Delete();
 283   if (--fgInstanceID==0) fgInstance = 0;
 284 }
 285
 286 ///////////////////////////////////////////////////////////////////////
 287
 288 //________________________________________________________________________________________________________
 289 TObjArray* AliITSAlignMille2::GetConfigRecord(FILE* stream, TString& recTitle, TString& recOpt, Bool_t rew)
 290 {
 291   // read new record from config file
 292   TString record;
 293   static TObjArray* recElems = 0;
 294   if (recElems) {delete recElems; recElems = 0;}
 295   recOpt = "";
 296   //
 297   TString keyws = recTitle;
 298   if (!keyws.IsNull()) {
 299     keyws.ToUpper();
 300     //    keyws += " ";
 301   }
 302   while (record.Gets(stream)) {
 303     int cmt=record.Index("#");
 304     if (cmt>=0) record.Remove(cmt);  // skip comment
 305     record.ReplaceAll("\t"," ");
 306     record.ReplaceAll("\r"," ");
 307     record.Remove(TString::kBoth,' ');
 308     if (record.IsNull()) continue;      // nothing to decode
 309     if (!keyws.IsNull() && !record.BeginsWith(keyws.Data())) continue; // specific record was requested
 310     //
 311     recElems = record.Tokenize(" ");
 312     recTitle = recElems->At(0)->GetName();
 313     recTitle.ToUpper();
 314     recOpt = recElems->GetLast()>0 ? recElems->At(1)->GetName() : "";
 315     break;
 316   }
 317   if (rew || !recElems) rewind(stream);
 318   return recElems;
 319 }
 320
 321 //________________________________________________________________________________________________________
 322 Int_t AliITSAlignMille2::CheckConfigRecords(FILE* stream)
 323 {
 324   TString record,recTitle;
 325   int lineCnt = 0;
 326   rewind(stream);
 327   while (record.Gets(stream)) {
 328     int cmt=record.Index("#");
 329     lineCnt++;
 330     if (cmt>=0) record.Remove(cmt);  // skip comment
 331     record.ReplaceAll("\t"," ");
 332     record.ReplaceAll("\r"," ");
 333     record.Remove(TString::kBoth,' ');
 334     if (record.IsNull()) continue;   // nothing to decode
 335     // extract keyword
 336     int spc = record.Index(" ");
 337     if (spc>0) recTitle = record(0,spc);
 338     else     recTitle = record;
 339     recTitle.ToUpper();
 340     Bool_t strOK = kFALSE;
 341     for (int ik=kNKeyWords;ik--;) if (recTitle == fgkRecKeys[ik]) {strOK = kTRUE; break;}
 342     if (strOK) continue;
 343     //
 344     AliError(Form("Unknown keyword %s at line %d",
 345                   recTitle.Data(),lineCnt));
 346     return -1;
 347     //
 348   }
 349   //
 350   rewind(stream);
 351   return 0;
 352 }
 353
 354
 355 //________________________________________________________________________________________________________
 356 Int_t AliITSAlignMille2::LoadConfig(const Char_t *cfile)
 357 {
 358   // return 0 if success
 359   //        1 if error in module index or voluid
 360   //
 361   AliInfo(Form("Loading MillePede2 configuration from %s",cfile));
 362   AliCDBManager::Instance()->SetCacheFlag(kFALSE);
 363   FILE *pfc=fopen(cfile,"r");
 364   if (!pfc) return -1;
 365   //
 366   TString record,recTitle,recOpt,recExt;
 367   Int_t nrecElems,irec;
 368   TObjArray *recArr=0;
 369   //
 370   fNModules = 0;
 371   Bool_t stopped = kFALSE;
 372   //
 373   if (CheckConfigRecords(pfc)<0) return -1;
 374   //
 375   while(1) {
 376     //
 377     // ============= 1: we read some important records in predefined order ================
 378     //
 379     recTitle = fgkRecKeys[kOCDBDefaultPath];
 380     if ( GetConfigRecord(pfc,recTitle,recOpt,1) && !recOpt.IsNull() ) {
 381       AliInfo(Form("Configuration sets OCDB default storage to %s",recOpt.Data()));
 382       AliCDBManager::Instance()->SetDefaultStorage( gSystem->ExpandPathName(recOpt.Data()) );
 383       TObjString* objStr = (TObjString*)AliCDBManager::Instance()->GetStorageMap()->GetValue("default");
 384       if (!objStr) {stopped = kTRUE; break;}
 385       objStr->SetUniqueID(1); // mark as user set
 386     }
 387     //
 388     if (fIniUserInfo && ProcessUserInfo(fIniUserInfo)) { AliError("Failed to process intial User Info"); stopped = kTRUE; break;}
 389     //
 390     recTitle = fgkRecKeys[kGeomFile];
 391     if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) fGeometryPath = gSystem->ExpandPathName(recOpt.Data());
 392     if ( LoadGeometry(fGeometryPath) ) { AliError("Failed to find/load target ideal Geometry"); stopped = kTRUE; break;}
 393     //
 394     // Do we use new TrackPointArray fitter ?
 395     recTitle = fgkRecKeys[kTPAFitter];
 396     if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) fTPAFitter = new AliITSTPArrayFit(kNLocal);
 397     //
 398     recTitle = fgkRecKeys[kSuperModileFile];
 399     if ( !GetConfigRecord(pfc,recTitle,recOpt,1) ||
 400          recOpt.IsNull()                         ||
 401          gSystem->ExpandPathName(recOpt)         ||
 402          gSystem->AccessPathName(recOpt.Data())  ||
 403          LoadSuperModuleFile(recOpt.Data()))
 404       { AliError("Failed to find/load SuperModules"); stopped = kTRUE; break;}
 405     //
 406     recTitle = fgkRecKeys[kConstrRefFile];      // LOCAL_CONSTRAINTS are defined wrt these deltas
 407     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
 408       if (recOpt.IsNull() || recOpt=="IDEAL") SetConstraintWrtRef( "IDEAL" );
 409       else {
 410         for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
 411         if ( SetConstraintWrtRef(recOpt.Data()) )
 412           { AliError("Failed to load reference deltas for local constraints"); stopped = kTRUE; break;}
 413       }
 414     }
 415     //
 416     recTitle = fgkRecKeys[kInitGeomFile];
 417     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1))  && !recOpt.IsNull() ) {
 418       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
 419       fIniGeomPath = recOpt;
 420       gSystem->ExpandPathName(fIniGeomPath);
 421       fUserProvided |= kSameInitGeomBit;
 422       AliInfo(Form("Configuration sets Production Geometry to %s",fIniGeomPath.Data()));
 423     }
 424     if (fIniGeomPath.IsNull()) fIniGeomPath = fGeometryPath;
 425     //
 426     recTitle = fgkRecKeys[kInitDeltaFile];
 427     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1))  && !recOpt.IsNull() ) {
 428       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
 429       fIniDeltaPath = recOpt;
 430       gSystem->ExpandPathName(fIniDeltaPath);
 431       fUserProvided |= kSameInitDeltasBit;
 432       AliInfo(Form("Configuration sets Production Deltas to %s",fIniDeltaPath.Data()));
 433     }
 434     //
 435     recTitle = fgkRecKeys[kPreDeltaFile];
 436     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
 437       if (!recOpt.IsNull()) {
 438         for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
 439         fPreDeltaPath = recOpt;
 440         gSystem->ExpandPathName(fPreDeltaPath);
 441       }
 442       else if (!fIniDeltaPath.IsNull()) {
 443         AliInfo("PreAlignment Deltas keyword is present but empty, will set to Init Deltas of the first tree");
 444         fPreDeltaPath = fIniDeltaPath;
 445         if (fIniGeomPath != fGeometryPath) fConvertPreDeltas = kTRUE; // production and target geometries are different, request conversion
 446       }
 447       AliInfo(Form("Configuration sets PreAlignment Deltas to %s",fPreDeltaPath.Data()));
 448     }
 449     //
 450     // if initial deltas were provided, load them, apply to geometry and store are "original" matrices
 451     if (CacheMatricesOrig()) {stopped = kTRUE; break;}
 452     //
 453     // then load prealignment deltas
 454     if (!fPreDeltaPath.IsNull()) {
 455       if (fConvertPreDeltas) ConvertDeltas();   // Prealignment deltas are the same as production ones, but need conversion to new geometry
 456       else if (LoadDeltas(fPreDeltaPath,fPrealignment)) {stopped = kTRUE; break;} // read deltas from the file
 457     }
 458     if (fPrealignment && ApplyToGeometry()) {stopped = kTRUE; break;}
 459     //
 460     recTitle = fgkRecKeys[ kInitCalSDDFile ];
 461     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
 462       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
 463       fIniSDDRespPath = recOpt;
 464       gSystem->ExpandPathName(fIniSDDRespPath);
 465       fUserProvided |= kSameInitSDDRespBit;
 466       AliInfo(Form("Configuration sets Production SDD Response to %s",fIniSDDRespPath.Data()));
 467     }
 468     if (LoadSDDResponse(fIniSDDRespPath, fIniRespSDD) ) {stopped = kTRUE; break;}
 469     //
 470     //
 471     recTitle = fgkRecKeys[ kInitCorrMapSDDFile ];
 472     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
 473       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
 474       fIniSDDCorrMapPath = recOpt;
 475       gSystem->ExpandPathName(fIniSDDCorrMapPath);
 476       fUserProvided |= kSameInitSDDCorrMapBit;
 477       AliInfo(Form("Configuration sets Production SDD Correction Map to %s",fIniSDDCorrMapPath.Data()));
 478     }
 479     if (LoadSDDCorrMap(fIniSDDCorrMapPath, fIniCorrMapSDD) ) {stopped = kTRUE; break;}
 480     //
 481     recTitle = fgkRecKeys[kPreCalSDDFile];
 482     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
 483       if (!recOpt.IsNull()) {
 484         for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
 485         fPreCalSDDRespPath = recOpt;
 486         gSystem->ExpandPathName(fPreCalSDDRespPath);
 487       }
 488       else if (!fIniSDDRespPath.IsNull()) {
 489         AliInfo("PreCalibration SDD response keyword is present but empty, will set to Init SDD repsonse");
 490         fPreCalSDDRespPath = fIniSDDRespPath;
 491       }
 492       AliInfo(Form("Configuration sets PreCalibration SDD Response to %s",fPreCalSDDRespPath.Data()));
 493     }
 494     //
 495     if (LoadSDDResponse(fPreCalSDDRespPath, fPreRespSDD) ) {stopped = kTRUE; break;}
 496     //
 497     recTitle = fgkRecKeys[kPreCorrMapSDDFile];
 498     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
 499       if (!recOpt.IsNull()) {
 500         for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
 501         fPreSDDCorrMapPath = recOpt;
 502         gSystem->ExpandPathName(fPreSDDCorrMapPath);
 503       }
 504       else if (!fIniSDDCorrMapPath.IsNull()) {
 505         AliInfo("PreCalibration SDD Correction Map keyword is present but empty, will set to Init SDD Correction Map");
 506         fPreSDDCorrMapPath = fIniSDDCorrMapPath;
 507       }
 508       AliInfo(Form("Configuration sets PreCalibration SDD Correction Map to %s",fPreSDDCorrMapPath.Data()));
 509     }
 510     //
 511     if (LoadSDDCorrMap(fPreSDDCorrMapPath, fPreCorrMapSDD) ) {stopped = kTRUE; break;}
 512     //    //
 513     recTitle = fgkRecKeys[ kInitVDriftSDDFile ];
 514     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
 515       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
 516       fIniSDDVDriftPath = recOpt;
 517       gSystem->ExpandPathName(fIniSDDVDriftPath);
 518       fUserProvided |= kSameInitSDDVDriftBit;
 519       AliInfo(Form("Configuration sets Production SDD VDrift to %s",fIniSDDVDriftPath.Data()));
 520     }
 521     if (LoadSDDVDrift(fIniSDDVDriftPath, fIniVDriftSDD) ) {stopped = kTRUE; break;}
 522     //
 523     recTitle = fgkRecKeys[ kPreVDriftSDDFile ];
 524     if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
 525       for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
 526       fPreSDDVDriftPath = recOpt;
 527       gSystem->ExpandPathName(fPreSDDVDriftPath);
 528       AliInfo(Form("Configuration sets PreCalibration SDD VDrift to %s",fPreSDDVDriftPath.Data()));
 529       if (LoadSDDVDrift(fPreSDDVDriftPath, fPreVDriftSDD) ) {stopped = kTRUE; break;}
 530     }
 531     //
 532     recTitle = fgkRecKeys[ kGlobalDeltas ];
 533     if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) SetUseGlobalDelta(kTRUE);
 534     //
 535     recTitle = fgkRecKeys[ kUseDiamond ];
 536     if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
 537       if (!GetUseGlobalDelta()) {
 538         AliError("Diamond constraint is supported only for Global Frame mode");
 539         stopped = kTRUE;
 540         break;
 541       }
 542       fUseDiamond = kTRUE;
 543       if (!recOpt.IsNull()) {
 544         fDiamondPath = recOpt;
 545         gSystem->ExpandPathName(fDiamondPath);
 546         fUserProvided |= kSameDiamondBit;
 547         AliInfo(Form("Configuration sets Diamond constraint to %s",fDiamondPath.Data()));
 548       }
 549     }
 550     //
 551     recTitle = fgkRecKeys[ kUseVertex ];
 552     if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
 553       if (!GetUseGlobalDelta()) {
 554         AliError("Vertex constraint is supported only for Global Frame mode");
 555         stopped = kTRUE;
 556         break;
 557       }
 558       fUseVertex = kTRUE;
 559       if (fUseDiamond) {
 560         AliError("Cannot use Vertex and Diamond constraints at the same time");
 561         stopped = kTRUE;
 562         break;
 563       }
 564       AliInfo("Will use Vertex constraint when available");
 565     }
 566     // =========== 2: see if there are local gaussian constraints defined =====================
 567     //            Note that they should be loaded before the modules declaration
 568     //
 569     recTitle = fgkRecKeys[ kConstrLocal ];
 570     while( (recArr=GetConfigRecord(pfc,recTitle,recOpt,0)) ) {
 571       nrecElems = recArr->GetLast()+1;
 572       if (recOpt.IsFloat()) {stopped = kTRUE; break;} // wrong name
 573       if (GetConstraint(recOpt.Data())) {
 574         AliError(Form("Existing constraint %s repeated",recOpt.Data()));
 575         stopped = kTRUE; break;
 576       }
 577       recExt = recArr->At(2)->GetName();
 578       if (!recExt.IsFloat()) {stopped = kTRUE; break;}
 579       double val = recExt.Atof();
 580       recExt = recArr->At(3)->GetName();
 581       if (!recExt.IsFloat()) {stopped = kTRUE; break;}
 582       double err = recExt.Atof();
 583       int nwgh = nrecElems - 4;
 584       double *wgh = new double[nwgh];
 585       for (nwgh=0,irec=4;irec<nrecElems;irec++) {
 586         recExt = recArr->At(irec)->GetName();
 587         if (!recExt.IsFloat()) {stopped = kTRUE; break;}
 588         wgh[nwgh++] = recExt.Atof();
 589       }
 590       if (stopped) {delete[] wgh; break;}
 591       //
 592       ConstrainLocal(recOpt.Data(),wgh,nwgh,val,err);
 593       delete[] wgh;
 594       //
 595     } // end while for loop over local constraints
 596     if (stopped) break;
 597     //
 598     // =========== 3: now read modules to align ===================================
 599     //
 600     rewind(pfc);
 601     // create fixed modules
 602     for (int j=0; j<fNSuperModules; j++) {
 603       AliITSAlignMille2Module* proto = GetSuperModule(j);
 604       if (!proto->IsAlignable()) continue;
 605       AliITSAlignMille2Module* mod = new AliITSAlignMille2Module(*proto);
 606       // the matrix might be updated in case some prealignment was applied, check
 607       TGeoHMatrix* mup = AliGeomManager::GetMatrix(mod->GetName());
 608       if (mup) *(mod->GetMatrix()) = *mup;
 609       fMilleModule.AddAtAndExpand(mod,fNModules);
 610       mod->SetGeomParamsGlobal(fUseGlobalDelta);
 611       mod->SetUniqueID(fNModules++);
 612       mod->SetNotInConf(kTRUE);
 613     }
 614     CreateVertexModule();
 615     //
 616     while( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0)) ) {
 617       if (!(recTitle==fgkRecKeys[ kModVolID ] || recTitle==fgkRecKeys[ kModIndex ])) continue;
 618       // Expected format: MODULE id tolX tolY tolZ tolPsi tolTh tolPhi [[sigX sigY sigZ]  extra params]
 619       // where tol* is the tolerance (sigma) for given DOF. 0 means fixed
 620       // sig* is the scaling parameters for the errors of the clusters of this module
 621       // extra params are defined for specific modules, e.g. t0 and vdrift corrections of SDD
 622       //
 623       nrecElems = recArr->GetLast()+1;
 624       if (nrecElems<2 || !recOpt.IsDigit()) {stopped = kTRUE; break;}
 625       int idx = recOpt.Atoi();
 626       UShort_t voluid =  (idx<=kMaxITSSensID) ? GetModuleVolumeID(idx) : idx;
 627       AliITSAlignMille2Module* mod = 0;
 628       //
 629       if (voluid>=kMinITSSupeModuleID) { // custom supermodule
 630         mod = GetMilleModuleByVID(voluid);
 631         if (!mod) { // need to create
 632           for (int j=0; j<fNSuperModules; j++) {
 633             if (voluid==GetSuperModule(j)->GetVolumeID()) {
 634               mod = new AliITSAlignMille2Module(*GetSuperModule(j));
 635               // the matrix might be updated in case some prealignment was applied, check
 636               TGeoHMatrix* mup = AliGeomManager::GetMatrix(mod->GetName());
 637               if (mup) *(mod->GetMatrix()) = *mup;
 638               fMilleModule.AddAtAndExpand(mod,fNModules);
 639               mod->SetGeomParamsGlobal(fUseGlobalDelta);
 640               mod->SetUniqueID(fNModules++);
 641               break;
 642             }
 643           }
 644         }
 645         mod->SetNotInConf(kFALSE);
 646       }
 647       else if (idx<=kMaxITSSensVID) {
 648         mod = new AliITSAlignMille2Module(voluid);
 649         fMilleModule.AddAtAndExpand(mod,fNModules);
 650         mod->SetGeomParamsGlobal(fUseGlobalDelta);
 651         mod->SetUniqueID(fNModules++);
 652       }
 653       if (!mod) {stopped = kTRUE; break;}  // bad volid
 654       //
 655       // geometry variation settings
 656       for (int i=0;i<AliITSAlignMille2Module::kMaxParGeom;i++) {
 657         irec = i+2;
 658         if (irec >= nrecElems) break;
 659         recExt = recArr->At(irec)->GetName();
 660         if (!recExt.IsFloat()) {stopped = kTRUE; break;}
 661         mod->SetFreeDOF(i, recExt.Atof() );
 662       }
 663       if (stopped) break;
 664       //
 665       // scaling factors for cluster errors
 666       // first set default ones
 667       for (int i=0;i<3;i++) mod->SetSigmaFactor(i, fSigmaFactor[i]);
 668       for (int i=0;i<3;i++) {
 669         irec = i+8;
 670         if (irec >= nrecElems) break;
 671         recExt = recArr->At(irec)->GetName();
 672         if (!recExt.IsFloat()) {stopped = kTRUE; break;}
 673         mod->SetSigmaFactor(i, recExt.Atof() );
 674       }
 675       if (stopped) break;
 676       //
 677       // now comes special detectors treatment
 678       if (mod->IsSDD()) {
 679         double vl = 0;
 680         if (nrecElems>11) {
 681           recExt = recArr->At(11)->GetName();
 682           if (recExt.IsFloat()) vl = recExt.Atof();
 683           else {stopped = kTRUE; break;}
 684           irec = 11;
 685         }
 686         mod->SetFreeDOF(AliITSAlignMille2Module::kDOFT0,vl);
 687         //
 688         Bool_t cstLR = kFALSE;
 689         for (int lr=0;lr<2;lr++) { // left right side vdrift corrections
 690           vl = 0;
 691           if (nrecElems>12+lr) {
 692             recExt = recArr->At(12+lr)->GetName();
 693             if (recExt.IsFloat()) vl = recExt.Atof();
 694             else {stopped = kTRUE; break;}
 695             irec = 12+lr;
 696           }
 697           mod->SetFreeDOF(lr==0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR,vl);
 698           if (lr==1 && vl>=10) cstLR = kTRUE;  // the right side should be constrained to left one
 699         }
 700         if (cstLR) mod->SetVDriftLRSame();
 701       }
 702       //
 703       mod->EvaluateDOF();
 704       //
 705       // now check if there are local constraints on this module
 706       for (++irec;irec<nrecElems;irec++) {
 707         recExt = recArr->At(irec)->GetName();
 708         if (recExt.IsFloat()) {stopped=kTRUE;break;}
 709         AliITSAlignMille2ConstrArray* cstr = (AliITSAlignMille2ConstrArray*)GetConstraint(recExt.Data());
 710         if (!cstr) {
 711           AliInfo(Form("No Local constraint %s was declared",recExt.Data()));
 712           stopped=kTRUE;
 713           break;
 714         }
 715         cstr->AddModule(mod);
 716       }
 717       if (stopped) break;
 718     } // end while for loop over modules
 719     if (stopped) break;
 720     //
 721     if (fNModules==0) {AliError("Failed to find any MODULE"); stopped = kTRUE; break;}
 722     BuildHierarchy();  // preprocess loaded modules
 723     //
 724     // =========== 4: the rest may come in arbitrary order =======================================
 725     rewind(pfc);
 726     while ( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0))!=0 ) {
 727       //
 728       nrecElems = recArr->GetLast()+1;
 729       //
 730       // some simple flags -----------------------------------------------------------------------
 731       //
 732       if      (recTitle == fgkRecKeys[ kPseudoParents ])  SetAllowPseudoParents(kTRUE);
 733       //
 734       // some optional parameters ----------------------------------------------------------------
 735       else if (recTitle == fgkRecKeys[ kTrackFitMethod ]) {
 736         if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
 737         SetInitTrackParamsMeth(recOpt.Atoi());
 738       }
 739       //
 740       else if (recTitle == fgkRecKeys[ kMinPntTrack ]) {
 741         if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
 742         fMinNPtsPerTrack = recOpt.Atoi();
 743       }
 744       //
 745       else if (recTitle == fgkRecKeys[ kNStDev ]) {
 746         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
 747         fNStdDev = (Int_t)recOpt.Atof();
 748       }
 749       //
 750       else if (recTitle == fgkRecKeys[ kResCutInit  ]) {
 751         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
 752         fResCutInitial = recOpt.Atof();
 753       }
 754       //
 755       else if (recTitle == fgkRecKeys[ kResCutOther ]) {
 756         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
 757         fResCut = recOpt.Atof();
 758       }
 759       //
 760       else if (recTitle == fgkRecKeys[ kLocalSigFactor ]) { //-------------------------
 761         for (irec=0;irec<3;irec++) if (nrecElems>irec+1) {
 762             fSigmaFactor[irec] = ((TObjString*)recArr->At(irec+1))->GetString().Atof();
 763             if (fSigmaFactor[irec]<=0.) stopped = kTRUE;
 764           }
 765         if (stopped) break;
 766       }
 767       //
 768       else if (recTitle == fgkRecKeys[ kStartFactor ]) {        //-------------------------
 769         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
 770         fStartFac = recOpt.Atof();
 771       }
 772       //
 773       else if (recTitle == fgkRecKeys[ kFinalFactor ]) {        //-------------------------
 774         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
 775         fFinalFac = recOpt.Atof();
 776       }
 777       //
 778       // pepo2708909
 779       else if (recTitle == fgkRecKeys[ kExtraClustersMode ]) {        //-------------------------
 780         if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
 781         fExtraClustersMode = recOpt.Atoi();
 782       }
 783       // endpepo270809
 784       //
 785       else if (recTitle == fgkRecKeys[ kBField ]) {         //-------------------------
 786         if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
 787         SetBField( recOpt.Atof() );
 788       }
 789       //
 790       else if (recTitle == fgkRecKeys[ kSDDVDCorrMult ]) {         //-------------------------
 791         SetSDDVDCorrMult( recOpt.IsNull() || (recOpt.IsFloat() && (recOpt.Atof())>-0.5) );
 792       }
 793       //
 794       else if (recTitle == fgkRecKeys[ kWeightPt ]) {         //-------------------------
 795         double wgh = 1;
 796         if (!recOpt.IsNull()) {
 797           if (!recOpt.IsFloat()) {stopped = kTRUE; break;}
 798           else wgh = recOpt.Atof();
 799         }
 800         SetWeightPt(wgh);
 801       }
 802       //
 803       else if (recTitle == fgkRecKeys[ kSparseMatrix ]) {   // matrix solver type
 804         //
 805         AliMillePede2::SetGlobalMatSparse(kTRUE);
 806         if (recOpt.IsNull()) continue;
 807         // solver type and settings
 808         if      (recOpt == "MINRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolMinRes );
 809         else if (recOpt == "FGMRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolFGMRes );
 810         else {stopped = kTRUE; break;}
 811         //
 812         if (nrecElems>=3) { // preconditioner type
 813           recExt = recArr->At(2)->GetName();
 814           if (!recExt.IsDigit()) {stopped = kTRUE; break;}
 815           AliMillePede2::SetMinResPrecondType( recExt.Atoi() );
 816         }
 817         //
 818         if (nrecElems>=4) { // tolerance
 819           recExt = recArr->At(3)->GetName();
 820           if (!recExt.IsFloat()) {stopped = kTRUE; break;}
 821           AliMillePede2::SetMinResTol( recExt.Atof() );
 822         }
 823         //
 824         if (nrecElems>=5) { // maxIter
 825           recExt = recArr->At(4)->GetName();
 826           if (!recExt.IsDigit()) {stopped = kTRUE; break;}
 827           AliMillePede2::SetMinResMaxIter( recExt.Atoi() );
 828         }
 829       }
 830       //
 831       else if (recTitle == fgkRecKeys[ kRequirePoint ]) {       //-------------------------
 832         // syntax:   REQUIRE_POINT where ndet updw nreqpts
 833         //    where = LAYER or DETECTOR
 834         //    ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
 835         //    updw = 1 for Y>0, -1 for Y<0, 0 if not specified
 836         //    nreqpts = minimum number of points of that type
 837         if (nrecElems>=5) {
 838           recOpt.ToUpper();
 839           int lr = ((TObjString*)recArr->At(2))->GetString().Atoi() - 1;
 840           int hb = ((TObjString*)recArr->At(3))->GetString().Atoi();
 841           int np = ((TObjString*)recArr->At(4))->GetString().Atoi();
 842           //
 843           int rtp = -1; // use for run type
 844           if (nrecElems>5) {
 845             TString tpstr = ((TObjString*)recArr->At(5))->GetString();
 846             if ( tpstr.Contains("cosmics",TString::kIgnoreCase) ) rtp = kCosmics;
 847             else if ( tpstr.Contains("collision",TString::kIgnoreCase) ) rtp = kCollision;
 848             else {stopped = kTRUE; break;}
 849           }
 850           //
 851           int tpmn= rtp<0 ? 0 : rtp;
 852           int tpmx= rtp<0 ? kNDataType-1 : rtp;
 853           for (int itp=tpmn;itp<=tpmx;itp++) {
 854             fRequirePoints[itp]=kTRUE;
 855             if (recOpt == "LAYER") {
 856               if (lr<0 || lr>5) {stopped = kTRUE; break;}
 857               if (hb>0) fNReqLayUp[itp][lr]=np;
 858               else if (hb<0) fNReqLayDown[itp][lr]=np;
 859               else fNReqLay[itp][lr]=np;
 860             }
 861             else if (recOpt == "DETECTOR") {
 862               if (lr<0 || lr>2) {stopped = kTRUE; break;}
 863               if (hb>0) fNReqDetUp[itp][lr]=np;
 864               else if (hb<0) fNReqDetDown[itp][lr]=np;
 865               else fNReqDet[itp][lr]=np;
 866             }
 867             else {stopped = kTRUE; break;}
 868           }
 869           if (stopped) break;
 870         }
 871         else {stopped = kTRUE; break;}
 872       }
 873       //
 874       // global constraints on the subunits/orphans
 875       else if (recTitle == fgkRecKeys[ kConstrOrphans ]) {    //------------------------
 876         // expect CONSTRAINT_ORPHANS MEAN/MEDIAN Value parID0 ... parID1 ...
 877         if (nrecElems<4) {stopped = kTRUE; break;}
 878         recExt = recArr->At(2)->GetName();
 879         if (!recExt.IsFloat()) {stopped = kTRUE; break;}
 880         double val = recExt.Atof();
 881         UInt_t pattern = 0;
 882         for (irec=3;irec<nrecElems;irec++) { // read params to constraint
 883           recExt = recArr->At(irec)->GetName();
 884           if (!recExt.IsDigit()) {stopped = kTRUE; break;}
 885           pattern |= 0x1 << recExt.Atoi();
 886         }
 887         if (stopped) break;
 888         if      (recOpt == "MEAN")   ConstrainOrphansMean(val,pattern);
 889         else if (recOpt == "MEDIAN") ConstrainOrphansMedian(val,pattern);
 890         else {stopped = kTRUE; break;}
 891       }
 892       //
 893       else if (recTitle == fgkRecKeys[ kConstrSubunits ]) {    //------------------------
 894         // expect CONSTRAINT_SUBUNITS MEAN/MEDIAN Value parID0 ... parID1 ... VolID1 ... VolIDn - VolIDm
 895         if (nrecElems<5) {stopped = kTRUE; break;}
 896         recExt = recArr->At(2)->GetName();
 897         if (!recExt.IsFloat()) {stopped = kTRUE; break;}
 898         double val = recExt.Atof();
 899         UInt_t pattern = 0;
 900         for (irec=3;irec<nrecElems;irec++) { // read params to constraint
 901           recExt = recArr->At(irec)->GetName();
 902           if (!recExt.IsDigit()) {stopped = kTRUE; break;}
 903           int parid = recExt.Atoi();
 904           if (parid<kMaxITSSensID) pattern |= 0x1 << recExt.Atoi();
 905           else break;           // list of params is over
 906         }
 907         if (stopped) break;
 908         //
 909         Bool_t meanC;
 910         if      (recOpt == "MEAN")   meanC = kTRUE;
 911         else if (recOpt == "MEDIAN") meanC = kFALSE;
 912         else    {stopped = kTRUE; break;}
 913         //
 914         int curID = -1;
 915         int rangeStart = -1;
 916         for (;irec<nrecElems;irec++) { // read modules to apply this constraint
 917           recExt = recArr->At(irec)->GetName();
 918           if (recExt == "-") {rangeStart = curID; continue;}  // range is requested
 919           else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
 920           else curID = recExt.Atoi();
 921           //
 922           if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
 923           // this was a range start or single
 924           int start;
 925           if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
 926           else start = curID;  // create constraint either for single module (or 1st in the range)
 927           for (int id=start;id<=curID;id++) {
 928             int id0 = IsVIDDefined(id);
 929             if (id0<0) {AliDebug(3,Form("Undefined module %d requested in the SubUnits constraint, skipping",id)); continue;}
 930             if (meanC) ConstrainModuleSubUnitsMean(id0,val,pattern);
 931             else       ConstrainModuleSubUnitsMedian(id0,val,pattern);
 932           }
 933         }
 934         if (rangeStart>=0) stopped = kTRUE; // unfinished range
 935         if (stopped) break;
 936       }
 937       //
 938       // association of modules with local constraints
 939       else if (recTitle == fgkRecKeys[ kApplyConstr ]) {            //------------------------
 940         // expect APPLY_CONSTRAINT NAME [NAME1...] [VolID1 ... VolIDn - VolIDm]
 941         if (nrecElems<3) {stopped = kTRUE; break;}
 942         int nmID0=-1,nmID1=-1;
 943         for (irec=1;irec<nrecElems;irec++) { // find the range of constraint names
 944           recExt = recArr->At(irec)->GetName();
 945           if (recExt.IsFloat()) break;
 946           // check if such a constraint was declared
 947           if (!GetConstraint(recExt.Data())) {
 948             AliInfo(Form("No Local constraint %s was declared",recExt.Data()));
 949             stopped=kTRUE;
 950             break;
 951           }
 952           if (nmID0<0) nmID0 = irec;
 953           nmID1 = irec;
 954         }
 955         if (stopped) break;
 956         //
 957         if (irec>=nrecElems) {stopped = kTRUE; break;} // no modules provided
 958         //
 959         // now read the list of modules to constrain
 960         int curID = -1;
 961         int rangeStart = -1;
 962         for (;irec<nrecElems;irec++) { // read modules to apply this constraint
 963           recExt = recArr->At(irec)->GetName();
 964           if (recExt == "-") {rangeStart = curID; continue;}  // range is requested
 965           else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
 966           else curID = recExt.Atoi();
 967           //
 968           if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
 969           //
 970           // this was a range start or single
 971           int start;
 972           if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
 973           else start = curID;  // create constraint either for single module (or 1st in the range)
 974           for (int id=start;id<=curID;id++) {
 975             AliITSAlignMille2Module *md = GetMilleModuleByVID(id);
 976             if (!md) {AliDebug(3,Form("Undefined module %d requested in the Local constraint, skipping",id)); continue;}
 977             for (int nmid=nmID0;nmid<=nmID1;nmid++)
 978               ((AliITSAlignMille2ConstrArray*)GetConstraint(recArr->At(nmid)->GetName()))->AddModule(md);
 979           }
 980         }
 981         if (rangeStart>=0) stopped = kTRUE; // unfinished range
 982         if (stopped) break;
 983       }
 984       //
 985       // request of the same T0 for group of SDD modules
 986       else if (recTitle == fgkRecKeys[ kSameSDDT0 ]) {            //------------------------
 987         // expect SET_SAME_SDDT0 [SensID1 ... SensIDn - SensIDm]
 988         if (nrecElems<3) {stopped = kTRUE; break;}
 989         //
 990         // now read the list of modules to constrain
 991         int curID = -1;
 992         int rangeStart = -1;
 993         AliITSAlignMille2ConstrArray *cstrT0 = new AliITSAlignMille2ConstrArray("SDDT0",0,0,0,0);
 994         int naddM = 0;
 995         cstrT0->SetPattern(BIT(AliITSAlignMille2Module::kDOFT0));
 996         for (irec=1;irec<nrecElems;irec++) { // read modules to apply this constraint
 997           recExt = recArr->At(irec)->GetName();
 998           if (recExt == "-") {rangeStart = curID; continue;}  // range is requested
 999           else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
1000           else curID = recExt.Atoi();
1001           //
1002           if (curID<kSDDoffsID || curID>=kSDDoffsID+kNSDDmod) {stopped = kTRUE; break;}
1003           //
1004           // this was a range start or single
1005           int start;
1006           if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
1007           else start = curID;  // create constraint either for single module (or 1st in the range)
1008           for (int id=start;id<=curID;id++) {
1009             int vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
1010             if (vid<=1) {AliDebug(3,Form("Undefined module index %d requested in the SAME_SDDT0 constraint, skipping",id)); continue;}
1011             AliITSAlignMille2Module *md = GetMilleModuleByVID(vid);
1012             if (!md) {AliDebug(3,Form("Undefined module %d requested in the Local constraint, skipping",id)); continue;}
1013             cstrT0->AddModule(md,kFALSE);
1014             naddM++;
1015           }
1016         }
1017         if (rangeStart>=0) stopped = kTRUE; // unfinished range
1018         if (stopped) break;
1019         if (naddM<2) delete cstrT0;
1020         else {
1021           cstrT0->SetConstraintID(GetNConstraints());
1022           fConstraints.Add(cstrT0);
1023         }
1024       }
1025       //
1026       // Do we use new local Y errors?
1027       else if (recTitle == fgkRecKeys[ kUseLocalYErr ]) {
1028         // expect SET_TPAFITTER
1029         fUseLocalYErr = kTRUE;
1030       }
1031       //
1032       else if (recTitle == fgkRecKeys[ kMinPointsSens ]) {         //-------------------------
1033         if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
1034         SetMinPointsPerSensor( recOpt.Atoi() );
1035       }
1036       //
1037       else if (recTitle == fgkRecKeys[ kOCDBSpecificPath ]) {         //-------------------------
1038         if (recOpt.IsNull() || nrecElems<3 ) {stopped = kTRUE; break;}
1039         AliCDBManager::Instance()->SetSpecificStorage(recOpt.Data(), gSystem->ExpandPathName(recArr->At(2)->GetName()));
1040         AliInfo(Form("Configuration sets OCDB specific storage %s to %s",recOpt.Data(),recArr->At(2)->GetName()));
1041         TObjString *pths = (TObjString*)AliCDBManager::Instance()->GetStorageMap()->GetValue(recOpt.Data());
1042         if (!pths) { stopped = kTRUE; break; }
1043         pths->SetUniqueID(1); // mark as set by user
1044       }
1045       //
1046       else if (recTitle == fgkRecKeys[ kCorrectDiamond ] && fUseDiamond) {         //-------------------------
1047         if (nrecElems<4) {stopped = kTRUE; break;}
1048         for (int i=0;i<3;i++) fCorrDiamond[i] = ((TObjString*)recArr->At(i+1))->GetString().Atof();
1049         AliInfo(Form("Correction %+.4f %+.4f %+.4f will be applied to diamond",fCorrDiamond[0],fCorrDiamond[1],fCorrDiamond[2]));
1050       }
1051       //
1052       else continue; // already processed record
1053       //
1054     } // end of while loop 4 over the various params
1055     //
1056     break;
1057   } // end of while(1) loop
1058   //
1059   fclose(pfc);
1060   if (!fDiamondPath.IsNull() && IsDiamondUsed() && LoadDiamond(fDiamondPath) ) stopped = kTRUE;
1061   if (stopped) {
1062     AliError(Form("Failed on record %s %s ...\n",recTitle.Data(),recOpt.Data()));
1063     return -1;
1064   }
1065   //
1066   if (CacheMatricesCurr()) return -1;
1067   SetUseLocalYErrors(fUseLocalYErr); // YErr used only with TPAFitter
1068   fSegmentationSDD = new AliITSsegmentationSDD();
1069   //
1070   fIsConfigured = kTRUE;
1071   return 0;
1072 }
1073
1074 //________________________________________________________________________________________________________
1075 void AliITSAlignMille2::BuildHierarchy()
1076 {
1077   // build the hieararhy of the modules to align
1078   //
1079   if (!GetUseGlobalDelta() && PseudoParentsAllowed()) {
1080     AliInfo("PseudoParents mode is allowed only when the deltas are global\n"
1081             "Since Deltas are local, switching to NoPseudoParents");
1082     SetAllowPseudoParents(kFALSE);
1083   }
1084   // set parent/child relationship for modules to align
1085   AliInfo("Setting parent/child relationships\n");
1086   //
1087   // 1) child -> parent reference
1088   for (int ipar=0;ipar<fNModules;ipar++) {
1089     AliITSAlignMille2Module* parent = GetMilleModule(ipar);
1090     if (parent->IsSensor()) continue; // sensor cannot be a parent
1091     //
1092     for (int icld=0;icld<fNModules;icld++) {
1093       if (icld==ipar) continue;
1094       AliITSAlignMille2Module* child = GetMilleModule(icld);
1095       if (!child->BelongsTo(parent)) continue;
1096       // child cannot have more sensors than the parent
1097       if (child->GetNSensitiveVolumes() > parent->GetNSensitiveVolumes()) continue;
1098       //
1099       AliITSAlignMille2Module* parOld = child->GetParent();
1100       // is this parent candidate closer than the old parent ?
1101       if (parOld && parOld->GetNSensitiveVolumes()<parent->GetNSensitiveVolumes()) continue; // parOld is closer
1102       child->SetParent(parent);
1103     }
1104     //
1105   }
1106   //
1107   // add parent -> children reference
1108   for (int icld=0;icld<fNModules;icld++) {
1109     AliITSAlignMille2Module* child = GetMilleModule(icld);
1110     AliITSAlignMille2Module* parent = child->GetParent();
1111     if (parent) parent->AddChild(child);
1112   }
1113   //
1114   // reorder the modules in such a way that parents come first
1115   for (int icld=0;icld<fNModules;icld++) {
1116     AliITSAlignMille2Module* child  = GetMilleModule(icld);
1117     AliITSAlignMille2Module* parent;
1118     while ( (parent=child->GetParent()) &&  (parent->GetUniqueID()>child->GetUniqueID()) ) {
1119       // swap
1120       fMilleModule[icld] = parent;
1121       fMilleModule[parent->GetUniqueID()] = child;
1122       child->SetUniqueID(parent->GetUniqueID());
1123       parent->SetUniqueID(icld);
1124       child = parent;
1125     }
1126     //
1127   }
1128   //
1129   // Go over the child->parent chain and mark modules with explicitly provided sensors.
1130   // If the sensors of the unit are explicitly declared, all undeclared sensors are
1131   // suppresed in this unit.
1132   for (int icld=fNModules;icld--;) {
1133     AliITSAlignMille2Module* child = GetMilleModule(icld);
1134     AliITSAlignMille2Module* parent = child->GetParent();
1135     if (!parent) continue;
1136     //
1137     // check if this parent was already processed
1138     if (!parent->AreSensorsProvided()) {
1139       parent->DelSensitiveVolumes();
1140       parent->SetSensorsProvided(kTRUE);
1141     }
1142     // reattach sensors to parent
1143     for (int isc=child->GetNSensitiveVolumes();isc--;) {
1144       UShort_t senVID = child->GetSensVolVolumeID(isc);
1145       if (!parent->IsIn(senVID)) parent->AddSensitiveVolume(senVID);
1146     }
1147   }
1148   //
1149 }
1150
1151 // pepo
1152 //________________________________________________________________________________________________________
1153 void AliITSAlignMille2::SetCurrentModule(Int_t id)
1154 {
1155   // set the current supermodule
1156   // new meaning
1157   if (fMilleVersion>=2) {
1158     fCurrentModule = GetMilleModule(id);
1159     return;
1160   }
1161   // old meaning
1162   if (fMilleVersion<=1) {
1163     Int_t index=id;
1164     /// set as current the SuperModule that contains the 'index' sens.vol.
1165     if (index<0 || index>2197) {
1166       AliInfo("index does not correspond to a sensitive volume!");
1167       return;
1168     }
1169     UShort_t voluid=AliITSAlignMille2Module::GetVolumeIDFromIndex(index);
1170     Int_t k=IsContained(voluid);
1171     if (k>=0){
1172       fCurrentSensID = index;
1173       fCluster.SetVolumeID(voluid);
1174       fCluster.SetXYZ(0,0,0);
1175       InitModuleParams();
1176     }
1177     else
1178       AliInfo(Form("module %d not defined\n",index));
1179   }
1180 }
1181 // endpepo
1182 //________________________________________________________________________________________________________
1183 void AliITSAlignMille2::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts,Int_t runtype)
1184 {
1185   // set minimum number of points in specific detector or layer
1186   // where = LAYER or DETECTOR
1187   // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
1188   // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
1189   // nreqpts = minimum number of points of that type
1190   ndet--;
1191   int tpmn= runtype<0 ? 0 : runtype;
1192   int tpmx= runtype<0 ? kNDataType-1 : runtype;
1193   //
1194   for (int itp=tpmn;itp<=tpmx;itp++) {
1195     fRequirePoints[itp]=kTRUE;
1196     if (strstr(where,"LAYER")) {
1197       if (ndet<0 || ndet>5) return;
1198       if (updw>0) fNReqLayUp[itp][ndet]=nreqpts;
1199       else if (updw<0) fNReqLayDown[itp][ndet]=nreqpts;
1200       else fNReqLay[itp][ndet]=nreqpts;
1201     }
1202     else if (strstr(where,"DETECTOR")) {
1203       if (ndet<0 || ndet>2) return;
1204       if (updw>0) fNReqDetUp[itp][ndet]=nreqpts;
1205       else if (updw<0) fNReqDetDown[itp][ndet]=nreqpts;
1206       else fNReqDet[itp][ndet]=nreqpts;
1207     }
1208   }
1209 }
1210
1211 //________________________________________________________________________________________________________
1212 Int_t AliITSAlignMille2::GetModuleIndex(const Char_t *symname)
1213 {
1214   /// index from symname
1215   if (!symname) return -1;
1216   for (Int_t i=0;i<=kMaxITSSensID; i++) {
1217     if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i;
1218   }
1219   return -1;
1220 }
1221
1222 //________________________________________________________________________________________________________
1223 Int_t AliITSAlignMille2::GetModuleIndex(UShort_t voluid)
1224 {
1225   /// index from volume ID
1226   AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
1227   if (lay<1|| lay>6) return -1;
1228   Int_t idx=Int_t(voluid)-2048*lay;
1229   if (idx>=AliGeomManager::LayerSize(lay)) return -1;
1230   for (Int_t ilay=1; ilay<lay; ilay++)
1231     idx += AliGeomManager::LayerSize(ilay);
1232   return idx;
1233 }
1234
1235 //________________________________________________________________________________________________________
1236 UShort_t AliITSAlignMille2::GetModuleVolumeID(const Char_t *symname)
1237 {
1238   /// volume ID from symname
1239   /// works for sensitive volumes only
1240   if (!symname) return 0;
1241
1242   for (UShort_t voluid=2000; voluid<13300; voluid++) {
1243     Int_t modId;
1244     AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
1245     if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
1246       if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
1247     }
1248   }
1249
1250   return 0;
1251 }
1252
1253 //________________________________________________________________________________________________________
1254 UShort_t AliITSAlignMille2::GetModuleVolumeID(Int_t index)
1255 {
1256   /// volume ID from index
1257   if (index<0) return 0;
1258   if (index<2198)
1259     return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index));
1260   else {
1261     for (int i=0; i<fNSuperModules; i++) {
1262       if (GetSuperModule(i)->GetIndex()==index) return GetSuperModule(i)->GetVolumeID();
1263     }
1264   }
1265   return 0;
1266 }
1267
1268 //________________________________________________________________________________________________________
1269 Int_t AliITSAlignMille2::LoadGeometry(TString& path)
1270 {
1271   // initialize ideal geometry
1272   AliInfo(Form("Loading ideal geometry %s",path.Data()));
1273   if (path.IsNull()) {
1274     AliError("Path to geometry is not provided");
1275     return -1;
1276   }
1277   //
1278   AliCDBEntry *entry = 0;
1279   TGeoManager *gm = 0;
1280   while(1) {
1281     if (path.BeginsWith("path: ")) { // must load from OCDB
1282       entry = GetCDBEntry(path.Data());
1283       if (!entry) break;
1284       gm = (TGeoManager*) entry->GetObject();
1285       entry->SetObject(NULL);
1286       entry->SetOwner(kTRUE);
1287       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
1288       //      delete cdbId;
1289       //      delete entry;
1290       break;
1291     }
1292     //
1293     if (gSystem->AccessPathName(path.Data())) break;
1294     TFile* precf = TFile::Open(path.Data());
1295     if (precf->FindKey("ALICE")) gm = (TGeoManager*)precf->Get("ALICE");
1296     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
1297       gm = (TGeoManager*) entry->GetObject();
1298       if (gm && gm->InheritsFrom(TGeoManager::Class())) entry->SetObject(NULL);
1299       else gm = 0;
1300       entry->SetObject(NULL);
1301       entry->SetOwner(kTRUE);
1302       delete entry;
1303     }
1304     precf->Close();
1305     delete precf;
1306     break;
1307   }
1308   //
1309   if (!gm) {AliError(Form("Failed to load geometry from %s",path.Data())); return -1;}
1310   AliGeomManager::SetGeometry(gm);
1311   fGeoManager = AliGeomManager::GetGeometry();
1312   if (!fGeoManager) {
1313     AliInfo("Couldn't initialize geometry");
1314     return -1;
1315   }
1316   return 0;
1317 }
1318
1319 //________________________________________________________________________________________________________
1320 Int_t AliITSAlignMille2::SetConstraintWrtRef(const char* reffname)
1321 {
1322   // Load the global deltas from this file. The local gaussian constraints on some modules
1323   // will be defined with respect to the deltas from this reference file, converted to local
1324   // delta format. Note: conversion to local format requires reloading the geometry!
1325   //
1326   AliInfo(Form("Loading reference deltas for local constraints from %s",reffname));
1327   if (!fGeoManager) return -1;
1328   fConstrRefPath = reffname;
1329   if (fConstrRefPath == "IDEAL") { // the reference is the ideal geometry, just create dummy reference array
1330     fConstrRef = new TClonesArray("AliAlignObjParams",1);
1331     return 0;
1332   }
1333   if (LoadDeltas(fConstrRefPath,fConstrRef)) return -1;
1334   //
1335   // we need ideal geometry to convert global deltas to local ones
1336   if (fUsePreAlignment) {
1337     AliError("The call of SetConstraintWrtRef must be done before application of the prealignment");
1338     return -1;
1339   }
1340   //
1341   AliInfo("Converting global reference deltas to local ones");
1342   Int_t nprea = fConstrRef->GetEntriesFast();
1343   for (int ix=0; ix<nprea; ix++) {
1344     AliAlignObjParams *preo=(AliAlignObjParams*) fConstrRef->At(ix);
1345     if (!preo->ApplyToGeometry()) return -1;
1346   }
1347   //
1348   // now convert the global reference deltas to local ones
1349   for (int i=fConstrRef->GetEntriesFast();i--;) {
1350     AliAlignObjParams *preo = (AliAlignObjParams*)fConstrRef->At(i);
1351     TGeoHMatrix * mupd = AliGeomManager::GetMatrix(preo->GetSymName());
1352     if (!mupd) {  // this is not alignable entry, need to look in the supermodules
1353       for (int im=fNSuperModules;im--;) {
1354         AliITSAlignMille2Module* mod = GetSuperModule(im);
1355         if ( strcmp(mod->GetName(), preo->GetSymName()) ) continue;
1356         mupd = mod->GetMatrix();
1357         break;
1358       }
1359       if (!mupd) {
1360         AliError(Form("Failed to find the volume for reference %s",preo->GetSymName()));
1361         return -1;
1362       }
1363     }
1364     TGeoHMatrix preMat;
1365     preo->GetMatrix(preMat);                     //  Delta_Glob
1366     TGeoHMatrix tmpMat    = *mupd;               //  Delta_Glob * Delta_Glob_Par * M
1367     preMat.MultiplyLeft( &tmpMat.Inverse() );    //  M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
1368     tmpMat.MultiplyLeft( &preMat );              //  (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
1369     preo->SetMatrix(tmpMat);     // local corrections
1370   }
1371   //
1372   // we need to reload the geometry spoiled by this reference deltas...
1373   delete fGeoManager;
1374   AliInfo("Reloading target ideal geometry");
1375   return LoadGeometry(fGeometryPath);
1376   //
1377 }
1378
1379 //________________________________________________________________________________________________________
1380 void AliITSAlignMille2::Init()
1381 {
1382   // perform global initialization
1383   //
1384   if (fIsMilleInit) {
1385     AliInfo("Millepede has been already initialized!");
1386     return;
1387   }
1388   // range constraints in such a way that the childs are constrained before their parents
1389   // orphan constraints come last
1390   for (int ic=0;ic<GetNConstraints();ic++) {
1391     for (int ic1=ic+1;ic1<GetNConstraints();ic1++) {
1392       AliITSAlignMille2Constraint *cst0 = GetConstraint(ic);
1393       AliITSAlignMille2Constraint *cst1 = GetConstraint(ic1);
1394       if (cst0->GetModuleID()<cst1->GetModuleID()) {
1395         // swap
1396         fConstraints[ic] = cst1;
1397         fConstraints[ic1] = cst0;
1398       }
1399     }
1400   }
1401   //
1402   if (!GetUseGlobalDelta()) {
1403     AliInfo("ATTENTION: The parameters are defined in the local frame, no check for degeneracy will be done");
1404     for (int imd=fNModules;imd--;) {
1405       AliITSAlignMille2Module* mod = GetMilleModule(imd);
1406       int npar = mod->GetNParTot();
1407       // the parameter may have max 1 free instance, otherwise the equations are underdefined
1408       for (int ipar=0;ipar<npar;ipar++) {
1409         if (!mod->IsFreeDOF(ipar)) continue;
1410         mod->SetParOffset(ipar,fNGlobal++);
1411       }
1412     }
1413   }
1414   else {
1415     // init millepede, decide which parameters are to be fitted explicitly
1416     for (int imd=fNModules;imd--;) {
1417       AliITSAlignMille2Module* mod = GetMilleModule(imd);
1418       if (mod->IsNotInConf()) continue; // dummy module
1419       int npar = mod->GetNParTot();
1420       // the parameter may have max 1 free instance, otherwise the equations are underdefined
1421       for (int ipar=0;ipar<npar;ipar++) {
1422         if (!mod->IsFreeDOF(ipar)) continue;  // fixed
1423         //
1424         int nFreeInstances = 0;
1425         //
1426         AliITSAlignMille2Module* parent = mod;
1427         Bool_t cstMeanMed=kFALSE,cstGauss=kFALSE;
1428         //
1429         Bool_t addToFit = kFALSE;
1430         // the parameter may be ommitted from explicit fit (if PseudoParentsAllowed is true) if
1431         // 1) it is not explicitly constrained or its does not participate in Gaussian constraint
1432         // 2) the same applies to all of its parents
1433         // 3) it has at least 1 unconstrained direct child
1434         while(parent) {
1435           if (parent->IsNotInConf()) {parent = parent->GetParent(); continue;}
1436           if (!parent->IsFreeDOF(ipar)) {parent = parent->GetParent(); continue;}
1437           nFreeInstances++;
1438           if (IsParModConstrained(parent,ipar, cstMeanMed, cstGauss)) nFreeInstances--;
1439           if (cstGauss) addToFit = kTRUE;
1440           parent = parent->GetParent();
1441         }
1442         if (nFreeInstances>1) {
1443           AliError(Form("Parameter#%d of module %s\nhas %d free instances in the "
1444                         "unconstrained parents\nSystem is undefined",ipar,mod->GetName(),nFreeInstances));
1445           exit(1);
1446         }
1447         //
1448         // i) Are PseudoParents allowed?
1449         if (!PseudoParentsAllowed()) addToFit = kTRUE;
1450         // ii) check if this module has no child with such a free parameter. Since the order of this check
1451         // goes from child to parent, by this moment such a parameter must have been already added
1452         else if (!IsParModFamilyVaried(mod,ipar))  addToFit = kTRUE;  // no varied children at all
1453         else if (!IsParFamilyFree(mod,ipar,1))     addToFit = kTRUE;  // no unconstrained direct children
1454         // otherwise the value of this parameter can be extracted from simple contraint and the values of
1455         // the relevant parameters of its children the fit is done. Hence it is not included
1456         if (!addToFit) continue;
1457         //
1458         // shall add this parameter to explicit fit
1459         //      printf("Adding %s %d -> %d\n",mod->GetName(), ipar, fNGlobal);
1460         mod->SetParOffset(ipar,fNGlobal++);
1461       }
1462     }
1463   }
1464   //
1465   AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, kNLocal, fNStdDev));
1466   fGlobalDerivatives = new Double_t[fNGlobal];
1467   memset(fGlobalDerivatives,0,fNGlobal*sizeof(Double_t));
1468   //
1469   fMillepede->InitMille(fNGlobal,kNLocal,fNStdDev,fResCut,fResCutInitial);
1470   fMillepede->SetMinPntValid(fMinPntPerSens);
1471   fIsMilleInit = kTRUE;
1472   //
1473   ResetLocalEquation();
1474   AliInfo("Parameters initialized to zero");
1475   //
1476   /// Fix non free parameters
1477   for (Int_t i=0; i<fNModules; i++) {
1478     AliITSAlignMille2Module* mod = GetMilleModule(i);
1479     for (Int_t j=0; j<mod->GetNParTot(); j++) {
1480       if (mod->GetParOffset(j)<0) continue; // not varied
1481       FixParameter(mod->GetParOffset(j),mod->GetParConstraint(j));
1482       fMillepede->SetParamGrID(i, mod->GetParOffset(j));
1483     }
1484   }
1485   //
1486   ResetCovIScale();
1487   // Set iterations
1488   if (fStartFac>1) fMillepede->SetIterations(fStartFac);
1489   if (fFinalFac>1) fMillepede->SetChi2CutRef(fFinalFac);
1490   fTrackBuff.Expand(24);
1491   //
1492 }
1493
1494 //________________________________________________________________________________________________________
1495 void AliITSAlignMille2::AddConstraint(Double_t *par, Double_t value, Double_t sigma)
1496 {
1497   /// Constrain equation defined by par to value
1498   if (!fIsMilleInit) Init();
1499   fMillepede->SetGlobalConstraint(par, value, sigma);
1500   AliInfo("Adding constraint");
1501 }
1502
1503 //________________________________________________________________________________________________________
1504 void AliITSAlignMille2::InitGlobalParameters(Double_t *par)
1505 {
1506   /// Initialize global parameters with par array
1507   if (!fIsMilleInit) Init();
1508   fMillepede->SetGlobalParameters(par);
1509   AliInfo("Init Global Parameters");
1510 }
1511
1512 //________________________________________________________________________________________________________
1513 void AliITSAlignMille2::FixParameter(Int_t iPar, Double_t value)
1514 {
1515   /// Parameter iPar is encourage to vary in [-value;value].
1516   /// If value == 0, parameter is fixed
1517   if (!fIsMilleInit) {
1518     AliInfo("Millepede has not been initialized!");
1519     return;
1520   }
1521   fMillepede->SetParSigma(iPar, value);
1522   if (IsZero(value)) AliInfo(Form("Parameter %i Fixed", iPar));
1523 }
1524
1525 //________________________________________________________________________________________________________
1526 void AliITSAlignMille2::ResetLocalEquation()
1527 {
1528   /// Reset the derivative vectors
1529   for(int i=kNLocal;i--;)  fLocalDerivatives[i] = 0.0;
1530   memset(fGlobalDerivatives, 0, fNGlobal*sizeof(double) );
1531 }
1532
1533 //________________________________________________________________________________________________________
1534 Int_t AliITSAlignMille2::ApplyToGeometry()
1535 {
1536   // apply prealignment to ideal geometry
1537   Int_t nprea = fPrealignment->GetEntriesFast();
1538   AliInfo(Form("Array of prealignment deltas: %d entries",nprea));
1539   //
1540   for (int ix=0; ix<nprea; ix++) {
1541     AliAlignObjParams *preo=(AliAlignObjParams*) fPrealignment->At(ix);
1542     Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(preo->GetVolUID());
1543     if (index>=0) {
1544       if (index>=fPreAlignQF.GetSize()) fPreAlignQF.Set(index+10);
1545       fPreAlignQF[index] = (int) preo->GetUniqueID()+1;
1546     }
1547     if (!preo->ApplyToGeometry()) {
1548       AliError(Form("Failed on ApplyToGeometry at %s",preo->GetSymName()));
1549       return -1;
1550     }
1551   }
1552   //
1553   fUsePreAlignment = kTRUE;
1554   return 0;
1555 }
1556
1557 //________________________________________________________________________________________________________
1558 Int_t AliITSAlignMille2::GetPreAlignmentQualityFactor(Int_t index) const
1559 {
1560   // quality factors from prealignment
1561   if (!fUsePreAlignment || index<0 || index>=fPreAlignQF.GetSize()) return -1;
1562   return fPreAlignQF[index]-1;
1563 }
1564
1565 //________________________________________________________________________________________________________
1566 AliTrackPointArray *AliITSAlignMille2::PrepareTrack(const AliTrackPointArray *atp)
1567 {
1568   /// create a new AliTrackPointArray keeping only defined modules
1569   /// move points according to a given prealignment, if any
1570   /// sort alitrackpoints w.r.t. global Y direction, if cosmics
1571   const Double_t kRad2L[6] = {5*5,10*10,18*18,30*30,40*40,60*60};
1572   const Float_t kSensSigY2[6] = {200e-4*200e-4/12, 200e-4*200e-4/12,
1573                                  300e-4*300e-4/12, 300e-4*300e-4/12,
1574                                  300e-4*300e-4/12, 300e-4*300e-4/12}; // thickness^2/12
1575   //
1576   fTrack = NULL;
1577   Int_t   idx[20] = {0};
1578   Short_t lrID[20] = {0};
1579   Int_t npts=atp->GetNPoints();
1580   if (npts<fMinNPtsPerTrack) return NULL;
1581   TGeoHMatrix hcov;
1582   //
1583   /// checks if AliTrackPoints belong to defined modules
1584   Int_t ngoodpts=0;
1585   Int_t intidx[20];
1586   for (int j=0; j<npts; j++) {
1587     intidx[j] = GetRequestedModID(atp->GetVolumeID()[j]);
1588     if (intidx[j]<0) continue;
1589     ngoodpts++;
1590     Float_t xx=atp->GetX()[j];
1591     Float_t yy=atp->GetY()[j];
1592     Float_t r=xx*xx + yy*yy;
1593     int lay;
1594     for (lay=0;lay<6;lay++) if (r<kRad2L[lay]) break;
1595     if (lay>5) continue;
1596     lrID[j] = lay;
1597   }
1598   //
1599   AliDebug(3,Form("Number of points in defined modules: %d out of %d",ngoodpts,npts));
1600
1601   // pepo270809
1602   Int_t nextra=0;
1603   // extra clusters selection mode
1604   if (fExtraClustersMode) {
1605     // 1 = keep one cluster, remove randomly the extra
1606     // 2 = keep one cluster, remove the internal one
1607     // 10 = keep tracks only if at least one extra is present
1608
1609     int iextra1[20],iextra2[20],layovl[20];
1610     // extra clusters mapping
1611     for (Int_t ipt=0; ipt<npts; ipt++) {
1612       if (intidx[ipt]<0) continue; // looks only defined modules...
1613       float p1x=atp->GetX()[ipt];
1614       float p1y=atp->GetY()[ipt];
1615       float p1z=atp->GetZ()[ipt];
1616       int lay1=int(AliGeomManager::VolUIDToLayer(atp->GetVolumeID()[ipt]));
1617       float r1 = p1x*p1x + p1y*p1y;
1618       UShort_t volid1=atp->GetVolumeID()[ipt];
1619
1620       for (int ik=ipt+1; ik<npts; ik++) {
1621         if (intidx[ik]<0) continue;
1622         // compare point ipt with next ones
1623         int lay2=int(AliGeomManager::VolUIDToLayer(atp->GetVolumeID()[ik]));
1624         // check if same layer
1625         if (lay2 != lay1) continue;
1626         UShort_t volid2=atp->GetVolumeID()[ik];
1627         // check if different module
1628         if (volid1 == volid2) continue;
1629
1630         float p2x=atp->GetX()[ik];
1631         float p2y=atp->GetY()[ik];
1632         float p2z=atp->GetZ()[ik];
1633         float r2 = p2x*p2x + p2y*p2y;
1634         float dr= (p1x-p2x)*(p1x-p2x) + (p1y-p2y)*(p1y-p2y) + (p1z-p2z)*(p1z-p2z);
1635
1636         // looks for pairs with dr<1 cm, same layer but different module
1637         if (dr<1.0) {
1638           // extra1 is the one with smaller radius in rphi plane
1639           if (r1<r2) {
1640             iextra1[nextra]=ipt;
1641             iextra2[nextra]=ik;
1642           }
1643           else {
1644             iextra1[nextra]=ik;
1645             iextra2[nextra]=ipt;
1646           }
1647           layovl[nextra]=lay1;
1648           nextra++;
1649         }
1650       }
1651     } // end overlaps mapping
1652
1653     // mode=1: keep only one clusters and remove the other randomly
1654     if (fExtraClustersMode==1 && nextra) {
1655       for (int ie=0; ie<nextra; ie++) {
1656         if (gRandom->Rndm()<0.5)
1657           intidx[iextra1[ie]]=-1;
1658         else
1659           intidx[iextra2[ie]]=-1;
1660       }
1661     }
1662
1663     // mode=2: keep only one clusters and remove the other...
1664     if (fExtraClustersMode==2 && nextra) {
1665       for (int ie=0; ie<nextra; ie++) {
1666         if (layovl[ie]==1) intidx[iextra2[ie]]=-1;
1667         else if (layovl[ie]==2) intidx[iextra1[ie]]=-1;
1668         else intidx[iextra1[ie]]=-1;
1669       }
1670     }
1671
1672     // mode=10: reject track if no overlaps are present
1673     if (fExtraClustersMode==10 && nextra==0) {
1674       AliInfo("Track with no extra clusters: rejected!");
1675       return NULL;
1676     }
1677
1678     // recalculate ngoodpts
1679     ngoodpts=0;
1680     for (int i=0; i<npts; i++) {
1681       if (intidx[i]>=0) ngoodpts++;
1682     }
1683   }
1684   // endpepo270809
1685
1686   // reject track if not enough points are left
1687   if (ngoodpts<fMinNPtsPerTrack) {
1688     AliDebug(2,"Track with not enough points!");
1689     return NULL;
1690   }
1691   // >> RS
1692   AliTrackPoint p;
1693   // check points in specific places
1694   if (fRequirePoints[fDataType]) {
1695     Int_t nlayup[6],nlaydown[6],nlay[6];
1696     Int_t ndetup[3],ndetdown[3],ndet[3];
1697     for (Int_t j=0; j<6; j++) {nlayup[j]=0; nlaydown[j]=0; nlay[j]=0;}
1698     for (Int_t j=0; j<3; j++) {ndetup[j]=0; ndetdown[j]=0; ndet[j]=0;}
1699
1700     for (int i=0; i<npts; i++) {
1701       // skip not defined points
1702       if (intidx[i]<0) continue;
1703       //
1704       Float_t yy=atp->GetY()[i];
1705       int lay = lrID[i];
1706       int det=lay/2;
1707       //printf("Point %d - x=%f  y=%f  R=%f  lay=%d  det=%d\n",i,xx,yy,r,lay,det);
1708
1709       if (yy>=0.0) { // UP point
1710         nlayup[lay]++;
1711         nlay[lay]++;
1712         ndetup[det]++;
1713         ndet[det]++;
1714       }
1715       else {
1716         nlaydown[lay]++;
1717         nlay[lay]++;
1718         ndetdown[det]++;
1719         ndet[det]++;
1720       }
1721     }
1722     //
1723     // checks minimum values
1724     Bool_t isok=kTRUE;
1725     for (Int_t j=0; j<6; j++) {
1726       if (nlayup[j]<fNReqLayUp[fDataType][j]) isok=kFALSE;
1727       if (nlaydown[j]<fNReqLayDown[fDataType][j]) isok=kFALSE;
1728       if (nlay[j]<fNReqLay[fDataType][j]) isok=kFALSE;
1729     }
1730     for (Int_t j=0; j<3; j++) {
1731       if (ndetup[j]<fNReqDetUp[fDataType][j]) isok=kFALSE;
1732       if (ndetdown[j]<fNReqDetDown[fDataType][j]) isok=kFALSE;
1733       if (ndet[j]<fNReqDet[fDataType][j]) isok=kFALSE;
1734     }
1735     if (!isok) {
1736       AliDebug(2,Form("Track does not meet all location point requirements!"));
1737       return NULL;
1738     }
1739   }
1740   // build a new track with (sorted) (prealigned) good points
1741   // pepo200709
1742   //fTrack = (AliTrackPointArray*)fTrackBuff[ngoodpts-fMinNPtsPerTrack];
1743   Int_t addVertex = IsTypeCollision()&&((fUseDiamond&&(fCheckDiamondPoint!=kDiamondIgnore))||(fUseVertex&&fVertexSet)) ? 1 : 0;
1744   fTrack = (AliTrackPointArray*)fTrackBuff[ngoodpts + addVertex ];
1745   if (!fTrack) {
1746     fTrack = new AliTrackPointArray(ngoodpts + addVertex);
1747     //    fTrackBuff.AddAtAndExpand(fTrack,ngoodpts-fMinNPtsPerTrack);
1748     fTrackBuff.AddAtAndExpand(fTrack,ngoodpts + addVertex);
1749   }
1750   //  fTrack = new AliTrackPointArray(ngoodpts);
1751   // endpepo200709
1752   //
1753   //
1754   for (int i=0; i<npts; i++) idx[i]=i;
1755   // sort track if required
1756   if (IsTypeCosmics()) TMath::Sort(npts,atp->GetY(),idx); // sort descending...
1757   //
1758   Int_t npto=0;
1759   if (fClusLoc.GetSize()<3*npts)    fClusLoc.Set(3*npts);
1760   if (fClusGlo.GetSize()<3*npts)    fClusGlo.Set(3*npts);
1761   if (fClusSigLoc.GetSize()<3*npts) fClusSigLoc.Set(3*npts);
1762   //
1763   for (int i=0; i<npts; i++) {
1764     // skip not defined points
1765     if (intidx[idx[i]]<0) continue;
1766     atp->GetPoint(p,idx[i]);
1767     int sid = AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());
1768     //
1769     // prealign point if required
1770     // get matrix used to produce the digits
1771     AliITSAlignMille2Module *mod = GetMilleModule(intidx[idx[i]]);
1772     TGeoHMatrix *svOrigMatrix = GetSensorOrigMatrixSID(sid); //mod->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
1773     // get back real local coordinate
1774     fMeasLoc  = fClusLoc.GetArray() + npto*3;
1775     fMeasGlo  = fClusGlo.GetArray() + npto*3;
1776     fSigmaLoc = fClusSigLoc.GetArray() + npto*3;
1777     fMeasGlo[0]=p.GetX();
1778     fMeasGlo[1]=p.GetY();
1779     fMeasGlo[2]=p.GetZ();
1780     AliDebug(3,Form("Global coordinates of measured point : X=%+f  Y=%+f  Z=%+f \n",fMeasGlo[0],fMeasGlo[1],fMeasGlo[2]));
1781     svOrigMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
1782     AliDebug(3,Form("Local coordinates of measured point : X=%+f  Y=%+f  Z=%+f \n",fMeasLoc[0],fMeasLoc[1],fMeasLoc[2]));
1783     //
1784     if (p.GetDriftTime()>0) ProcessSDDPointInfo(&p,sid, npto);     // for SDD points extract vdrift
1785     //
1786     // update covariance matrix
1787     Double_t hcovel[9];
1788     hcovel[0]=double(p.GetCov()[0]);
1789     hcovel[1]=double(p.GetCov()[1]);
1790     hcovel[2]=double(p.GetCov()[2]);
1791     hcovel[3]=double(p.GetCov()[1]);
1792     hcovel[4]=double(p.GetCov()[3]);
1793     hcovel[5]=double(p.GetCov()[4]);
1794     hcovel[6]=double(p.GetCov()[2]);
1795     hcovel[7]=double(p.GetCov()[4]);
1796     hcovel[8]=double(p.GetCov()[5]);
1797     hcov.SetRotation(hcovel);
1798     //
1799     if (AliLog::GetGlobalDebugLevel()>=2) {
1800       AliInfo("Original Global Cov Matrix");
1801       printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovel[0],hcovel[1],hcovel[2],hcovel[4],hcovel[5],hcovel[8]);
1802     }
1803     //
1804     // now rotate in local system
1805     hcov.Multiply(svOrigMatrix);
1806     hcov.MultiplyLeft(&svOrigMatrix->Inverse());
1807     // now hcov is LOCAL COVARIANCE MATRIX
1808     // apply sigma scaling
1809     Double_t *hcovscl = hcov.GetRotationMatrix();
1810     /*
1811     const float *cv = p.GetCov();
1812     printf("## %d %d  %+.3e %+.3e %+.3e   %+.3e %+.3e %+.3e   %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e   %+.3e %+.3e %+.3e\n",
1813            sid,p.GetClusterType(),
1814            fMeasGlo[0],fMeasGlo[1],fMeasGlo[2],
1815            fMeasLoc[0],fMeasLoc[1],fMeasLoc[2],
1816            cv[0],cv[1],cv[2],cv[3],cv[4],cv[5],
1817            hcovscl[0],hcovscl[4],hcovscl[8]);
1818
1819     */
1820     if (AliLog::GetGlobalDebugLevel()>=2) {
1821       AliInfo("Original Local Cov Matrix");
1822       printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1823     }
1824     hcovscl[4] = fUseLocalYErr ? kSensSigY2[lrID[idx[i]]] : 1E-8; // error due to the sensor thickness
1825     //
1826     for (int ir=3;ir--;) for (int ic=3;ic--;) {
1827         if (ir==ic) {
1828           if ( IsZero(hcovscl[ir*3+ic],1e-8) ) hcovscl[ir*3+ic] = 1E-8;
1829           else hcovscl[ir*3+ic] *= mod->GetSigmaFactor(ir)*mod->GetSigmaFactor(ic); //RRR
1830           fSigmaLoc[ir] = TMath::Sqrt(hcovscl[ir*3+ic]);
1831         }
1832         else hcovscl[ir*3+ic]  = 0;
1833       }
1834     //
1835     if (AliLog::GetGlobalDebugLevel()>=2) {
1836       AliInfo("Modified Local Cov Matrix");
1837       printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1838     }
1839     //
1840     if (fBug==1) {
1841       // correzione bug LAYER 5  SSD temporanea..
1842       int ssdidx=AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());
1843       if (ssdidx>=500 && ssdidx<1248) {
1844         int ladder=(ssdidx-500)%22;
1845         if (ladder==18) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx+1));
1846         if (ladder==19) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx-1));
1847       }
1848     }
1849     /// get (evenctually prealigned) matrix of sens. vol.
1850     TGeoHMatrix *svMatrix = GetSensorCurrMatrixSID(sid);    //mod->GetSensitiveVolumeMatrix(p.GetVolumeID());
1851     // modify global coordinates according with pre-aligment
1852     svMatrix->LocalToMaster(fMeasLoc,fMeasGlo);
1853     // now rotate in local system
1854     hcov.Multiply(&svMatrix->Inverse());
1855     hcov.MultiplyLeft(svMatrix);         // hcov is back in GLOBAL RF
1856     // cure once more
1857     for (int ir=3;ir--;) for (int ic=3;ic--;) if (IsZero(hcovscl[ir*3+ic])) hcovscl[ir*3+ic] = 0.;
1858     //    printf("\nErrMatGlob: after\n"); hcov.Print(""); //RRR
1859     //
1860     if (AliLog::GetGlobalDebugLevel()>=2) {
1861       AliInfo("Modified Global Cov Matrix");
1862       printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1863     }
1864     //
1865     Float_t pcov[6];
1866     pcov[0]=hcovscl[0];
1867     pcov[1]=hcovscl[1];
1868     pcov[2]=hcovscl[2];
1869     pcov[3]=hcovscl[4];
1870     pcov[4]=hcovscl[5];
1871     pcov[5]=hcovscl[8];
1872     //
1873     // make sure the matrix is positive definite
1874     {
1875       enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5};
1876       if (pcov[kXX]*pcov[kYY]*0.999<pcov[kXY]*pcov[kXY]) pcov[kXY] = 0.999*TMath::Sign((float)TMath::Sqrt(pcov[kXX]*pcov[kYY]),pcov[kXY]);
1877       if (pcov[kXX]*pcov[kZZ]*0.999<pcov[kXZ]*pcov[kXZ]) pcov[kXZ] = 0.999*TMath::Sign((float)TMath::Sqrt(pcov[kXX]*pcov[kZZ]),pcov[kXZ]);
1878       if (pcov[kYY]*pcov[kZZ]*0.999<pcov[kYZ]*pcov[kYZ]) pcov[kYZ] = 0.999*TMath::Sign((float)TMath::Sqrt(pcov[kYY]*pcov[kZZ]),pcov[kYZ]);
1879     }
1880     //
1881     p.SetXYZ(fMeasGlo[0],fMeasGlo[1],fMeasGlo[2],pcov);
1882     //    printf("New Gl coordinates of measured point : X=%f  Y=%f  Z=%f \n",fMeasGlo[0],fMeasGlo[1],fMeasGlo[2]);
1883     AliDebug(3,Form("New global coordinates of measured point : X=%+f  Y=%+f  Z=%+f \n",fMeasGlo[0],fMeasGlo[1],fMeasGlo[2]));
1884     fTrack->AddPoint(npto,&p);
1885     AliDebug(2,Form("Adding point[%d] = ( %+f , %+f , %+f )     volid = %d",npto,fTrack->GetX()[npto],
1886                     fTrack->GetY()[npto],fTrack->GetZ()[npto],fTrack->GetVolumeID()[npto] ));
1887     //    printf("Adding %d %d %f\n",npto, p.GetVolumeID(), p.GetY());
1888     npto++;
1889   }
1890   //
1891   fDiamondPointID = -1;
1892   if (addVertex) {
1893     fTrack->AddPoint(npto,&fDiamond);
1894     fMeasLoc  = fClusLoc.GetArray() + npto*3;
1895     fMeasGlo  = fClusGlo.GetArray() + npto*3;
1896     fSigmaLoc = fClusSigLoc.GetArray() + npto*3;
1897     fMeasLoc[0] = fMeasGlo[0] = fDiamond.GetX();
1898     fMeasLoc[1] = fMeasGlo[1] = fDiamond.GetY();
1899     fMeasLoc[2] = fMeasGlo[2] = fDiamond.GetZ();
1900     fSigmaLoc[0] = TMath::Sqrt(fDiamond.GetCov()[0]);
1901     fSigmaLoc[1] = TMath::Sqrt(fDiamond.GetCov()[3]);
1902     fSigmaLoc[2] = TMath::Sqrt(fDiamond.GetCov()[5]);
1903     fDiamondPointID = npto++;
1904   }
1905   //
1906   return fTrack;
1907 }
1908
1909 //________________________________________________________________________________________________________
1910 AliTrackPointArray *AliITSAlignMille2::SortTrack(const AliTrackPointArray *atp)
1911 {
1912   /// sort alitrackpoints w.r.t. global Y direction
1913   AliTrackPointArray *atps=NULL;
1914   Int_t idx[20];
1915   Int_t npts=atp->GetNPoints();
1916   AliTrackPoint p;
1917   atps=new AliTrackPointArray(npts);
1918
1919   TMath::Sort(npts,atp->GetY(),idx);
1920
1921   for (int i=0; i<npts; i++) {
1922     atp->GetPoint(p,idx[i]);
1923     atps->AddPoint(i,&p);
1924     AliDebug(2,Form("Point[%d] = ( %+f , %+f , %+f )     volid = %d",i,atps->GetX()[i],atps->GetY()[i],atps->GetZ()[i],atps->GetVolumeID()[i] ));
1925   }
1926   return atps;
1927 }
1928
1929 //________________________________________________________________________________________________________
1930 Int_t AliITSAlignMille2::GetCurrentLayer() const
1931 {
1932   // get current layer id
1933   if (!fGeoManager) {
1934     AliInfo("ITS geometry not initialized!");
1935     return -1;
1936   }
1937   return (Int_t)AliGeomManager::VolUIDToLayer(fCluster.GetVolumeID());
1938 }
1939
1940 //________________________________________________________________________________________________________
1941 Int_t AliITSAlignMille2::InitModuleParams()
1942 {
1943   /// initialize geometry parameters for a given detector
1944   /// for current cluster (fCluster)
1945   /// fGlobalInitParam[] is set as:
1946   ///    [tx,ty,tz,psi,theta,phi]
1947   ///    (old was [tx,ty,tz,theta,psi,phi] ROOT's angles...)
1948   /// *** At the moment: using Raffalele's angles definition ***
1949   ///
1950   /// return 0 if success
1951   /// If module is found but has no parameters to vary, return 1
1952
1953   if (!fGeoManager) {
1954     AliInfo("ITS geometry not initialized!");
1955     return -1;
1956   }
1957
1958   // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster)
1959
1960   // set the internal index (index in module list)
1961   UShort_t voluid=fCluster.GetVolumeID();
1962   fCurrentSensID = AliITSAlignMille2Module::GetIndexFromVolumeID(voluid);
1963   //
1964   if (fCurrentSensID==-1) { // this is a special "vertex" module
1965     fCurrentModule = GetMilleModuleByVID(voluid);
1966     fCurrentSensID = fCurrentModule->GetIndex();
1967
1968   }
1969   else {
1970     //
1971     // IT IS VERY IMPORTANT: start from the end of the list, where the childs are located !!!
1972     Int_t k=fNModules-1;
1973     fCurrentModule = 0;
1974     // VERY IMPORTANT: if the sensors were explicitly provided, don't look in the supermodules
1975     while (k>=0 && ! (fCurrentModule=GetMilleModule(k))->IsIn(voluid)) k--;
1976     if (k<0) return -3;
1977   }
1978   //
1979   for (int i=AliITSAlignMille2Module::kMaxParTot;i--;) fModuleInitParam[i] = 0.0;
1980   //
1981   int clID = fCluster.GetUniqueID()-1;
1982   if (clID<0) { // external cluster
1983     fMeasGlo  = &fExtClusterPar[0];
1984     fMeasLoc  = &fExtClusterPar[3];
1985     fSigmaLoc = &fExtClusterPar[6];
1986     fExtClusterPar[0] = fCluster.GetX();
1987     fExtClusterPar[1] = fCluster.GetY();
1988     fExtClusterPar[2] = fCluster.GetZ();
1989     //
1990     TGeoHMatrix *svMatrix = fCurrentModule->GetSensitiveVolumeMatrix(voluid);
1991     svMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
1992     TGeoHMatrix hcov;
1993     Double_t hcovel[9];
1994     hcovel[0]=double(fCluster.GetCov()[0]);
1995     hcovel[1]=double(fCluster.GetCov()[1]);
1996     hcovel[2]=double(fCluster.GetCov()[2]);
1997     hcovel[3]=double(fCluster.GetCov()[1]);
1998     hcovel[4]=double(fCluster.GetCov()[3]);
1999     hcovel[5]=double(fCluster.GetCov()[4]);
2000     hcovel[6]=double(fCluster.GetCov()[2]);
2001     hcovel[7]=double(fCluster.GetCov()[4]);
2002     hcovel[8]=double(fCluster.GetCov()[5]);
2003     hcov.SetRotation(hcovel);
2004     // now rotate in local system
2005     hcov.Multiply(svMatrix);
2006     hcov.MultiplyLeft(&svMatrix->Inverse());
2007     if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
2008     if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
2009     //
2010   }
2011   else {
2012     int offs = 3*clID;
2013     fMeasGlo  = fClusGlo.GetArray()  + offs;
2014     fMeasLoc  = fClusLoc.GetArray()  + offs;
2015     fSigmaLoc = fClusSigLoc.GetArray() + offs;
2016   }
2017   //
2018   // set minimum value for SigmaLoc to 10 micron
2019   if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
2020   if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
2021   if (fCurrentSensID==kVtxSensID || fUseLocalYErr) if (fSigmaLoc[1]<0.0010) fSigmaLoc[1]=0.0010;
2022   //
2023   AliDebug(2,Form("Local coordinates of measured point : X=%+f  Y=%+f  Z=%+f \n",fMeasLoc[0] ,fMeasLoc[1] ,fMeasLoc[2] ));
2024   AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%g  fSigmaLocY=%g fSigmaLocZ=%g \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
2025   //
2026   return 0;
2027 }
2028
2029 //________________________________________________________________________________________________________
2030 void AliITSAlignMille2::Print(Option_t*) const
2031 {
2032   // print current status
2033   printf("*** AliMillepede for ITS ***\n");
2034   printf("    Number of defined super modules: %d\n",fNModules);
2035   printf("    Obtained parameters refer to %s Deltas\n",fUseGlobalDelta ? "GLOBAL":"LOCAL");
2036   //
2037   if (fGeoManager)
2038     printf("    geometry loaded from %s\n",fGeometryPath.Data());
2039   else
2040     printf("    geometry not loaded\n");
2041   //
2042   if (fUsePreAlignment)
2043     printf("    using prealignment from %s \n",fPreDeltaPath.Data());
2044   else
2045     printf("    prealignment not used\n");
2046   //
2047   //
2048   if (fBOn)
2049     printf("    B Field set to %+f T - using helices\n",fBField);
2050   else
2051     printf("    B Field OFF - using straight lines \n");
2052   //
2053   if (fTPAFitter)
2054     printf("    Using AliITSTPArrayFit class for track fitting\n");
2055   else
2056     printf("    Using StraightLine/Riemann fitter for track fitting\n");
2057   //
2058   printf("Using local Y error due to the sensor thickness: %s\n",(fUseLocalYErr && fTPAFitter) ? "ON":"OFF");
2059   //
2060   for (int itp=0;itp<kNDataType;itp++) {
2061     if (fRequirePoints[itp]) printf("    Required points in %s tracks:\n",itp==kCosmics? "cosmics" : "collisions");
2062     for (Int_t i=0; i<6; i++) {
2063       if (fNReqLayUp[itp][i]>0) printf("        Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[itp][i]);
2064       if (fNReqLayDown[itp][i]>0) printf("        Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[itp][i]);
2065       if (fNReqLay[itp][i]>0) printf("        Layer %d : %d points \n",i+1,fNReqLay[itp][i]);
2066     }
2067     for (Int_t i=0; i<3; i++) {
2068       if (fNReqDetUp[itp][i]>0) printf("        Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[itp][i]);
2069       if (fNReqDetDown[itp][i]>0) printf("        Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[itp][i]);
2070       if (fNReqDet[itp][i]>0) printf("        Detector %d : %d points \n",i+1,fNReqDet[itp][i]);
2071     }
2072   }
2073   printf("        SDD VDrift correction         : %s",fIsSDDVDriftMult ? "Mult":"Add");
2074   printf("        Weight acc. to pT in power    : %f",fWeightPt);
2075   //
2076   printf("\n    Millepede configuration parameters:\n");
2077   printf("        init factor for chi2 cut      : %.4f\n",fStartFac);
2078   printf("        final factor for chi2 cut     : %.4f\n",fFinalFac);
2079   printf("        first iteration cut value     : %.4f\n",fResCutInitial);
2080   printf("        other iterations cut value    : %.4f\n",fResCut);
2081   printf("        number of stddev for chi2 cut : %d\n",fNStdDev);
2082   printf("        def.scaling for local sigmas  : %.4f %.4f %.4f\n",fSigmaFactor[0],fSigmaFactor[1],fSigmaFactor[2]);
2083   printf("        min.tracks per module         : %d\n",fMinPntPerSens);
2084   //
2085   printf("List of defined modules:\n");
2086   printf("  intidx\tindex\tvoluid\tname\n");
2087   for (int i=0; i<fNModules; i++) {
2088     AliITSAlignMille2Module* md = GetMilleModule(i);
2089     printf("  %d\t%d\t%d\t%s\n",i,md->GetIndex(),md->GetVolumeID(),md->GetName());
2090   }
2091 }
2092
2093 //________________________________________________________________________________________________________
2094 AliITSAlignMille2Module  *AliITSAlignMille2::GetMilleModuleByVID(UShort_t voluid) const
2095 {
2096   // return pointer to a defined supermodule
2097   // return NULL if error
2098   Int_t i=IsVIDDefined(voluid);
2099   if (i<0) return NULL;
2100   return GetMilleModule(i);
2101 }
2102
2103 //________________________________________________________________________________________________________
2104 AliITSAlignMille2Module  *AliITSAlignMille2::GetMilleModuleBySymName(const Char_t* symname) const
2105 {
2106   // return pointer to a defined supermodule
2107   // return NULL if error
2108   UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
2109   if (vid>0) return GetMilleModuleByVID(vid);
2110   else {    // this is not alignable module, need to look within defined supermodules
2111     int i = IsSymDefined(symname);
2112     if (i>=0) return  GetMilleModule(i);
2113   }
2114   return 0;
2115 }
2116
2117 //________________________________________________________________________________________________________
2118 AliITSAlignMille2Module  *AliITSAlignMille2::GetMilleModuleIfContained(const Char_t* symname) const
2119 {
2120   // return pointer to a defined/contained supermodule
2121   // return NULL otherwise
2122   int i = IsSymContained(symname);
2123   return i<0 ? 0 : GetMilleModule(i);
2124 }
2125
2126 //________________________________________________________________________________________________________
2127 AliAlignObjParams* AliITSAlignMille2::GetPrealignedObject(const Char_t* symname) const
2128 {
2129   // get delta from prealignment for given volume
2130   if (!fPrealignment) return 0;
2131   for (int ipre=fPrealignment->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
2132     AliAlignObjParams* preob = (AliAlignObjParams*)fPrealignment->At(ipre);
2133     if (!strcmp(preob->GetSymName(),symname)) return preob;
2134   }
2135   return 0;
2136 }
2137
2138 //________________________________________________________________________________________________________
2139 AliAlignObjParams* AliITSAlignMille2::GetConstrRefObject(const Char_t* symname) const
2140 {
2141   // get delta with respect to which the constraint is declared
2142   if (!fConstrRef) return 0;
2143   for (int ipre=fConstrRef->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
2144     AliAlignObjParams* preob = (AliAlignObjParams*)fConstrRef->At(ipre);
2145     if (!strcmp(preob->GetSymName(),symname)) return preob;
2146   }
2147   return 0;
2148 }
2149
2150 //________________________________________________________________________________________________________
2151 Bool_t AliITSAlignMille2::InitRiemanFit()
2152 {
2153   // Initialize Riemann Fitter for current track
2154   // return kFALSE if error
2155
2156   if (!fBOn) return kFALSE;
2157
2158   Int_t npts=0;
2159   AliTrackPoint ap;
2160   npts = fTrack->GetNPoints();
2161   AliDebug(3,Form("Fitting track with %d points",npts));
2162   if (!fRieman) fRieman = new AliTrackFitterRieman();
2163   fRieman->Reset();
2164   fRieman->SetTrackPointArray(fTrack);
2165
2166   TArrayI ai(npts);
2167   for (Int_t ipt=0; ipt<npts; ipt++) ai[ipt]=fTrack->GetVolumeID()[ipt];
2168
2169   // fit track with 5 params in his own tracking-rotated reference system
2170   // xc = -p[1]/p[0];
2171   // yc = 1/p[0];
2172   // R  = sqrt( x0*x0 + y0*y0 - y0*p[2]);
2173   if (!fRieman->Fit(&ai,NULL,(AliGeomManager::ELayerID)1,(AliGeomManager::ELayerID)6)) {
2174     return kFALSE;
2175   }
2176
2177   for (int i=0; i<5; i++)
2178     fLocalInitParam[i] = fRieman->GetParam()[i];
2179
2180   return kTRUE;
2181 }
2182
2183 //________________________________________________________________________________________________________
2184 void trackFit2D(Int_t &, Double_t *, double &chi2, double *par, int flag)
2185 {
2186   // local function for minuit
2187   const double kTiny = 1.e-14;
2188   chi2 = 0;
2189   static AliTrackPoint pnt;
2190   static Bool_t fullErr2D;
2191   //
2192   if (flag==1) fullErr2D = kFALSE;//kTRUE;
2193   //  fullErr2D = kTRUE;
2194   enum {kAX,kAZ,kBX,kBZ};
2195   enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5};
2196   //
2197   AliITSAlignMille2* alig = AliITSAlignMille2::GetInstance();
2198   AliTrackPointArray* track = alig->GetCurrentTrack();
2199   //
2200   int npts = track->GetNPoints();
2201   for (int ip=0;ip<npts;ip++) {
2202     track->GetPoint(pnt,ip);
2203     const float *cov = pnt.GetCov();
2204     double y  = pnt.GetY();
2205     double dx = pnt.GetX() - (par[kAX]+y*par[kBX]);
2206     double dz = pnt.GetZ() - (par[kAZ]+y*par[kBZ]);
2207     double xxe = cov[kXX];
2208     double zze = cov[kZZ];
2209     double xze = cov[kXZ];
2210     //
2211     if (fullErr2D) {
2212       xxe += par[kBX]*par[kBX]*cov[kYY]-2.*par[kBX]*cov[kXY];
2213       zze += par[kBZ]*par[kBZ]*cov[kYY]-2.*par[kBZ]*cov[kZY];
2214       xze += par[kBX]*par[kBZ]*cov[kYY]-cov[kYZ]*par[kBZ]-cov[kXY]*par[kBX];
2215     }
2216     //
2217     double det = xxe*zze - xze*xze;
2218     if (det<kTiny) {
2219       printf("Negative diag. error (det=%+e) |sxx:%+e szz:%+e sxz:%+e| bx:%+e bz:%+e|\n"
2220              "Discarding correlation term\n",det,xxe,zze,xze,par[kBX],par[kBZ]);
2221       xxe = cov[kXX];
2222       zze = cov[kZZ];
2223       xze = cov[kXZ];
2224       fullErr2D = kFALSE;
2225     }
2226     double xxeI = zze/det;
2227     double zzeI = xxe/det;
2228     double xzeI =-xze/det;
2229     //
2230     chi2 += dx*dx*xxeI + dz*dz*zzeI + 2.*dx*dz*xzeI;
2231     //
2232     //    printf("%d | %+e %+e %+e %+e %+e -> %+e\n",ip,dx,dz,xxeI,zzeI,xzeI,  chi2);
2233   }
2234   //
2235 }
2236
2237 //________________________________________________________________________________________________________
2238 void AliITSAlignMille2::InitTrackParams(int meth)
2239 {
2240   /// initialize local parameters with different methods
2241   /// for current track (fTrack)
2242   Int_t npts=0;
2243   AliTrackPoint ap;
2244   double sX=0,sXY=0,sZ=0,sZY=0,sY=0,sYY=0,det=0;
2245   // simple linear interpolation
2246   // get local starting parameters (to be substituted by ESD track parms)
2247   // local parms (fLocalInitParam[]) are:
2248   //      [0] = global x coord. of straight line intersection at y=0 plane
2249   //      [1] = global z coord. of straight line intersection at y=0 plane
2250   //      [2] = px/py
2251   //      [3] = pz/py
2252   // test #1: linear fit in x(y) and z(y)
2253   npts = fTrack->GetNPoints();
2254   AliDebug(3,Form("*** initializing track with %d points ***",npts));
2255   for (int i=npts;i--;) {
2256     sY  += fTrack->GetY()[i];
2257     sYY += fTrack->GetY()[i]*fTrack->GetY()[i];
2258     sX  += fTrack->GetX()[i];
2259     sXY += fTrack->GetX()[i]*fTrack->GetY()[i];
2260     sZ  += fTrack->GetZ()[i];
2261     sZY += fTrack->GetZ()[i]*fTrack->GetY()[i];
2262   }
2263   det = sYY*npts-sY*sY;
2264   if (IsZero(det)) det = 1E-16;
2265   fLocalInitParam[0] = (sX*sYY-sY*sXY)/det;
2266   fLocalInitParam[2] = (sXY*npts-sY*sX)/det;
2267   //
2268   fLocalInitParam[1] = (sZ*sYY-sY*sZY)/det;
2269   fLocalInitParam[3] = (sZY*npts-sY*sZ)/det;
2270   // pepo200709
2271   fLocalInitParam[4] = 0.0;
2272   // endpepo200709
2273
2274   AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %+f    ugx = %+f\n",fLocalInitParam[0],fLocalInitParam[2]));
2275   //
2276   if (meth==1) return;
2277   //
2278   // perform full fit accounting for cov.matrix
2279   static TVirtualFitter *minuit = 0;
2280   static Double_t step[5]   = {1E-3,1E-3,1E-4,1E-4,1E-5};
2281   static Double_t arglist[10];
2282   //
2283   if (!minuit) {
2284     minuit = TVirtualFitter::Fitter(0,4);
2285     minuit->SetFCN(trackFit2D);
2286     arglist[0] = 1;
2287     minuit->ExecuteCommand("SET ERR",arglist, 1);
2288     //
2289     arglist[0] = -1;
2290     minuit->ExecuteCommand("SET PRINT",arglist,1);
2291     //
2292   }
2293   //
2294   minuit->SetParameter(0, "ax",   fLocalInitParam[0], step[0], 0,0);
2295   minuit->SetParameter(1, "az",   fLocalInitParam[1], step[1], 0,0);
2296   minuit->SetParameter(2, "bx",   fLocalInitParam[2], step[2], 0,0);
2297   minuit->SetParameter(3, "bz",   fLocalInitParam[3], step[3], 0,0);
2298   //
2299   arglist[0] = 1000; // number of function calls
2300   arglist[1] = 0.001; // tolerance
2301   minuit->ExecuteCommand("MIGRAD",arglist,2);
2302   //
2303   for (int i=0;i<4;i++) fLocalInitParam[i] = minuit->GetParameter(i);
2304   for (int i=0;i<4;i++) for (int j=0;j<4;j++) fLocalInitParEr[i][j] = minuit->GetCovarianceMatrixElement(i,j);
2305   /*
2306   double amin,edm,errdef;
2307   int nvpar,nparx;
2308   minuit->GetStats(amin,edm,errdef,nvpar,nparx);
2309   amin /= (2*npts - 4);
2310   printf("Mchi2: %+e\n",amin);
2311   */
2312   //
2313 }
2314
2315 //________________________________________________________________________________________________________
2316 Int_t AliITSAlignMille2::IsSymDefined(const Char_t* symname) const
2317 {
2318   // checks if supermodule with this symname is defined and return the internal index
2319   // return -1 if not.
2320   for (int k=fNModules;k--;) if (!strcmp(symname,GetMilleModule(k)->GetName())) return k;
2321   return -1;
2322 }
2323
2324 //________________________________________________________________________________________________________
2325 Int_t AliITSAlignMille2::IsSymContained(const Char_t* symname) const
2326 {
2327   // checks if module with this symname is defined and return the internal index
2328   // return -1 if not.
2329   UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
2330   if (vid>0) return IsVIDContained(vid);
2331   // only sensors have real vid, but maybe we have a supermodule with fake vid?
2332   // IMPORTANT: always start from the end to start from the sensors
2333   return IsSymDefined(symname);
2334 }
2335
2336 //________________________________________________________________________________________________________
2337 Int_t AliITSAlignMille2::IsVIDDefined(UShort_t voluid) const
2338 {
2339   // checks if supermodule 'voluid' is defined and return the internal index
2340   // return -1 if not.
2341   for (int k=fNModules;k--;) if (voluid==GetMilleModule(k)->GetVolumeID()) return k;
2342   return -1;
2343 }
2344
2345 //________________________________________________________________________________________________________
2346 Int_t AliITSAlignMille2::IsVIDContained(UShort_t voluid) const
2347 {
2348   // checks if the sensitive module 'voluid' is contained inside a supermodule
2349   // and return the internal index of the last identified supermodule
2350   // return -1 if error
2351   // IMPORTANT: always start from the end to start from the sensors
2352   if (AliITSAlignMille2Module::GetIndexFromVolumeID(voluid)<0) return -1;
2353   for (int k=fNModules;k--;) if (GetMilleModule(k)->IsIn(voluid)) return k;
2354   return -1;
2355 }
2356
2357 //________________________________________________________________________________________________________
2358 Int_t AliITSAlignMille2::GetRequestedModID(UShort_t voluid) const
2359 {
2360   // checks if the sensitive module 'voluid' is contained inside a supermodule
2361   // and return the internal index of the last identified supermodule
2362   // return -1 if error
2363   // IMPORTANT: always start from the end to start from the sensors
2364   if (AliITSAlignMille2Module::GetIndexFromVolumeID(voluid)<0) return -1;
2365   int k;
2366   for (k=fNModules;k--;) if (GetMilleModule(k)->IsIn(voluid)) break;
2367   if (k<0) return -1;
2368   AliITSAlignMille2Module* md = GetMilleModule(k);
2369   while (md && md->IsNotInConf()) md = md->GetParent();
2370   return md ? md->GetUniqueID() : -1;
2371 }
2372
2373 //________________________________________________________________________________________________________
2374 Int_t AliITSAlignMille2::CheckCurrentTrack()
2375 {
2376   /// checks if AliTrackPoints belongs to defined modules
2377   /// return number of good poins
2378   /// return 0 if not enough points
2379
2380   Int_t npts = fTrack->GetNPoints();
2381   Int_t ngoodpts=0;
2382   // debug points
2383   for (int j=0; j<npts; j++) if (IsVIDContained(fTrack->GetVolumeID()[j])>=0) ngoodpts++;
2384   //
2385   if (ngoodpts<fMinNPtsPerTrack) return 0;
2386
2387   return ngoodpts;
2388 }
2389
2390 //________________________________________________________________________________________________________
2391 Int_t AliITSAlignMille2::ProcessTrack(const AliTrackPointArray *track, Double_t wgh)
2392 {
2393   /// Process track; Loop over hits and set local equations
2394   /// here 'track' is a AliTrackPointArray
2395   /// return 0 if success;
2396   //
2397   if (!fIsMilleInit) Init();
2398   //
2399   Int_t npts = track->GetNPoints();
2400   AliDebug(2,Form("*** Input track with %d points ***",npts));
2401
2402   // preprocessing of the input track: keep only points in defined volumes,
2403   // move points if prealignment is set, sort by Yglo if required
2404   fTrackWeight = wgh;
2405   fTrack=PrepareTrack(track);
2406   if (!fTrack) {
2407     RemoveHelixFitConstraint();
2408     RemoveVertexConstraint();
2409     return -1;
2410   }
2411   npts = fTrack->GetNPoints();
2412   if (npts>kMaxPoints) {
2413     AliError(Form("Compiled with kMaxPoints=%d, current track has %d points",kMaxPoints,npts));
2414   }
2415   AliDebug(2,Form("*** Processing prepared track with %d points ***",npts));
2416   //
2417   npts = FitTrack();
2418   if (npts<0) return npts;
2419   //
2420   //  printf("Params: "); for (int i=0;i<fNLocal;i++) printf("%+.2e ",fLocalInitParam[i]); printf("\n");//RRR
2421   Int_t nloceq=0;
2422   Int_t ngloeq=0;
2423   static Mille2Data md[kMaxPoints];
2424   //
2425   for (Int_t ipt=0; ipt<npts; ipt++) {
2426     fTrack->GetPoint(fCluster,ipt);
2427     fCluster.SetUniqueID(ipt+1);
2428     AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
2429
2430     // set geometry parameters for the the current module
2431     if (InitModuleParams()) continue;
2432     AliDebug(2,Form("    VolID=%d  Index=%d  InternalIdx=%d  symname=%s\n",
2433                     track->GetVolumeID()[ipt], fCurrentModule->GetIndex(),
2434                     fCurrentModule->GetUniqueID(), AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
2435     AliDebug(2,Form("    Preprocessed Point = ( %+f , %+f , %+f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
2436     int res = fTPAFitter ? AddLocalEquationTPA(md[nloceq]) : AddLocalEquation(md[nloceq]);
2437     if (res<0) {fTotBadLocEqPoints++; nloceq = 0; break;}
2438     else if (res==0) nloceq++;
2439     else {nloceq++; ngloeq++;}
2440   } // end loop over points
2441   //
2442   fTrack=NULL;
2443   // not enough good points?
2444   if (nloceq<fMinNPtsPerTrack || ngloeq<1) return -1;
2445   //
2446   // finally send local equations to millepede
2447   SetLocalEquations(md,nloceq);
2448   fMillepede->SaveRecordData(); // RRR
2449   //
2450   return 0;
2451 }
2452
2453 //________________________________________________________________________________________________________
2454 Int_t AliITSAlignMille2::FitTrack()
2455 {
2456   // Fit the track with selected constraints
2457   //
2458   const Double_t kfDiamondTolerance = 0.1;  //diamond tolerance on top of the MS error
2459   if (!fTrack) return -1;
2460   int npts = fTrack->GetNPoints();
2461   //
2462   if (fTPAFitter) {  // use dediacted fitter
2463     //
2464     // if the diamond point is attached, for the moment don't include it in the fit
2465     fTPAFitter->AttachPoints(fTrack,0, npts-1);
2466     fTPAFitter->SetBz(fBField);
2467     fTPAFitter->SetTypeCosmics(IsTypeCosmics());
2468     if (fIniTrackParamsMeth==1) fTPAFitter->SetIgnoreCov();
2469     //
2470     double chi2;
2471     double chi2f = 0;
2472     double dca2err;
2473     double dca2 = 0.;
2474     Bool_t fitIsDone = kFALSE;
2475     if (fUseDiamond && fDiamondPointID>0 && fCheckDiamondPoint==kDiamondCheckIfPrompt) { // diamond constraint was added, check if the track looks like prompt
2476       fTPAFitter->SetFirstLast(0,fDiamondPointID-1);
2477       if (IsCovIScaleTouched()) for (int i=npts;i--;) fTPAFitter->SetCovIScale(i,GetCovIScale(i));
2478       //
2479       chi2f = fTPAFitter->Fit(fConstrCharge,fConstrPT,fConstrPTErr);
2480       if ( chi2f<0 || (chi2f>fNStdDev*fStartFac && fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations()) ) { //RRR
2481         AliInfo(Form("Track fit failed on checking if it is prompt! skipping this track... Chi2:%+e",chi2f));
2482         fTPAFitter->Reset();
2483         //      fTrack = NULL;
2484         return -5;
2485       }
2486       double xyzRes[3];
2487       fTPAFitter->GetResiduals(xyzRes,&fDiamondI,kTRUE);
2488       dca2 = xyzRes[0]*xyzRes[0] + xyzRes[1]*xyzRes[1];
2489       double pT = IsFieldON() ? fTPAFitter->GetPt() : 0.45;
2490       if (pT<0.1) pT = 0.1;
2491       dca2err = kfDiamondTolerance + 0.02/pT;
2492       if (dca2>dca2err*dca2err) { // this is secondary
2493         int* clst = (int*) fTrack->GetClusterType();
2494         clst[fDiamondPointID] = -1;;
2495         fDiamondPointID = -1;
2496         fitIsDone = kTRUE;
2497         npts--;
2498       }
2499       else fTPAFitter->SetFirstLast(0,fDiamondPointID); // fit with diamond
2500     }
2501     //    fTPAFitter->SetParAxis(1);
2502     if (!fitIsDone) {
2503       if (IsCovIScaleTouched()) for (int i=npts;i--;) fTPAFitter->SetCovIScale(i,GetCovIScale(i));
2504       chi2 = fTPAFitter->Fit(fConstrCharge,fConstrPT,fConstrPTErr);
2505     }
2506     //
2507     RemoveHelixFitConstraint();  // suppress eventual constraints to not affect fit of the next track
2508     RemoveVertexConstraint(); // same ...
2509     //
2510     if ( !fitIsDone && (chi2<0 || (chi2>fNStdDev*fStartFac && fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations())) ) { //RRR
2511       AliInfo(Form("Track fit failed! skipping this track... Chi2:%+e",chi2));
2512       if (fUseDiamond && fDiamondPointID>0 && fCheckDiamondPoint==kDiamondCheckIfPrompt) AliInfo(Form("VertexFree fit gave Chi2:%+e with residual %+e",chi2f,TMath::Sqrt(dca2)));
2513       /*
2514         fTrack->Print("");
2515         fTPAFitter->FitHelixCrude();
2516         fTPAFitter->SetFitDone();
2517         fTPAFitter->Print();
2518       */
2519       fTPAFitter->Reset();
2520       //      fTrack = NULL;
2521       return -5;
2522     }
2523     fNLocal = fTPAFitter->IsFieldON() ? 5:4; // Attention: the fitter might have decided to work in line mode
2524     npts  = fTPAFitter->GetLast() - fTPAFitter->GetFirst() + 1; // actual number of points
2525     /*
2526       double *pr = fTPAFitter->GetParams();
2527       printf("FtPar: %+.5e  %+.5e  %+.5e  %+.5e | chi2:%.3e\n",pr[2],pr[0],pr[3],pr[1],chi2); // RRR
2528     */
2529   }
2530   else {
2531     //
2532     if (!fBOn) { // straight lines
2533       // set local starting parameters (to be substituted by ESD track parms)
2534       // local parms (fLocalInitParam[]) are:
2535       //      [0] = global x coord. of straight line intersection at y=0 plane
2536       //      [1] = global z coord. of straight line intersection at y=0 plane
2537       //      [2] = px/py
2538       //      [3] = pz/py
2539       InitTrackParams(fIniTrackParamsMeth);
2540       /*
2541       double *pr = fLocalInitParam;
2542       printf("FtPar: %+.5e  %+.5e  %+.5e  %+.5e |\n",pr[0],pr[1],pr[2],pr[3]); // RRR
2543       */
2544     }
2545     else {
2546       // local parms (fLocalInitParam[]) are the Riemann Fitter params
2547       if (!InitRiemanFit()) {
2548         AliInfo("Riemann fit failed! skipping this track...");
2549         fTrack=NULL;
2550         return -5;
2551       }
2552     }
2553   }
2554   return npts;
2555   //
2556 }
2557
2558 //________________________________________________________________________________________________________
2559 Int_t AliITSAlignMille2::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar)
2560 {
2561   /// calculate track intersection point in local coordinates
2562   /// according with a given set of parameters (local(4) and global(6))
2563   /// and fill fPintLoc/Glo
2564   ///    local are:   pgx0, pgz0, ugx, ugz   OR   riemann fitters pars
2565   ///    global are:  tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
2566   /// return 0 if success
2567
2568   AliDebug(3,Form("lpar = %g %g %g %g %g\ngpar= %g %g %g %g %g %g\n",lpar[0],lpar[1],lpar[2],lpar[3],lpar[4],gpar[0],gpar[1],gpar[2],gpar[3],gpar[4],gpar[5]));
2569   AliDebug(3,Form("deltalpar = %g %g %g %g %g\n",lpar[0]-fLocalInitParam[0],lpar[1]-fLocalInitParam[1],lpar[2]-fLocalInitParam[2],lpar[3]-fLocalInitParam[3],lpar[4]-fLocalInitParam[4]));
2570
2571
2572   // prepare the TGeoHMatrix
2573   TGeoHMatrix *tempHMat = fCurrentModule->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar,
2574                                                                            !fUseGlobalDelta);
2575   if (!tempHMat) return -1;
2576
2577   Double_t v0g[3]; // vector with straight line direction in global coord.
2578   Double_t p0g[3]; // point of the straight line (glo)
2579
2580   if (fBOn) { // B FIELD!
2581     AliTrackPoint prf;
2582     for (int ip=0; ip<5; ip++)
2583       fRieman->SetParam(ip,lpar[ip]);
2584
2585     if (!fRieman->GetPCA(fCluster,prf))  {
2586       AliInfo(Form("error in GetPCA for point %d",fCluster.GetVolumeID()));
2587       return -3;
2588     }
2589     // now determine straight line passing tangent to fit curve at prf
2590     // ugx = dX/dY_glo = DeltaX/DeltaY_glo
2591     // mo' P1=(X,Y,Z)_glo_prf
2592     //       => (x,y,Z)_trk_prf ruotando di alpha...
2593     Double_t alpha=fRieman->GetAlpha();
2594     Double_t x1g = prf.GetX();
2595     Double_t y1g = prf.GetY();
2596     Double_t z1g = prf.GetZ();
2597     Double_t x1t =  x1g*TMath::Cos(alpha) + y1g*TMath::Sin(alpha);
2598     Double_t y1t = -x1g*TMath::Sin(alpha) + y1g*TMath::Cos(alpha);
2599     Double_t z1t =  z1g;
2600
2601     Double_t x2t = x1t+1.0;
2602     Double_t y2t = y1t+fRieman->GetDYat(x1t);
2603     Double_t z2t = z1t+fRieman->GetDZat(x1t);
2604     Double_t x2g =  x2t*TMath::Cos(alpha) - y2t*TMath::Sin(alpha);
2605     Double_t y2g =  x2t*TMath::Sin(alpha) + y2t*TMath::Cos(alpha);
2606     Double_t z2g =  z2t;
2607
2608     AliDebug(3,Form("Riemann frame:  fAlpha = %+f  =  %+f  ",alpha,alpha*180./TMath::Pi()));
2609     AliDebug(3,Form("   prf_glo=( %+f , %+f , %+f )  prf_rf=( %+f , %+f , %+f )\n", x1g,y1g,z1g, x1t,y1t,z1t));
2610     AliDebug(3,Form("   mov_glo=( %+f , %+f , %+f )      rf=( %+f , %+f , %+f )\n",x2g,y2g,z2g, x2t,y2t,z2t));
2611
2612     if (TMath::Abs(y2g-y1g)<1e-15) {
2613       AliInfo("DeltaY=0! Cannot proceed...");
2614       return -1;
2615     }
2616     // ugx,1,ugz
2617     v0g[0] = (x2g-x1g)/(y2g-y1g);
2618     v0g[1]=1.0;
2619     v0g[2] = (z2g-z1g)/(y2g-y1g);
2620
2621     // point: just keep prf
2622     p0g[0]=x1g;
2623     p0g[1]=y1g;
2624     p0g[2]=z1g;
2625   }
2626   else { // staight line
2627     // vector of initial straight line direction in glob. coord
2628     v0g[0]=lpar[2];
2629     v0g[1]=1.0;
2630     v0g[2]=lpar[3];
2631
2632     // intercept in yg=0 plane in glob coord
2633     p0g[0]=lpar[0];
2634     p0g[1]=0.0;
2635     p0g[2]=lpar[1];
2636   }
2637   AliDebug(3,Form("Line vector: ( %+f , %+f , %+f )  point:( %+f , %+f , %+f )\n",v0g[0],v0g[1],v0g[2],p0g[0],p0g[1],p0g[2]));
2638
2639   // same in local coord.
2640   Double_t p0l[3],v0l[3];
2641   tempHMat->MasterToLocalVect(v0g,v0l);
2642   tempHMat->MasterToLocal(p0g,p0l);
2643
2644   if (TMath::Abs(v0l[1])<1e-15) {
2645     AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
2646     return -1;
2647   }
2648
2649   // local intersection point
2650   fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
2651   fPintLoc[1] = 0;
2652   fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
2653
2654   // global intersection point
2655   tempHMat->LocalToMaster(fPintLoc,fPintGlo);
2656   AliDebug(3,Form("Intesect. point: L( %+f , %+f , %+f )  G( %+f , %+f , %+f )\n",fPintLoc[0],fPintLoc[1],fPintLoc[2],fPintGlo[0],fPintGlo[1],fPintGlo[2]));
2657
2658   return 0;
2659 }
2660
2661 //________________________________________________________________________________________________________
2662 Int_t AliITSAlignMille2::CalcDerivatives(Int_t paridx, Bool_t islpar)
2663 {
2664   /// calculate numerically (ROOT's style) the derivatives for
2665   /// local X intersection  and local Z intersection
2666   /// parlist: local  (islpar=kTRUE)  pgx0, pgz0, ugx0, ugz0  OR riemann's params
2667   ///          global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
2668   /// return 0 if success
2669
2670   // copy initial parameters
2671   Double_t lpar[kNLocal];
2672   Double_t gpar[kNParCh];
2673   Double_t *derivative;
2674   for (Int_t i=0; i<kNLocal; i++) lpar[i]=fLocalInitParam[i];
2675   for (Int_t i=0; i<kNParCh; i++) gpar[i]=fModuleInitParam[i];
2676
2677   // trial with fixed dpar...
2678   Double_t dpar = 0.;
2679
2680   if (islpar) { // track parameters
2681     //dpar=fLocalInitParam[paridx]*0.001;
2682     // set minimum dpar
2683     derivative = fDerivativeLoc[paridx];
2684     if (!fBOn) {
2685       if (paridx<3) dpar=1.0e-4; // translations
2686       else dpar=1.0e-6; // direction
2687     }
2688     else { // B Field
2689       // pepo: proviamo con 1/1000, poi evenctually 1/100...
2690       Double_t dfrac=0.01;
2691       switch(paridx) {
2692       case 0:
2693         // RMS cosmics: 1e-4
2694         dpar = TMath::Max(1.0e-6,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2695         break;
2696       case 1:
2697         // RMS cosmics: 0.2
2698         dpar = TMath::Max(0.002,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2699         break;
2700       case 2:
2701         // RMS cosmics: 9
2702         dpar = TMath::Max(0.09,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2703         break;
2704       case 3:
2705         // RMS cosmics: 7
2706         dpar = TMath::Max(0.07,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2707         break;
2708       case 4:
2709         // RMS cosmics: 0.3
2710         dpar = TMath::Max(0.003,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2711         break;
2712       }
2713     }
2714   }
2715   else { // alignment global parameters
2716     derivative = fDerivativeGlo[paridx];
2717     //dpar=fModuleInitParam[paridx]*0.001;
2718     if (paridx<3) dpar=1.0e-4; // translations
2719     else dpar=1.0e-2; // angles
2720   }
2721
2722   AliDebug(3,Form("+++ using dpar=%g",dpar));
2723
2724   // calculate derivative ROOT's like:
2725   //  using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
2726   Double_t pintl1[3]; // f(x-h)
2727   Double_t pintl2[3]; // f(x-h/2)
2728   Double_t pintl3[3]; // f(x+h/2)
2729   Double_t pintl4[3]; // f(x+h)
2730
2731   // first values
2732   if (islpar) lpar[paridx] -= dpar;
2733   else gpar[paridx] -= dpar;
2734   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2735   for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
2736
2737   // second values
2738   if (islpar) lpar[paridx] += dpar/2;
2739   else gpar[paridx] += dpar/2;
2740   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2741   for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
2742
2743   // third values
2744   if (islpar) lpar[paridx] += dpar;
2745   else gpar[paridx] += dpar;
2746   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2747   for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
2748
2749   // fourth values
2750   if (islpar) lpar[paridx] += dpar/2;
2751   else gpar[paridx] += dpar/2;
2752   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2753   for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
2754
2755   Double_t h2 = 1./(2.*dpar);
2756   Double_t d0 = pintl4[0]-pintl1[0];
2757   Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
2758   derivative[0] = h2*(4*d2 - d0)/3.;
2759   if (TMath::Abs(derivative[0]) < 1.0e-9) derivative[0] = 0.0;
2760
2761   d0 = pintl4[2]-pintl1[2];
2762   d2 = 2.*(pintl3[2]-pintl2[2]);
2763   derivative[2] = h2*(4*d2 - d0)/3.;
2764   if (TMath::Abs(derivative[2]) < 1.0e-9) derivative[2]=0.0;
2765
2766   AliDebug(3,Form("\n+++ derivatives +++ \n"));
2767   AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",derivative[0]));
2768   AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",derivative[2]));
2769
2770   return 0;
2771 }
2772
2773 //________________________________________________________________________________________________________
2774 Int_t AliITSAlignMille2::AddLocalEquation(Mille2Data &m)
2775 {
2776   /// Define local equation for current cluster in X and Z coor.
2777   /// and store them to memory
2778   /// return -1 in case of failure to build some equation
2779   ///         0 if no free global parameters were found but local eq is built
2780   ///         1 if both local and global eqs are built
2781   //
2782   // store first intersection point
2783   if (CalcIntersectionPoint(fLocalInitParam, fModuleInitParam)) return -1;
2784   for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
2785
2786   AliDebug(2,Form("Intersect. point: L( %+f , %+f , %+f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
2787
2788   // calculate local derivatives numerically
2789   Bool_t zeroX = kTRUE;
2790   Bool_t zeroZ = kTRUE;
2791   //
2792   for (Int_t i=0; i<fNLocal; i++) {
2793     if (CalcDerivatives(i,kTRUE)) return -1;
2794     m.fDerLoc[i][kX] = fDerivativeLoc[i][0];
2795     m.fDerLoc[i][kZ] = fDerivativeLoc[i][2];
2796     if (zeroX) zeroX = IsZero(fDerivativeLoc[i][0]);
2797     if (zeroZ) zeroZ = IsZero(fDerivativeLoc[i][2]);
2798   }
2799   //  for (Int_t i=0; i<fNLocal; i++) AliDebug(2,Form("Local parameter %d - dXdpar = %g  - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
2800   //
2801   if (zeroX) {AliInfo("Skipping: zero local X derivatives!"); return -1;}
2802   if (zeroZ) {AliInfo("Skipping: zero local Z derivatives!"); return -1;}
2803   //
2804   int status = 0;
2805   int ifill = 0;
2806   //
2807   AliITSAlignMille2Module* endModule = fCurrentModule;
2808   //
2809   zeroX = zeroZ = kTRUE;
2810   Bool_t dfDone[kNParCh];
2811   for (int i=kNParCh;i--;) dfDone[i] = kFALSE;
2812   m.fNModFilled = 0;
2813   //
2814   // special block for SDD derivatives
2815   Double_t jacobian[kNParChGeom];
2816   Int_t nmodTested = 0;
2817   //
2818   do {
2819     if (fCurrentModule->GetNParFree()==0) continue;
2820     nmodTested++;
2821     for (Int_t i=0; i<kNParChGeom; i++) {   // common for all sensors: derivatives over geom params
2822       //
2823       if (!fUseGlobalDelta) dfDone[i] = kFALSE; // for global deltas the derivatives at diff. levels are different
2824       if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2825       if (!dfDone[i]) {
2826         if (CalcDerivatives(i,kFALSE)) return -1;
2827         else {
2828           dfDone[i] = kTRUE;
2829           if (zeroX) zeroX = IsZero(fDerivativeGlo[i][0]);
2830           if (zeroZ) zeroZ = IsZero(fDerivativeGlo[i][2]);
2831         }
2832       }
2833       //
2834       m.fDerGlo[ifill][kX] = fDerivativeGlo[i][0];
2835       m.fDerGlo[ifill][kZ] = fDerivativeGlo[i][2];
2836       m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
2837     }
2838     //
2839     // specific for special sensors
2840     Int_t sddLR = -1;
2841     if ( fCurrentModule->IsSDD() &&
2842          (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0  ||
2843           //      fCurrentModule->GetParOffset(sddLR = fMeasLoc[kX]>0 ?
2844           fCurrentModule->GetParOffset(sddLR = GetVDriftSDD()>0 ?
2845                                        AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR)>=0)
2846          ) {
2847       //
2848       // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
2849       // where V0 and T are the nominal drift velocity, time and time0
2850       // and the dT0 and dV are the corrections:
2851       // dX/dT0 = dX/dxloc * dxloc/dT0 = dX/dxloc * V0
2852       // dX/dV  = dX/dxloc * dxloc/dV =  dX/dxloc * (T-T0)
2853       // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
2854       //
2855       if (!dfDone[AliITSAlignMille2Module::kDOFT0] ||  !dfDone[sddLR]) {
2856         //
2857         double dXdxlocsens=0., dZdxlocsens=0.;
2858         //
2859         // if the current module is the sensor itself and we work with local params, then
2860         // we can directly take dX/dxloc_sens dZ/dxloc_sens
2861         if (!fUseGlobalDelta && fCurrentModule->GetVolumeID()==fCluster.GetVolumeID()) {
2862           if (!dfDone[AliITSAlignMille2Module::kDOFTX]) {
2863             CalcDerivatives(AliITSAlignMille2Module::kDOFTX,kFALSE);
2864             dfDone[AliITSAlignMille2Module::kDOFTX] = kTRUE;
2865           }
2866           dXdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][0];
2867           dZdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][2];
2868         }
2869         //
2870         else { // need to perform some transformations
2871           // fetch the jacobian of the transformation from the sensors local frame to the frame
2872           // where the parameters are defined:
2873           // Global: dX/dxloc_sens = dX/dxgl*dxgl/dxloc_sens + ...dX/dphigl*dphigl/dxloc_sens
2874           if (fUseGlobalDelta) fCurrentModule->CalcDerivGloLoc(fCluster.GetVolumeID(),
2875                                                                AliITSAlignMille2Module::kDOFTX, jacobian);
2876           // Local:  dX/dxloc_sens = dX/dxcurr*dxcurr/dxloc_sens +..+dX/dphicurr * dphicurr/dxloc_sens
2877           else                 fCurrentModule->CalcDerivCurLoc(fCluster.GetVolumeID(),
2878                                                                AliITSAlignMille2Module::kDOFTX, jacobian);
2879           //
2880           for (int j=0;j<kNParChGeom;j++) {
2881             // need global derivative even if the j-th param is locked
2882             if (!dfDone[j]) {CalcDerivatives(j,kFALSE); dfDone[j] = kTRUE;}
2883             dXdxlocsens += fDerivativeGlo[j][0] * jacobian[j];
2884             dZdxlocsens += fDerivativeGlo[j][2] * jacobian[j];
2885           }
2886         }
2887         //
2888         if (zeroX) zeroX = IsZero(dXdxlocsens);
2889         if (zeroZ) zeroZ = IsZero(dZdxlocsens);
2890         //
2891         double vdrift = GetVDriftSDD();
2892         double tdrift = GetTDriftSDD();
2893         //
2894         fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0] = dXdxlocsens*vdrift;
2895         fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2] = dZdxlocsens*vdrift;
2896         dfDone[AliITSAlignMille2Module::kDOFT0] = kTRUE;
2897         //
2898         double mltCorr = fIsSDDVDriftMult ? TMath::Abs(vdrift) : 1;
2899         fDerivativeGlo[sddLR][0] = -dXdxlocsens*mltCorr*TMath::Sign(tdrift,vdrift);
2900         fDerivativeGlo[sddLR][2] = -dZdxlocsens*mltCorr*TMath::Sign(tdrift,vdrift);
2901         dfDone[sddLR] = kTRUE;
2902         //
2903       }
2904       //
2905       if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
2906         m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0];
2907         m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2];
2908         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);
2909       }
2910       //
2911       if (fCurrentModule->GetParOffset(sddLR)>=0) {
2912         m.fDerGlo[ifill][kX] = fDerivativeGlo[sddLR][0];
2913         m.fDerGlo[ifill][kZ] = fDerivativeGlo[sddLR][2];
2914         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(sddLR);
2915       }
2916     }
2917     //
2918     m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
2919   } while( (fCurrentModule=fCurrentModule->GetParent()) );
2920   //
2921   if (nmodTested>0 && zeroX) {AliInfo("Skipping: zero global X derivatives!");return -1;}
2922   if (nmodTested>0 && zeroZ) {AliInfo("Skipping: zero global Z derivatives!");return -1;}
2923   //
2924   // ok, can copy to m
2925   AliDebug(2,Form("Adding local equation X with fMeas=%.6f  and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
2926   m.fMeas[kX] = fMeasLoc[0]-fPintLoc0[0];
2927   m.fSigma[kX] = fSigmaLoc[0];
2928   //
2929   AliDebug(2,Form("Adding local equation Z with fMeas=%.6f  and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
2930   m.fMeas[kZ] = fMeasLoc[2]-fPintLoc0[2];
2931   m.fSigma[kZ] = fSigmaLoc[2];
2932   //
2933   m.fNGlobFilled = ifill;
2934   fCurrentModule = endModule;
2935   //
2936   status += Int_t(!zeroX && !zeroZ); // 0 - only locals, 1 locals + globals
2937   return status;
2938 }
2939
2940 //________________________________________________________________________________________________________
2941 Int_t AliITSAlignMille2::AddLocalEquationTPA(Mille2Data &m)
2942 {
2943   /// Define local equation for current cluster in X Y and Z coor.
2944   /// and store them to memory
2945   /// return -1 in case of failure to build some equation
2946   ///         0 if no free global parameters were found but local eq is built
2947   ///         1 if both local and global eqs are built
2948   //
2949   static int cnt = 0;
2950   Bool_t dbg = kFALSE;//kTRUE;
2951   if (++cnt>100000) dbg = kFALSE;
2952
2953   int curpoint = fCluster.GetUniqueID()-1;
2954   TGeoHMatrix *tempHMat = GetSensorCurrMatrixSID(fCurrentSensID);// fCurrentModule->GetSensitiveVolumeMatrix(fCluster.GetVolumeID());
2955   //
2956   fTPAFitter->GetDResDParams(&fDerivativeLoc[0][0], curpoint);    // resid. derivatives over the track parameters
2957   for (Int_t i=fNLocal; i--;) tempHMat->MasterToLocalVect(fDerivativeLoc[i],m.fDerLoc[i]);
2958   //
2959   int status = 0;
2960   // derivatives over the global parameters ---------------------------------------->>>
2961   Double_t dGL[3];     // derivative of global position vs local X (for SDD)
2962   Double_t dRdP[3][3]; // derivative of local residuals vs local position
2963   Double_t dPdG[AliITSAlignMille2Module::kMaxParGeom][3]; // derivatives of local position vs global params
2964   fTPAFitter->GetDResDPos(&fDerivativeGlo[0][0], curpoint);
2965   if (fCurrentSensID!=kVtxSensID) for (int i=3;i--;) tempHMat->MasterToLocalVect(fDerivativeGlo[i],dRdP[i]);
2966   else for (int i=3;i--;) for (int j=3;j--;) dRdP[i][j] = fDerivativeGlo[i][j];  // vertex constraint is in lab
2967   //
2968   if (dbg) {
2969     printf("\nCurrentMod: %s Sens:%d\n",fCurrentModule->GetName(),fCurrentSensID); //RRR
2970     printf("Module Matrix: ");
2971     fCurrentModule->GetMatrix()->Print(); //RRR
2972     for (int i=0;i<3;i++) {
2973       printf("dRdP[M%d][resI] ",i); for (int j=0;j<3;j++) printf(":[%d] %+.3e ",j,dRdP[i][j]); printf("\n");
2974     }//RRR
2975     printf("Sensor Matrix: "); tempHMat->Print();
2976   }
2977   UInt_t ifill=0, dfDone = 0;
2978   m.fNModFilled = 0;
2979   //
2980   AliITSAlignMille2Module* endModule = fCurrentModule;
2981   //
2982   m.fModuleID[0] = fCurrentModule->GetUniqueID(); // always register id of the base module, even if it has no DOF
2983   //
2984   do {
2985     if (fCurrentModule->GetNParFree()==0) continue;
2986     status = 1;
2987     if (!fUseGlobalDelta) dfDone = 0; // for local deltas the derivatives at diff. levels are different
2988     Bool_t jacobOK = kFALSE;
2989     //
2990     for (Int_t i=0; i<kNParChGeom; i++) {              // common for all sensors: derivatives over geom params
2991       if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2992       //
2993       if (!TestWordBit(dfDone,i)) {                    // need to calculate new derivative
2994         if (!jacobOK) {
2995           if (fCurrentSensID!=kVtxSensID) {
2996             fCurrentModule->CalcDerivDPosDPar(fCluster.GetVolumeID(),fMeasLoc,&dPdG[0][0]);
2997             if (dbg) {
2998               for (int i1=0;i1<3;i1++) {
2999                 printf("Jacob:dPdG[gpar%d][Mj]",i1); for (int j1=0;j1<3;j1++) printf(":[%d] %+.3e ",j1,dPdG[i1][j1]); printf("\n");//RRR
3000               }
3001             }
3002           }
3003           else {
3004             // this is a vertex constraint: only lateral shifts are allowed, no rotations
3005             for (int ip=AliITSAlignMille2Module::kMaxParGeom;ip--;) for (int jp=3;jp--;) dPdG[ip][jp] = (ip==jp) ? 1:0;
3006           }
3007           jacobOK = kTRUE;
3008         }
3009         // dRes_j/dGlo_i = \sum_{k=1:3}  dRes_j/dPos_k * dPos_k/dGlo_i
3010         fDerivativeGlo[i][kX] = dRdP[kX][kX]*dPdG[i][kX] + dRdP[kY][kX]*dPdG[i][kY] + dRdP[kZ][kX]*dPdG[i][kZ];
3011         fDerivativeGlo[i][kY] = dRdP[kX][kY]*dPdG[i][kX] + dRdP[kY][kY]*dPdG[i][kY] + dRdP[kZ][kY]*dPdG[i][kZ];
3012         fDerivativeGlo[i][kZ] = dRdP[kX][kZ]*dPdG[i][kX] + dRdP[kY][kZ]*dPdG[i][kY] + dRdP[kZ][kZ]*dPdG[i][kZ];
3013         SetWordBit(dfDone,i);
3014       }
3015       //
3016       if (dbg) {
3017         printf("Level %s DGlob[par%d][resJ] ",fCurrentModule->GetName(),i); //RRR
3018         for (int k=0;k<3;k++) printf(":[%d] %+.3e ",k, fDerivativeGlo[i][k]); printf("\n");//RRR
3019       }
3020       m.fDerGlo[ifill][kX] = fDerivativeGlo[i][kX];
3021       m.fDerGlo[ifill][kY] = fDerivativeGlo[i][kY];
3022       m.fDerGlo[ifill][kZ] = fDerivativeGlo[i][kZ];
3023       m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
3024       //
3025     }
3026     //
3027     if ( fCurrentModule->IsSDD() ) {     // specific for SDD
3028       //
3029       // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
3030       // where V0 and T are the nominal drift velocity, time and time0
3031       // and the dT0 and dV are the corrections:
3032       // drloc_i/dT0 = sum_j drloc_i/dMeasGlo_j * dMeasGlo_j/dT0 =
3033       //             = sum_j drloc_i/dMeasGlo_j sum_k dMeasGlo_j/dMeasLoc_k * dMeasLoc_k/dT0
3034       //             = sum_j drloc_i/dMeasGlo_j dMeasGlo_j/dMeasLoc_X * V0
3035       //
3036       // drloc_i/dV0 = sum_j drloc_i/dMeasGlo_j * dMeasGlo_j/dV0 =
3037       //             = sum_j drloc_i/dMeasGlo_j sum_k dMeasGlo_j/dMeasLoc_k * dMeasLoc_k/dV0
3038       //             = sum_j drloc_i/dMeasGlo_j dMeasGlo_j/dMeasLoc_X * T0
3039
3040       // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
3041       //
3042       Bool_t jacOK = kFALSE;
3043       //Int_t sddLR = fMeasLoc[kX]>0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR;
3044       Int_t sddLR = GetVDriftSDD()>0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR;
3045       if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
3046         if (!TestWordBit(dfDone, AliITSAlignMille2Module::kDOFT0)) {
3047           double vdrift = GetVDriftSDD();
3048           JacobianPosGloLoc(kX,dGL);
3049           jacOK = kTRUE;
3050           fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kX] =
3051             vdrift*(dRdP[kX][kX]*dGL[kX] + dRdP[kY][kX]*dGL[kY] + dRdP[kZ][kX]*dGL[kZ]);
3052           fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kY] =
3053             vdrift*(dRdP[kX][kY]*dGL[kX] + dRdP[kY][kY]*dGL[kY] + dRdP[kZ][kY]*dGL[kZ]);
3054           fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kZ] =
3055             vdrift*(dRdP[kX][kZ]*dGL[kX] + dRdP[kY][kZ]*dGL[kY] + dRdP[kZ][kZ]*dGL[kZ]);
3056           //
3057           SetWordBit(dfDone, AliITSAlignMille2Module::kDOFT0);
3058         }
3059         m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kX];
3060         m.fDerGlo[ifill][kY] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kY];
3061         m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kZ];
3062         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);
3063       }
3064       //
3065       if (fCurrentModule->GetParOffset(sddLR)>=0) {
3066         if (!TestWordBit(dfDone, sddLR)) {
3067           double tdrift = TMath::Sign(GetTDriftSDD(), GetVDriftSDD());
3068           double vdrift = fIsSDDVDriftMult ? TMath::Abs(GetVDriftSDD()) : 1;
3069           if (!jacOK) JacobianPosGloLoc(kX,dGL);
3070           fDerivativeGlo[sddLR][kX] =
3071             -tdrift*vdrift*(dRdP[kX][kX]*dGL[kX] + dRdP[kY][kX]*dGL[kY] + dRdP[kZ][kX]*dGL[kZ]);
3072           fDerivativeGlo[sddLR][kY] =
3073             -tdrift*vdrift*(dRdP[kX][kY]*dGL[kX] + dRdP[kY][kY]*dGL[kY] + dRdP[kZ][kY]*dGL[kZ]);
3074           fDerivativeGlo[sddLR][kZ] =
3075             -tdrift*vdrift*(dRdP[kX][kZ]*dGL[kX] + dRdP[kY][kZ]*dGL[kY] + dRdP[kZ][kZ]*dGL[kZ]);
3076           SetWordBit(dfDone, sddLR);
3077         }
3078         m.fDerGlo[ifill][kX] = fDerivativeGlo[sddLR][kX];
3079         m.fDerGlo[ifill][kY] = fDerivativeGlo[sddLR][kY];
3080         m.fDerGlo[ifill][kZ] = fDerivativeGlo[sddLR][kZ];
3081         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(sddLR);
3082       }
3083     }
3084     //
3085     m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
3086   } while( (fCurrentModule=fCurrentModule->GetParent()) );
3087   //
3088   // store first local residuals
3089   fTPAFitter->GetResiduals(fPintLoc , curpoint);       // lab residuals
3090   for (int i=3;i--;) fPintLoc[i] = -fPintLoc[i];
3091   if   (fCurrentSensID!=kVtxSensID) tempHMat->MasterToLocalVect(fPintLoc,m.fMeas);       // local residuals
3092   else for (int i=3;i--;) m.fMeas[i] = fPintLoc[i];
3093   if (dbg) {
3094     printf("res(meas-loc) "); for (int k=0;k<3;k++) printf(":[%d] %+.3e ",k,m.fMeas[k]); printf("\n");
3095     printf("Fin:%s %+e %+e\n",endModule->GetName(), fDerivativeGlo[kZ][kZ], fPintLoc[kZ]);
3096   }//RRR
3097   m.fSigma[kX] = fSigmaLoc[kX];
3098   m.fSigma[kY] = fSigmaLoc[kY];
3099   m.fSigma[kZ] = fSigmaLoc[kZ];
3100   //
3101   m.fNGlobFilled = ifill;
3102   fCurrentModule = endModule;
3103   //
3104   return status;
3105 }
3106
3107 //________________________________________________________________________________________________________
3108 void AliITSAlignMille2::SetLocalEquations(const Mille2Data *marr, Int_t neq)
3109 {
3110   /// Set local equations with data stored in m
3111   /// return 0 if success
3112   //
3113   for (Int_t j=0; j<neq; j++) {
3114     //
3115     const Mille2Data &m = marr[j];
3116     //
3117     Bool_t filled = kFALSE;
3118     for (int ic=3;ic--;) {
3119       // for the diamond point (if any) the Y residual is accounted
3120       if (ic==kY && !fUseLocalYErr && !(m.fModuleID[0]==fDiamondModID)) continue;
3121       AliDebug(2,Form("setting local equation %c with fMeas=%.6f  and fSigma=%.6f",fgkXYZ[ic],m.fMeas[ic], m.fSigma[ic]));
3122       Int_t nzero = 0;
3123       for (int i=fNLocal; i--;) nzero += SetLocalDerivative(i,m.fDerLoc[i][ic] );
3124       if (nzero==fNLocal) {
3125         AliInfo(Form("Skipping %c residual due to the zero derivatives!",fgkXYZ[ic]));
3126         continue;
3127       }
3128       for (int i=m.fNGlobFilled;i--;) SetGlobalDerivative( m.fParMilleID[i] , m.fDerGlo[i][ic] );
3129       fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.fMeas[ic], m.fSigma[ic]);
3130       filled = kTRUE;
3131       //
3132     }
3133     //
3134     if (filled) for (int i=m.fNModFilled;i--;) GetMilleModule(m.fModuleID[i])->IncNProcessedPoints();
3135   }
3136   //
3137   double wgh = 1;
3138   if (GetWeightPt() && fTPAFitter) {
3139     wgh = fTPAFitter->GetPt();
3140     if (wgh>10) wgh = 10.;
3141     if (wgh<0) wgh = fTPAFitter->IsTypeCosmics() ? 7 : 0.5;
3142     if (GetWeightPt()>0) wgh = TMath::Power(wgh,GetWeightPt());
3143   }
3144   fMillepede->SetRecordWeight(wgh*fTrackWeight);
3145   fMillepede->SetRecordRun(fRunID);
3146   //
3147 }
3148
3149 //________________________________________________________________________________________________________
3150 Int_t AliITSAlignMille2::GlobalFit()
3151 {
3152   /// Call global fit; Global parameters are stored in parameters
3153   if (!fIsMilleInit) Init();
3154   //
3155   ApplyPreConstraints();
3156   int res = fMillepede->GlobalFit();
3157   AliInfo(Form("%s fitting global parameters!",res ? "Done":"Failed"));
3158   if (res) {
3159     // fetch the parameters
3160     for (int imd=fNModules;imd--;) {
3161       AliITSAlignMille2Module* mod = GetMilleModule(imd);
3162       int nprocp = 0;
3163       for (int ip=mod->GetNParTot();ip--;) {
3164         int idp = mod->GetParOffset(ip);
3165         if (idp<0) continue;    // was not in the explicit fit
3166         mod->SetParVal(ip,fMillepede->GetFinalParam(idp));
3167         mod->SetParErr(ip,fMillepede->GetFinalError(idp));
3168         int np = fMillepede->GetProcessedPoints(idp);
3169         if (TMath::Abs(np)>TMath::Abs(nprocp)) nprocp = np;
3170       }
3171       if (!mod->GetNProcessedPoints()) mod->SetNProcessedPoints(nprocp);
3172     }
3173
3174   }
3175   ApplyPostConstraints();
3176   return res;
3177 }
3178
3179 //________________________________________________________________________________________________________
3180 void AliITSAlignMille2::PrintGlobalParameters()
3181 {
3182   /// Print global parameters
3183   if (!fIsMilleInit) {
3184     AliInfo("Millepede has not been initialized!");
3185     return;
3186   }
3187   fMillepede->PrintGlobalParameters();
3188 }
3189
3190 //________________________________________________________________________________________________________
3191 Int_t AliITSAlignMille2::LoadSuperModuleFile(const Char_t *sfile)
3192 {
3193   // load definitions of supermodules from a root file
3194   // return 0 if success
3195   AliInfo(Form("Loading SuperModule definitions from %s",sfile));
3196   TFile *smf=TFile::Open(sfile);
3197   if (!smf->IsOpen()) {
3198     AliInfo(Form("Cannot open supermodule file %s",sfile));
3199     return -1;
3200   }
3201
3202   TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules");
3203   if (!sma) {
3204     AliInfo(Form("Cannot find ITSMilleSuperModules array in file"));
3205     return -2;
3206   }
3207   Int_t nsma=sma->GetEntriesFast();
3208   AliInfo(Form("Array of SuperModules with %d entries\n",nsma));
3209   //
3210   // pepo200709
3211   Char_t st[2048];
3212   char symname[250];
3213   // end pepo200709
3214
3215   UShort_t volid;
3216   TGeoHMatrix m;
3217   //
3218   for (Int_t i=0; i<nsma; i++) {
3219     AliAlignObjParams *a = (AliAlignObjParams*)sma->UncheckedAt(i);
3220     volid=a->GetVolUID();
3221     strcpy(st,a->GetSymName());
3222     a->GetMatrix(m);
3223     //
3224     sscanf(st,"%s",symname);
3225     //
3226     // decode module list
3227     char *stp=strstr(st,"ModuleList:");
3228     if (!stp) return -3;
3229     stp += 11;
3230     int idx[2200];
3231     char spp[200]; int jp=0;
3232     char cl[20];
3233     strcpy(st,stp);
3234     int l=strlen(st);
3235     int j=0;
3236     int n=0;
3237     //
3238     while (j<=l) {
3239       if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) {
3240         spp[jp]=0;
3241         jp=0;
3242         if (strlen(spp)) {
3243           int k=strcspn(spp,"-");
3244           if (k<int(strlen(spp))) { // c'e' il -
3245             strcpy(cl,&(spp[k+1]));
3246             spp[k]=0;
3247             int ifrom=atoi(spp); int ito=atoi(cl);
3248             for (int b=ifrom; b<=ito; b++) {
3249               idx[n]=b;
3250               n++;
3251             }
3252           }
3253           else { // numerillo singolo
3254             idx[n]=atoi(spp);
3255             n++;
3256           }
3257         }
3258       }
3259       else {
3260         spp[jp]=st[j];
3261         jp++;
3262       }
3263       j++;
3264     }
3265     UShort_t volidsv[2198];
3266     for (j=0;j<n;j++) {
3267       volidsv[j]=AliITSAlignMille2Module::GetVolumeIDFromIndex(idx[j]);
3268       if (!volidsv[j]) {
3269         AliInfo(Form("Index %d not valid (range 0->%d)",idx[j],kMaxITSSensID));
3270         return -5;
3271       }
3272     }
3273     Int_t smindex=int(2198+volid-14336); // virtual index
3274     //
3275     fSuperModule.AddAtAndExpand(new AliITSAlignMille2Module(smindex,volid,symname,&m,n,volidsv),fNSuperModules);
3276     //
3277     fNSuperModules++;
3278   }
3279   //
3280   smf->Close();
3281   //
3282   return 0;
3283 }
3284
3285 //________________________________________________________________________________________________________
3286 void AliITSAlignMille2::ConstrainModuleSubUnitsMean(Int_t idm, Double_t val, UInt_t pattern)
3287 {
3288   // require that sum of modifications for the childs of this module is = val, i.e.
3289   // the internal corrections moves the module as a whole by fixed value (0 by default).
3290   // pattern is the bit pattern for the parameters to constrain
3291   //
3292   if (fIsMilleInit) {
3293     AliInfo("Millepede has been already initialized: no constrain may be added!");
3294     return;
3295   }
3296   if (!GetMilleModule(idm)->GetNChildren()) return;
3297   TString nm = "cstrSUMean";
3298   nm += GetNConstraints();
3299   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
3300                                                                       idm,val,pattern);
3301   cstr->SetConstraintID(GetNConstraints());
3302   fConstraints.Add(cstr);
3303 }
3304
3305 //________________________________________________________________________________________________________
3306 void AliITSAlignMille2::ConstrainModuleSubUnitsMedian(Int_t idm, Double_t val, UInt_t pattern)
3307 {
3308   // require that median of the modifications for the childs of this module is = val, i.e.
3309   // the internal corrections moves the module as a whole by fixed value (0 by default)
3310   // module the outliers.
3311   // pattern is the bit pattern for the parameters to constrain
3312   // The difference between the mean and the median will be transfered to the parent
3313   if (fIsMilleInit) {
3314     AliInfo("Millepede has been already initialized: no constrain may be added!");
3315     return;
3316   }
3317   if (!GetMilleModule(idm)->GetNChildren()) return;
3318   TString nm = "cstrSUMed";
3319   nm += GetNConstraints();
3320   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
3321                                                                       idm,val,pattern);
3322   cstr->SetConstraintID(GetNConstraints());
3323   fConstraints.Add(cstr);
3324 }
3325
3326 //________________________________________________________________________________________________________
3327 void AliITSAlignMille2::ConstrainOrphansMean(Double_t val, UInt_t pattern)
3328 {
3329   // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3330   // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3331   // pattern is the bit pattern for the parameters to constrain
3332   //
3333   if (fIsMilleInit) {
3334     AliInfo("Millepede has been already initialized: no constrain may be added!");
3335     return;
3336   }
3337   TString nm = "cstrOMean";
3338   nm += GetNConstraints();
3339   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
3340                                                                       -1,val,pattern);
3341   cstr->SetConstraintID(GetNConstraints());
3342   fConstraints.Add(cstr);
3343 }
3344
3345 //________________________________________________________________________________________________________
3346 void AliITSAlignMille2::ConstrainOrphansMedian(Double_t val, UInt_t pattern)
3347 {
3348   // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3349   // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3350   // pattern is the bit pattern for the parameters to constrain
3351   //
3352   if (fIsMilleInit) {
3353     AliInfo("Millepede has been already initialized: no constrain may be added!");
3354     return;
3355   }
3356   TString nm = "cstrOMed";
3357   nm += GetNConstraints();
3358   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
3359                                                                       -1,val,pattern);
3360   cstr->SetConstraintID(GetNConstraints());
3361   fConstraints.Add(cstr);
3362 }
3363
3364 //________________________________________________________________________________________________________
3365 void AliITSAlignMille2::ConstrainLocal(const Char_t* name,Double_t *parcf,Int_t npar,Double_t val,Double_t err)
3366 {
3367   // apply constraint on parameters in the local frame
3368   if (fIsMilleInit) {
3369     AliInfo("Millepede has been already initialized: no constrain may be added!");
3370     return;
3371   }
3372   AliITSAlignMille2ConstrArray *cstr = new AliITSAlignMille2ConstrArray(name,parcf,npar,val,err);
3373   cstr->SetConstraintID(GetNConstraints());
3374   fConstraints.Add(cstr);
3375 }
3376
3377 //________________________________________________________________________________________________________
3378 void AliITSAlignMille2::ApplyGaussianConstraint(const AliITSAlignMille2ConstrArray* cstr)
3379 {
3380   // apply the constraint on the local corrections of a list of modules
3381   int nmod = cstr->GetNModules();
3382   double jacobian[AliITSAlignMille2Module::kMaxParGeom][AliITSAlignMille2Module::kMaxParGeom];
3383   //
3384   // check if this not special SDDT0 constraint
3385   if (cstr->GetPattern()==BIT(AliITSAlignMille2Module::kDOFT0)) {
3386     for (int i=0;i<cstr->GetNModules()-1;i++) {
3387       AliITSAlignMille2Module *mdI = GetMilleModule(cstr->GetModuleID(i));
3388       if (!mdI->IsFreeDOF(AliITSAlignMille2Module::kDOFT0)) continue;
3389       for (int j=i+1;j<cstr->GetNModules();j++) {
3390         AliITSAlignMille2Module *mdJ = GetMilleModule(cstr->GetModuleID(j));
3391         if (!mdJ->IsFreeDOF(AliITSAlignMille2Module::kDOFT0)) continue;
3392         //
3393         ResetLocalEquation();
3394         fGlobalDerivatives[mdI->GetParOffset(AliITSAlignMille2Module::kDOFT0)] = 1;
3395         fGlobalDerivatives[mdJ->GetParOffset(AliITSAlignMille2Module::kDOFT0)] =-1;
3396         AddConstraint(fGlobalDerivatives, 0, 1.E-6);
3397       }
3398     }
3399     return;
3400   }
3401
3402   for (int imd=nmod;imd--;) {
3403     int modID = cstr->GetModuleID(imd);
3404     AliITSAlignMille2Module* mod = GetMilleModule(modID);
3405     ResetLocalEquation();
3406     int nadded = 0;
3407     double value = cstr->GetValue();
3408     double sigma = cstr->GetError();
3409     //
3410     // in case the reference (survey) deltas were imposed for Gaussian constraints
3411     // already accumulated corrections: they must be subtracted from the constraint value.
3412     if (IsConstraintWrtRef()) {
3413       //
3414       Double_t precal[AliITSAlignMille2Module::kMaxParTot];
3415       Double_t refcal[AliITSAlignMille2Module::kMaxParTot];
3416       for (int ip=AliITSAlignMille2Module::kMaxParTot;ip--;) {precal[ip]=0; refcal[ip] = 0.;}
3417       //
3418       // check if there was a reference delta provided for this module
3419       AliAlignObjParams* parref = GetConstrRefObject(mod->GetName());
3420       if (parref) parref->GetPars(refcal, refcal+3);    // found reference delta
3421       //
3422       // extract already applied local corrections for this module
3423       if (fPrealignment) {
3424         //
3425         AliAlignObjParams *preo = GetPrealignedObject(mod->GetName());
3426         if (preo) {
3427           TGeoHMatrix preMat,tmpMat = *mod->GetMatrix(); //  Delta_Glob * Delta_Glob_Par * M
3428           preo->GetMatrix(preMat);                       //  Delta_Glob
3429           preMat.MultiplyLeft( &tmpMat.Inverse() );      //  M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
3430           tmpMat.MultiplyLeft( &preMat );                //  (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
3431           AliAlignObjParams algob;
3432           algob.SetMatrix(tmpMat);
3433           algob.GetPars(precal,precal+3); // local corrections for geometry
3434         }
3435       }
3436       //
3437       // subtract the contribution to constraint from precalibration
3438       for (int ipar=cstr->GetNCoeffs();ipar--;) value += (refcal[ipar]-precal[ipar])*cstr->GetCoeff(ipar);
3439       //
3440     }
3441     //
3442     if (fUseGlobalDelta) mod->CalcDerivLocGlo(&jacobian[0][0]);
3443     //
3444     for (int ipar=cstr->GetNCoeffs();ipar--;) {
3445       double coef = cstr->GetCoeff(ipar);
3446       if (IsZero(coef)) continue;
3447       //
3448       if (!fUseGlobalDelta || ipar>= AliITSAlignMille2Module::kMaxParGeom) { //
3449         // we are working with local params or if the given param is not related to geometry,
3450         // apply the constraint directly
3451         int parPos = mod->GetParOffset(ipar);
3452         if (parPos<0) continue; // not in the fit
3453         fGlobalDerivatives[parPos] += coef;
3454         nadded++;
3455       }
3456       else { // we are working with global params, while the constraint is on local ones -> jacobian
3457         for (int jpar=AliITSAlignMille2Module::kMaxParGeom;jpar--;) {
3458           int parPos = mod->GetParOffset(jpar);
3459           if (parPos<0) continue;
3460           fGlobalDerivatives[parPos] += coef*jacobian[ipar][jpar];
3461           nadded++;
3462         }
3463       }
3464     }
3465     if (nadded) AddConstraint(fGlobalDerivatives, value, sigma);
3466   }
3467   //
3468 }
3469
3470 //________________________________________________________________________________________________________
3471 void AliITSAlignMille2::ApplyPreConstraints()
3472 {
3473   // apply constriants which cannot be imposed after the fit
3474   int nconstr = GetNConstraints();
3475   for (int i=0;i<nconstr;i++) {
3476     AliITSAlignMille2Constraint* cstr = GetConstraint(i);
3477     //
3478     if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) {
3479       ApplyGaussianConstraint( (AliITSAlignMille2ConstrArray*)cstr);
3480       continue;
3481     }
3482     //
3483     if (cstr->GetType() == AliITSAlignMille2Constraint::kTypeMedian) continue; // post type constraint
3484     //
3485     if (!fUseGlobalDelta) continue; // mean/med constraints must be applied to global deltas
3486     // apply constraint on the mean's before the fit
3487     int imd = cstr->GetModuleID();
3488     if (imd>=0) {
3489       AliITSAlignMille2Module* mod = GetMilleModule(imd);
3490       UInt_t pattern = 0;
3491       for (int ipar=mod->GetNParTot();ipar--;) {
3492         if (!cstr->IncludesParam(ipar)) continue;
3493         if (mod->GetParOffset(ipar)<0) continue; // parameter is not in the explicit fit -> post constraint
3494         pattern |= 0x1<<ipar;
3495         cstr->SetApplied(ipar);
3496       }
3497       ConstrainModuleSubUnits(imd,cstr->GetValue(),pattern);
3498       //
3499     }
3500     else if (!PseudoParentsAllowed()) {
3501       ConstrainOrphans(cstr->GetValue(),(UInt_t)cstr->GetPattern());
3502       cstr->SetApplied(-1);
3503     }
3504   }
3505   //
3506   // do we need to tie the SDD left/right VDrift corrections
3507   for (int imd=0;imd<fNModules;imd++) {
3508     AliITSAlignMille2Module* mod = GetMilleModule(imd);
3509     if (mod->IsSDD() && mod->IsVDriftLRSame()) TieSDDVDriftsLR(mod);
3510   }
3511   //
3512 }
3513
3514 //________________________________________________________________________________________________________
3515 void AliITSAlignMille2::ApplyPostConstraints()
3516 {
3517   // apply constraints which can be imposed after the fit
3518   int nconstr = GetNConstraints();
3519   Bool_t convGlo      = kFALSE;
3520   // check if there is something to do
3521   int ntodo = 0;
3522   for (int i=0;i<nconstr;i++) {
3523     AliITSAlignMille2Constraint* cstr = GetConstraint(i);
3524     if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
3525     if (cstr->GetRemainingPattern() == 0) continue;
3526     ntodo++;
3527   }
3528   if (!ntodo) return;
3529   //
3530   if (!fUseGlobalDelta) { // need to convert to global params
3531     ConvertParamsToGlobal();
3532     convGlo = kTRUE;
3533   }
3534   //
3535   for (int i=0;i<nconstr;i++) {
3536     AliITSAlignMille2Constraint* cstr = GetConstraint(i);
3537     if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
3538     //
3539     int imd = cstr->GetModuleID();
3540     //
3541     if (imd>=0) {
3542       AliITSAlignMille2Module* mod = GetMilleModule(imd);
3543       if (mod->IsNotInConf()) continue;
3544       UInt_t pattern = 0;
3545       for (int ipar=mod->GetNParTot();ipar--;) {
3546         if (cstr->IsApplied(ipar))      continue;
3547         if (!cstr->IncludesParam(ipar)) continue;
3548         if (!mod->IsFreeDOF(ipar))      continue; // parameter is fixed, will not apply constraint
3549         pattern |= 0x1<<ipar;
3550         cstr->SetApplied(ipar);
3551       }
3552       if (pattern) PostConstrainModuleSubUnits(cstr->GetType(),cstr->GetModuleID(),cstr->GetValue(),pattern);
3553       //
3554     }
3555     else if (PseudoParentsAllowed()) {
3556       UInt_t pattern = (UInt_t)cstr->GetRemainingPattern();
3557       PostConstrainOrphans(cstr->GetType(),cstr->GetValue(),pattern);
3558       cstr->SetApplied(-1);
3559     }
3560   }
3561   // if there was a conversion, rewind it
3562   if (convGlo) ConvertParamsToLocal();
3563   //
3564 }
3565
3566 //________________________________________________________________________________________________________
3567 void AliITSAlignMille2::ConstrainModuleSubUnits(Int_t idm, Double_t val, UInt_t pattern)
3568 {
3569   // require that sum of modifications for the childs of this module is = val, i.e.
3570   // the internal corrections moves the module as a whole by fixed value (0 by default).
3571   // pattern is the bit pattern for the parameters to constrain
3572   //
3573   //
3574   AliITSAlignMille2Module* mod = GetMilleModule(idm);
3575   //
3576   for (int ip=0;ip<kNParCh;ip++) {
3577     if ( !((pattern>>ip)&0x1) /*|| !parent->IsFreeDOF(ip)*/) continue;
3578     ResetLocalEquation();
3579     int nadd = 0;
3580     for (int ich=mod->GetNChildren();ich--;) {
3581       int idpar = ((AliITSAlignMille2Module*)mod->GetChild(ich))->GetParOffset(ip);
3582       if (idpar<0) continue;
3583       fGlobalDerivatives[idpar] = 1.0;
3584       nadd++;
3585     }
3586     //
3587     if (nadd>0) {
3588       AddConstraint(fGlobalDerivatives,val);
3589       AliInfo(Form("Constrained param %d for %d submodules of module #%d: %s",ip,nadd,idm,mod->GetName()));
3590     }
3591   }
3592   //
3593 }
3594
3595 //________________________________________________________________________________________________________
3596 void AliITSAlignMille2::ConstrainOrphans(Double_t val, UInt_t pattern)
3597 {
3598   // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3599   // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3600   // pattern is the bit pattern for the parameters to constrain
3601   //
3602   for (int ip=0;ip<kNParCh;ip++) {
3603     //
3604     if ( !((pattern>>ip)&0x1) ) continue;
3605     ResetLocalEquation();
3606     int nadd = 0;
3607     for (int imd=fNModules;imd--;) {
3608       AliITSAlignMille2Module* mod = GetMilleModule(imd);
3609       if (mod->IsNotInConf()) continue; // dummy module
3610       AliITSAlignMille2Module* par = mod->GetParent();
3611       while (par && par->IsNotInConf() ) par = par->GetParent(); // use only decalred parents
3612       if (par) continue; // this is not an orphan
3613       int idpar = mod->GetParOffset(ip);
3614       if (idpar<0) continue;
3615       fGlobalDerivatives[idpar] = 1.0;
3616       nadd++;
3617     }
3618     if (nadd>0) {
3619       AddConstraint(fGlobalDerivatives,val);
3620       AliInfo(Form("Constrained param %d for %d orphan modules",ip,nadd));
3621     }
3622   }
3623   //
3624   //
3625 }
3626
3627 //________________________________________________________________________________________________________
3628 void AliITSAlignMille2::PostConstrainModuleSubUnits(Int_t type,Int_t idm, Double_t val, UInt_t pattern)
3629 {
3630   // require that median or mean of the modifications for the childs of this module is = val, i.e.
3631   // the internal corrections moves the module as a whole by fixed value (0 by default)
3632   // module the outliers.
3633   // pattern is the bit pattern for the parameters to constrain
3634   // The difference between the mean and the median will be transfered to the parent
3635   //
3636   AliITSAlignMille2Module* parent = GetMilleModule(idm);
3637   int nc = parent->GetNChildren();
3638   //
3639   double *tmpArr = new double[nc];
3640   //
3641   for (int ip=0;ip<kNParCh;ip++) {
3642     int npc = 0;
3643     if ( !((pattern>>ip)&0x1) || !parent->IsFreeDOF(ip)) continue;
3644     // compute the mean and median of the deltas
3645     int nfree = 0;
3646     for (int ich=nc;ich--;) {
3647       AliITSAlignMille2Module* child = parent->GetChild(ich);
3648       //      if (!child->IsFreeDOF(ip)) continue;
3649       tmpArr[nfree++] = child->GetParVal(ip);
3650     }
3651     double median=0,mean=0;
3652     for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
3653       mean += tmpArr[ic0];
3654       for (int ic1=ic0+1;ic1<nfree;ic1++)
3655         if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
3656     }
3657     //
3658     int kmed = nfree/2;
3659     median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
3660     if (nfree>0) mean /= nfree;
3661     //
3662     double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
3663     //
3664     for (int ich=nc;ich--;) {
3665       AliITSAlignMille2Module* child = parent->GetChild(ich);
3666       //    if (!child->IsFreeDOF(ip)) continue;
3667       child->SetParVal(ip, child->GetParVal(ip) + shift);
3668       npc++;
3669     }
3670     //
3671     parent->SetParVal(ip, parent->GetParVal(ip) - shift);
3672     AliInfo(Form("%s constraint: added %+f shift to param[%d] of %d children of module %d: %s",
3673                  type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
3674                  ip,npc,idm,parent->GetName()));
3675   }
3676   delete[] tmpArr;
3677   //
3678   //
3679 }
3680
3681 //________________________________________________________________________________________________________
3682 void AliITSAlignMille2::PostConstrainOrphans(Int_t type,Double_t val, UInt_t pattern)
3683 {
3684   // require that median or mean of modifications for the supermodules which have no parents is = val, i.e.
3685   // the corrections moves the whole setup by fixed value (0 by default).
3686   // pattern is the bit pattern for the parameters to constrain
3687   //
3688   int nc = fNModules;
3689   //
3690   int norph = 0;
3691   for (int ich=nc;ich--;) {
3692     AliITSAlignMille2Module *par= GetMilleModule(ich)->GetParent();
3693     while (par && par->IsNotInConf()) par = par->GetParent(); // use only decalred parents
3694     if (!par) norph ++;
3695   }
3696   //
3697   if (!norph) return;
3698   double *tmpArr = new double[norph];
3699   for (int i=norph;i--;) tmpArr[i] = 0;
3700   //
3701   for (int ip=0;ip<kNParCh;ip++) {
3702     int npc = 0;
3703     if ( !((pattern>>ip)&0x1)) continue;
3704     // compute the mean and median of the deltas
3705     int nfree = 0;
3706     for (int ich=nc;ich--;) {
3707       AliITSAlignMille2Module* child = GetMilleModule(ich);
3708       if (child->IsNotInConf()) continue; // dummy module
3709       //      if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
3710       AliITSAlignMille2Module* par = child->GetParent();
3711       while (par && par->IsNotInConf()) par = par->GetParent(); // count only declared parents
3712       if (par) continue;
3713       tmpArr[nfree++] = child->GetParVal(ip);
3714     }
3715     double median=0,mean=0;
3716     for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
3717       mean += tmpArr[ic0];
3718       for (int ic1=ic0+1;ic1<nfree;ic1++)
3719         if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
3720     }
3721     //
3722     int kmed = nfree/2;
3723     median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
3724     if (nfree>0) mean /= nfree;
3725     //
3726     double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
3727     //
3728     for (int ich=nc;ich--;) {
3729       AliITSAlignMille2Module* child = GetMilleModule(ich);
3730       if (child->IsNotInConf()) continue; // dummy module
3731       //      if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
3732       AliITSAlignMille2Module* par = child->GetParent();
3733       while (par && par->IsNotInConf()) par = par->GetParent(); // count only declared parents
3734       if (par) continue;
3735       child->SetParVal(ip, child->GetParVal(ip) + shift);
3736       npc++;
3737     }
3738     //
3739     AliInfo(Form("%s constraint: added %+f shift to param[%d] of %d orphan modules",
3740                  type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
3741                  ip,npc));
3742   }
3743   delete[] tmpArr;
3744   //
3745 }
3746
3747 //________________________________________________________________________________________________________
3748 Bool_t AliITSAlignMille2::IsParModConstrained(const AliITSAlignMille2Module* mod,Int_t par, Bool_t &meanmed, Bool_t &gaussian) const
3749 {
3750   // check if par of the module participates in some constraint, and set the flag for their types
3751   meanmed = gaussian = kFALSE;
3752   //
3753   if ( mod->IsParConstrained(par) ) gaussian = kTRUE;     // direct constraint on this param
3754   //
3755   for (int icstr=GetNConstraints();icstr--;) {
3756     AliITSAlignMille2Constraint* cstr = GetConstraint(icstr);
3757     //
3758     if (!cstr->IncludesModPar(mod,par)) continue;
3759     if (cstr->GetType()==AliITSAlignMille2ConstrArray::kTypeGaussian) gaussian = kTRUE;
3760     else meanmed = kTRUE;
3761     //
3762     if (meanmed && gaussian) break; // no sense to check further
3763   }
3764   //
3765   return meanmed||gaussian;
3766 }
3767
3768 //________________________________________________________________________________________________________
3769 Bool_t AliITSAlignMille2::IsParModFamilyVaried(const AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3770 {
3771   // check if parameter par is varied for this module or its children up to the level depth
3772   if (depth<0) return kFALSE;
3773   if (mod->GetParOffset(par)>=0) return kTRUE;
3774   for (int icld=mod->GetNChildren();icld--;) {
3775     AliITSAlignMille2Module* child = mod->GetChild(icld);
3776     if (IsParModFamilyVaried(child, par, depth-1)) return kTRUE;
3777   }
3778   return kFALSE;
3779   //
3780 }
3781
3782 /*
3783 //________________________________________________________________________________________________________
3784 Bool_t AliITSAlignMille2::IsParFamilyFree(AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3785 {
3786   // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
3787   if (depth<0) return kTRUE;
3788   for (int icld=mod->GetNChildren();icld--;) {
3789     AliITSAlignMille2Module* child = mod->GetChild(icld);
3790     //if (child->GetParOffset(par)<0) continue;                  // fixed
3791     Bool_t cstMM=kFALSE,cstGS=kFALSE;
3792     // does this child have gaussian constraint ?
3793     if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
3794     // check its children
3795     if (!IsParFamilyFree(child,par,depth-1)) return kTRUE;
3796   }
3797   return kFALSE;
3798   //
3799 }
3800 */
3801
3802 //________________________________________________________________________________________________________
3803 Bool_t AliITSAlignMille2::IsParFamilyFree(const AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3804 {
3805   // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
3806   if (depth<0) return kFALSE;
3807   for (int icld=mod->GetNChildren();icld--;) {
3808     AliITSAlignMille2Module* child = mod->GetChild(icld);
3809     //if (child->GetParOffset(par)<0) continue;                  // fixed
3810     Bool_t cstMM=kFALSE,cstGS=kFALSE;
3811     // does this child have gaussian constraint ?
3812     if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
3813     // check its children
3814     if (IsParFamilyFree(child,par,depth-1)) return kTRUE;
3815   }
3816   return kFALSE;
3817   //
3818 }
3819
3820 //________________________________________________________________________________________________________
3821 Double_t AliITSAlignMille2::GetTDriftSDD() const
3822 {
3823   // obtain drift time corrected for t0
3824   double t = fCluster.GetDriftTime();
3825   return t - fDriftTime0[ fCluster.GetUniqueID()-1 ];
3826 }
3827
3828 //________________________________________________________________________________________________________
3829 Double_t AliITSAlignMille2::GetVDriftSDD() const
3830 {
3831   // obtain corrected drift speed
3832   return fDriftSpeed[ fCluster.GetUniqueID()-1 ];
3833 }
3834
3835 //________________________________________________________________________________________________________
3836 Bool_t AliITSAlignMille2::FixedOrphans() const
3837 {
3838   // are there fixed modules with no parent (normally in such a case
3839   // the constraints on the orphans should not be applied
3840   if (!IsConfigured()) {
3841     AliInfo("Still not configured");
3842     return kFALSE;
3843   }
3844   for (int i=0;i<fNModules;i++) {
3845     AliITSAlignMille2Module* md = GetMilleModule(i);
3846     if (md->IsNotInConf()) continue;
3847     if (md->GetParent()==0 && md->GetNParFree()==0) return kTRUE;
3848   }
3849   return kFALSE;
3850 }
3851
3852 //________________________________________________________________________________________________________
3853 void AliITSAlignMille2::ConvertParamsToGlobal()
3854 {
3855   // convert params in local frame to global one
3856   double pars[AliITSAlignMille2Module::kMaxParGeom];
3857   for (int imd=fNModules;imd--;) {
3858     AliITSAlignMille2Module* mod = GetMilleModule(imd);
3859     if (mod->GeomParamsGlobal()) continue;
3860     mod->GetGeomParamsGlo(pars);
3861     mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
3862     mod->SetGeomParamsGlobal(kTRUE);
3863   }
3864 }
3865
3866 //________________________________________________________________________________________________________
3867 void AliITSAlignMille2::ConvertParamsToLocal()
3868 {
3869   // convert params in global frame to local one
3870   double pars[AliITSAlignMille2Module::kMaxParGeom];
3871   for (int imd=fNModules;imd--;) {
3872     AliITSAlignMille2Module* mod = GetMilleModule(imd);
3873     if (!mod->GeomParamsGlobal()) continue;
3874     mod->GetGeomParamsLoc(pars);
3875     mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
3876     mod->SetGeomParamsGlobal(kFALSE);
3877   }
3878 }
3879
3880 //________________________________________________________________________________________________________
3881 void AliITSAlignMille2::SetBField(Double_t b)
3882 {
3883   // set Bz value
3884   if (IsZero(b,1e-5)) {
3885     fBField = 0.0;
3886     fBOn = kFALSE;
3887     fNLocal = 4;
3888   }
3889   else {
3890     fBField = b;
3891     fBOn = kTRUE;
3892     fNLocal = 5; // helices
3893   }
3894 }
3895
3896 //________________________________________________________________________________________________________
3897 Int_t AliITSAlignMille2::ProcessUserInfo(TList* userInfo)
3898 {
3899   // extract calibration information used for TrackPointArray creation from run info
3900   //
3901   if (!userInfo) { AliInfo("No UserInfo is provided"); return 0;}
3902   //
3903   TMap *cdbMap=0;
3904   TList* cdbList=0;
3905   TObjString *objStr,*objStr1,*keyStr;
3906   TString cdbStr;
3907   AliCDBManager* man = AliCDBManager::Instance();
3908   man->SetCacheFlag(kFALSE);
3909   //
3910   int run = userInfo->GetUniqueID();
3911   if (run>0) SetRunID(run);
3912   AliInfo(Form("UserInfo corresponds to run#%d",run));
3913   cdbMap  = (TMap*)userInfo->FindObject("cdbMap");
3914   const TMap *curMap = man->GetStorageMap();
3915   if (!cdbMap) {AliInfo("No CDB Map found in UserInfo");}
3916   else {
3917     if ((objStr=(TObjString*)cdbMap->GetValue("default"))) { // first set default CDB path
3918       if ((objStr1=(TObjString*)curMap->GetValue("default")) && objStr1->GetUniqueID()) {
3919         AliInfo(Form("OCDB default path from UserInfo: %s is overriden by user setting %s",objStr->GetName(),objStr1->GetName()));
3920       }
3921       else {
3922         cdbStr = objStr->GetString();
3923         man->UnsetDefaultStorage();
3924         if (man->GetRaw()) man->SetRaw(kFALSE);
3925         if (cdbStr.BeginsWith("raw://")) cdbStr = "raw://";
3926         AliInfo(Form("Default CDB Storage from UserInfo: %s",cdbStr.Data()));
3927         man->SetDefaultStorage( cdbStr.Data() ); // this may be overriden later by configuration file
3928       }
3929     }
3930     if (man->GetRaw() && run>0) man->SetRun(run);
3931     //
3932     // set specific paths relevant for alignment
3933     TIter itMap(cdbMap);
3934     while( (keyStr=(TObjString*)itMap.Next()) ) {
3935       TString keyS = keyStr->GetString();
3936       if ( keyS == "default" ) continue;
3937       //
3938       TObjString* curPath = (TObjString*)curMap->GetValue(keyStr->GetName());
3939       if (curPath && curPath->GetUniqueID()) {
3940         AliInfo(Form("Storage for %s from UserInfo\n is overriden by user setting %s",keyS.Data(),curPath->GetName()));
3941         continue;
3942       }
3943       man->SetSpecificStorage( keyS.Data(), cdbMap->GetValue(keyS)->GetName() );
3944     }
3945   }
3946   //
3947   cdbList = (TList*)userInfo->FindObject("cdbList");
3948   if (!cdbList) {AliInfo("No CDB List found in UserInfo");}
3949   else {
3950     // Objects used for TrackPointArray production
3951     GetPathFromUserInfo(cdbList,"GRP/Geometry/Data",fIniGeomPath ,kSameInitGeomBit);
3952     GetPathFromUserInfo(cdbList,"ITS/Align/Data"   ,fIniDeltaPath,kSameInitDeltasBit);
3953     GetPathFromUserInfo(cdbList,"ITS/Calib/RespSDD",fIniSDDRespPath,kSameInitSDDRespBit);
3954     GetPathFromUserInfo(cdbList,"ITS/Calib/DriftSpeedSDD",fIniSDDVDriftPath,kSameInitSDDVDriftBit);
3955     GetPathFromUserInfo(cdbList,"ITS/Calib/MapsTimeSDD",fIniSDDCorrMapPath,kSameInitSDDCorrMapBit);
3956     GetPathFromUserInfo(cdbList,"GRP/Calib/MeanVertexSPD",fDiamondPath,kSameDiamondBit);
3957   }
3958   //
3959   TList *bzlst = (TList*)userInfo->FindObject("BzkGauss");
3960   if (bzlst && bzlst->At(0)) {
3961     objStr = (TObjString*)bzlst->At(0);
3962     SetBField( objStr->GetString().Atof() );
3963     AliInfo(Form("Magnetic field from UserInfo: %+.2e",GetBField()));
3964   }
3965   return 0;
3966 }
3967
3968 //________________________________________________________________________________________________________
3969 Int_t AliITSAlignMille2::GetPathFromUserInfo(TList* cdbList,const char* calib,TString& path, Int_t useBit)
3970 {
3971   // extract the path for specific CDB path from user info. If it is the same as already loaded, set corresponing bit
3972   TIter itList(cdbList);
3973   if (useBit>=0) ResetBit(useBit);
3974   TObjString* objStr;
3975   while( (objStr=(TObjString*)itList.Next()) )
3976     if (objStr->GetString().Contains(calib)) {
3977       TString newpath = objStr->GetString();
3978       AliInfo(Form("Found path in UserInfo: %s",newpath.Data()));
3979       if ( useBit>=0 && (fUserProvided&useBit) ) {
3980         AliInfo(Form("Will use the one provided in config: %s",path.Data()));
3981         SetBit(useBit);
3982       }
3983       else if ( useBit>=0 && (newpath == path) ) {
3984         AliInfo(Form("Path %s is the same as already loaded",path.Data()));
3985         SetBit(useBit);
3986       }
3987       else path = newpath;
3988       //
3989       return 0;
3990     }
3991   AliInfo(Form("Did not find path for %s in UserInfo",calib));
3992   path = "";
3993   return -1;
3994 }
3995
3996 //________________________________________________________________________________________________________
3997 Int_t AliITSAlignMille2::LoadSDDResponse(TString& path, AliITSresponseSDD *&resp)
3998 {
3999   // load SDD response
4000   if (path.IsNull()) return 0;
4001   AliInfo(Form("Loading SDD response from %s",path.Data()));
4002   //
4003   AliCDBEntry *entry = 0;
4004   delete resp;
4005   resp = 0;
4006   //
4007   while(1) {
4008     if (path.BeginsWith("path: ")) { // must load from OCDB
4009       entry = GetCDBEntry(path.Data());
4010       if (!entry) break;
4011       resp = (AliITSresponseSDD*) entry->GetObject();
4012       entry->SetObject(NULL);
4013       entry->SetOwner(kTRUE);
4014       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4015       //      delete cdbId;
4016       //      delete entry;
4017       break;
4018     }
4019     //
4020     if (gSystem->AccessPathName(path.Data())) break;
4021     TFile* precf = TFile::Open(path.Data());
4022     if (precf->FindKey("AliITSresponseSDD")) resp = (AliITSresponseSDD*)precf->Get("AliITSresponseSDD");
4023     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4024       resp = (AliITSresponseSDD*) entry->GetObject();
4025       if (resp && resp->InheritsFrom(AliITSresponseSDD::Class())) entry->SetObject(NULL);
4026       else resp = 0;
4027       entry->SetObject(NULL);
4028       entry->SetOwner(kTRUE);
4029       delete entry;
4030     }
4031     //
4032     precf->Close();
4033     delete precf;
4034     break;
4035   }
4036   //
4037   if (!resp) {AliError(Form("Failed to load SDD response from %s",path.Data())); return -1;}
4038   return 0;
4039 }
4040
4041 //________________________________________________________________________________________________________
4042 Int_t AliITSAlignMille2::LoadSDDVDrift(TString& path, TObjArray *&arr)
4043 {
4044   // load VDrift object
4045   if (path.IsNull()) return 0;
4046   AliInfo(Form("Loading SDD VDrift from %s",path.Data()));
4047   //
4048   AliCDBEntry *entry = 0;
4049   delete arr;
4050   arr = 0;
4051   while(1) {
4052     if (path.BeginsWith("path: ")) { // must load from OCDB
4053       entry = GetCDBEntry(path.Data());
4054       if (!entry) break;
4055       arr = (TObjArray*) entry->GetObject();
4056       entry->SetObject(NULL);
4057       entry->SetOwner(kTRUE);
4058       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4059       //      delete cdbId;
4060       //      delete entry;
4061       break;
4062     }
4063     //
4064     if (gSystem->AccessPathName(path.Data())) break;
4065     TFile* precf = TFile::Open(path.Data());
4066     if (precf->FindKey("TObjArray")) arr = (TObjArray*)precf->Get("TObjArray");
4067     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4068       arr = (TObjArray*) entry->GetObject();
4069       if (arr && arr->InheritsFrom(TObjArray::Class())) entry->SetObject(NULL);
4070       else arr = 0;
4071       entry->SetObject(NULL);
4072       entry->SetOwner(kTRUE);
4073       delete entry;
4074     }
4075     //
4076     precf->Close();
4077     delete precf;
4078     break;
4079   }
4080   //
4081   if (!arr) {AliError(Form("Failed to load SDD vdrift from %s",path.Data())); return -1;}
4082   arr->SetOwner(kTRUE);
4083   return 0;
4084 }
4085
4086 //________________________________________________________________________________________________________
4087 Int_t AliITSAlignMille2::LoadSDDCorrMap(TString& path, AliITSCorrectSDDPoints *&map)
4088 {
4089   // Load SDD correction map
4090   //
4091   if (path.IsNull()) return 0;
4092   AliInfo(Form("Loading SDD Correction Maps from %s",path.Data()));
4093   //
4094   AliCDBEntry *entry = 0;
4095   delete map;
4096   map = 0;
4097   TObjArray* arr = 0;
4098   while(1) {
4099     if (path.BeginsWith("path: ")) { // must load from OCDB
4100       entry = GetCDBEntry(path.Data());
4101       if (!entry) break;
4102       arr = (TObjArray*) entry->GetObject();
4103       entry->SetObject(NULL);
4104       entry->SetOwner(kTRUE);
4105       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4106       //      delete cdbId;
4107       //      delete entry;
4108       break;
4109     }
4110     //
4111     if (gSystem->AccessPathName(path.Data())) break;
4112     TFile* precf = TFile::Open(path.Data());
4113     if (precf->FindKey("TObjArray")) arr = (TObjArray*)precf->Get("TObjArray");
4114     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4115       arr = (TObjArray*) entry->GetObject();
4116       if (arr && arr->InheritsFrom(TObjArray::Class())) entry->SetObject(NULL);
4117       else arr = 0;
4118       entry->SetObject(NULL);
4119       entry->SetOwner(kTRUE);
4120       delete entry;
4121     }
4122     //
4123     precf->Close();
4124     delete precf;
4125     break;
4126   }
4127   //
4128   if (!arr) {AliError(Form("Failed to load SDD Correction Map from %s",path.Data())); return -1;}
4129   arr->SetOwner(kTRUE);
4130   map = new AliITSCorrectSDDPoints(arr);
4131
4132   return 0;
4133 }
4134
4135 //________________________________________________________________________________________________________
4136 Int_t AliITSAlignMille2::LoadDiamond(TString& path)
4137 {
4138   // load vertex constraint
4139   if (path.IsNull()) return 0;
4140   AliInfo(Form("Loading Diamond Constraint from %s",path.Data()));
4141   //
4142   AliCDBEntry *entry = 0;
4143   AliESDVertex *vtx = 0;
4144   while(1) {
4145     if (path.BeginsWith("path: ")) { // must load from OCDB
4146       entry = GetCDBEntry(path.Data());
4147       if (!entry) break;
4148       vtx = (AliESDVertex*) entry->GetObject();
4149       entry->SetObject(NULL);
4150       entry->SetOwner(kTRUE);
4151       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4152       //      delete cdbId;
4153       //      delete entry;
4154       break;
4155     }
4156     //
4157     if (gSystem->AccessPathName(path.Data())) break;
4158     TFile* precf = TFile::Open(path.Data());
4159     if (precf->FindKey("AliESDVertex")) vtx = (AliESDVertex*)precf->Get("AliESDVertex");
4160     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4161       vtx = (AliESDVertex*) entry->GetObject();
4162       if (vtx && vtx->InheritsFrom(AliESDVertex::Class())) entry->SetObject(NULL);
4163       else vtx = 0;
4164       entry->SetObject(NULL);
4165       entry->SetOwner(kTRUE);
4166       delete entry;
4167     }
4168     //
4169     precf->Close();
4170     delete precf;
4171     break;
4172   }
4173   //
4174   if (!vtx) {AliError(Form("Failed to load Diamond constraint from %s",path.Data())); return -1;}
4175   //
4176   double vtxXYZ[3];
4177   vtx->GetXYZ(vtxXYZ);
4178   for (int i=3;i--;) vtxXYZ[i] -= fCorrDiamond[i];
4179   vtx->SetXYZ(vtxXYZ);
4180   SetVertexConstraint(vtx);
4181   AliInfo("Will use following Diamond Constraint (errors inverted):");
4182   fDiamondI.Print("");
4183   delete vtx;
4184   return 0;
4185 }
4186
4187 //________________________________________________________________________________________________________
4188 Int_t AliITSAlignMille2::LoadDeltas(TString& path, TClonesArray *&arr)
4189 {
4190   // load ITS geom deltas
4191   if (path.IsNull()) return 0;
4192   AliInfo(Form("Loading Alignment Deltas from %s",path.Data()));
4193   //
4194   AliCDBEntry *entry = 0;
4195   delete arr;
4196   arr = 0;
4197   while(1) {
4198     if (path.BeginsWith("path: ")) { // must load from OCDB
4199       entry = GetCDBEntry(path.Data());
4200       if (!entry) break;
4201       arr = (TClonesArray*) entry->GetObject();
4202       entry->SetObject(NULL);
4203       entry->SetOwner(kTRUE);
4204       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4205       //      delete cdbId;
4206       //      delete entry;
4207       break;
4208     }
4209     //
4210     if (gSystem->AccessPathName(path.Data())) break;
4211     TFile* precf = TFile::Open(path.Data());
4212     if (precf->FindKey("ITSAlignObjs")) arr = (TClonesArray*)precf->Get("ITSAlignObjs");
4213     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4214       arr = (TClonesArray*) entry->GetObject();
4215       if (arr && arr->InheritsFrom(TClonesArray::Class())) entry->SetObject(NULL);
4216       else arr = 0;
4217       entry->SetObject(NULL);
4218       entry->SetOwner(kTRUE);
4219       delete entry;
4220     }
4221     precf->Close();
4222     delete precf;
4223     break;
4224   }
4225   //
4226   if (!arr) {AliError(Form("Failed to load Deltas from %s",path.Data())); return -1;}
4227   //
4228   return 0;
4229 }
4230
4231 //________________________________________________________________________________________________________
4232 Int_t AliITSAlignMille2::CacheMatricesCurr()
4233 {
4234   // build arrays for the fast access to sensor matrices from their sensor ID
4235   //
4236   TGeoHMatrix mdel;
4237   AliInfo("Building sensors current matrices cache");
4238   //
4239   fCacheMatrixCurr.Delete();
4240   for (int idx=0;idx<=kMaxITSSensID;idx++) {
4241     int volID = AliITSAlignMille2Module::GetVolumeIDFromIndex(idx);
4242     TGeoHMatrix *mcurr = new TGeoHMatrix();
4243     AliITSAlignMille2Module::SensVolMatrix(volID, mcurr);
4244     fCacheMatrixCurr.AddAtAndExpand(mcurr,idx);
4245     //
4246   }
4247   //
4248   TGeoHMatrix *mcurr = new TGeoHMatrix();
4249   fCacheMatrixCurr.AddAtAndExpand(mcurr,kVtxSensID); // special unit matrix for diamond constraint
4250   //
4251   fCacheMatrixCurr.SetOwner(kTRUE);
4252   return 0;
4253 }
4254
4255 //________________________________________________________________________________________________________
4256 Int_t AliITSAlignMille2::CacheMatricesOrig()
4257 {
4258   // build arrays for the fast access to sensor original matrices (used for production)
4259   //
4260   TGeoHMatrix mdel;
4261   AliInfo("Building sensors original matrices cache");
4262   //
4263   /*if (fIniGeomPath!=fGeometryPath)*/ if (LoadGeometry(fIniGeomPath)) {AliInfo("Failed to re-load ideal geometry");exit(1);}
4264   //
4265   fCacheMatrixOrig.Delete();
4266   if (!fIniDeltaPath.IsNull()) {
4267     TClonesArray* prealSav = fPrealignment;
4268     fPrealignment = 0;
4269     if (LoadDeltas(fIniDeltaPath,fPrealignment) || ApplyToGeometry())
4270       { AliError("Failed to load/apply initial deltas used to produce points"); return -1;}
4271     delete fPrealignment;
4272     fPrealignment = prealSav;
4273   }
4274   //
4275   for (int idx=0;idx<=kMaxITSSensID;idx++) {
4276     int volID = AliITSAlignMille2Module::GetVolumeIDFromIndex(idx);
4277     TGeoHMatrix *morig = new TGeoHMatrix();
4278     AliITSAlignMille2Module::SensVolMatrix(volID,morig);
4279     fCacheMatrixOrig.AddAtAndExpand(morig,idx);
4280   }
4281   //
4282   if (fConvertPreDeltas) {
4283     // in order to convert deltas from old to new geometry we need the final matrices for all alignable objects
4284     int nmat = fGeoManager->GetNAlignable();
4285     fConvAlgMatOld.Delete();
4286     int nmatSel = 0;
4287     for (int i=0;i<nmat;i++) {
4288       TString nm = fGeoManager->GetAlignableEntry(i)->GetName();
4289       if (!nm.BeginsWith("ITS")) continue;
4290       TGeoHMatrix *mo = new TGeoHMatrix();
4291       (*mo) = *(AliGeomManager::GetMatrix(nm));
4292       fConvAlgMatOld.AddAtAndExpand(mo,nmatSel++);
4293       mo->SetTitle(nm);
4294       mo->SetName(nm);
4295     }
4296     ConvSortHierarchically(fConvAlgMatOld);
4297   }
4298   //
4299   TGeoHMatrix *mcurr = new TGeoHMatrix();
4300   fCacheMatrixOrig.AddAtAndExpand(mcurr,kVtxSensID); // special unit matrix for diamond constraint
4301   //
4302   fCacheMatrixOrig.SetOwner(kTRUE);
4303
4304   fUsePreAlignment = 0;
4305   LoadGeometry(fGeometryPath);   // reload target geometry
4306   //
4307   return 0;
4308 }
4309
4310 //________________________________________________________________________________________________________
4311 void AliITSAlignMille2::RemoveHelixFitConstraint()
4312 {
4313   // suppress constraint
4314   fConstrCharge = 0;
4315   fConstrPT = fConstrPTErr = -1;
4316 }
4317
4318 //________________________________________________________________________________________________________
4319 void AliITSAlignMille2::ConstrainHelixFitPT(Int_t q,Double_t pt,Double_t pterr)
4320 {
4321   // constrain q and pT of the helical fit of the track (should be set before process.track)
4322   //
4323   fConstrCharge = q==0 ? q:TMath::Sign(1,q);
4324   fConstrPT = pt;
4325   fConstrPTErr = pterr;
4326 }
4327
4328 //________________________________________________________________________________________________________
4329 void AliITSAlignMille2::ConstrainHelixFitCurv(Int_t q,Double_t crv,Double_t crverr)
4330 {
4331   // constrain charge and curvature of the helical fit of the track (should be set before process.track)
4332   //
4333   const double kCQConv = 0.299792458e-3;// R = PT/Bz/fgkCQConv with GeV,kGauss,cm
4334
4335   fConstrCharge = q==0 ? q:TMath::Sign(1,q);
4336   if (crv<0 || IsZero(crv)) {
4337     fConstrPT    = -1;
4338     fConstrPTErr = -1;
4339   }
4340   else {
4341     fConstrPT    = TMath::Abs(1./crv*fBField*kCQConv);
4342     fConstrPTErr = crverr>1e-10 ? TMath::Abs(fConstrPT/crv*crverr) : 0.;
4343   }
4344 }
4345
4346 //________________________________________________________________________________________________________
4347 TClonesArray* AliITSAlignMille2::CreateDeltas()
4348 {
4349   // Create \Deltas for every explicitly or implicitly (via non-alignable volumes) varied
4350   // or prealigned module.
4351   // If the module has inded J in the hierarchy of alignable volumes (0 - the top, most
4352   // coarse level), then its Delta is expressed via MP2 \deltas (in global frame) and
4353   // prealignment \DeltaP's as:
4354   // \Delta_J = Y X Y^-1
4355   // where X = \delta_J * \DeltaP_J
4356   // Y = Prod_{K=0,J-1} \delta_K
4357   // Note that \delta_L accounts not only for its own correction but also of all non-alignable
4358   // modules in the hierarchy chain from L up to the closest alignable:
4359   // while (parent && !parent->IsAlignable()) {
4360   //   \delta_L->MultiplyLeft( \delta_parent );
4361   //   parent = parent->GetParent();
4362   // }
4363   //
4364   Bool_t convLoc = kFALSE;
4365   if (!GetUseGlobalDelta()) {
4366     ConvertParamsToGlobal();
4367     convLoc = kTRUE;
4368   }
4369   //
4370   AliAlignObjParams tempAlignObj;
4371   TGeoHMatrix tempMatX,tempMatY,tempMat1;
4372   //
4373   TClonesArray *array = new TClonesArray("AliAlignObjParams",10);
4374   TClonesArray &alobj = *array;
4375   int idx = 0;
4376   //
4377   TGeoManager* geoManager = AliGeomManager::GetGeometry();
4378   int nalgtot = geoManager->GetNAlignable();
4379   //
4380   for (int ialg=0;ialg<nalgtot;ialg++) {             // loop over all alignable entries
4381     //
4382     const char* algname = geoManager->GetAlignableEntry(ialg)->GetName();
4383     //
4384     AliITSAlignMille2Module* md     = GetMilleModuleBySymName(algname); // explicitly varied?
4385     AliITSAlignMille2Module* parent = md ? md->GetParent(): GetMilleModuleIfContained(algname);
4386     if (md && parent) {
4387       TString mdName = md->GetName();
4388       TString prName = parent->GetName();
4389       // SPD Sector -> Layer parentship is fake, need special treatment
4390       if ( mdName.CountChar('/')==2 && mdName.BeginsWith("ITS/SPD") && // SPD sector
4391            prName.CountChar('/')==1 && mdName.BeginsWith("ITS/SPD") )  // SPD Layer
4392         parent = parent ? parent->GetParent(): GetMilleModuleIfContained(prName.Data());
4393     }
4394     //
4395     AliAlignObjParams*       preob  = GetPrealignedObject(algname);  // was it prealigned ?
4396     //
4397     if (!preob && !md && (!parent || parent->IsAlignable())) continue; // noting to do
4398     //
4399     // create matrix X (see comment) ------------------------------------------------->>>
4400     // start from unity matrix
4401     tempMatX.Clear();
4402     if (preob) {   // account prealigngment
4403       preob->GetMatrix(tempMat1);
4404       tempMatX.MultiplyLeft(&tempMat1);
4405     }
4406     //
4407     if (md) {
4408       tempAlignObj.SetTranslation( md->GetParVal(0),md->GetParVal(1),md->GetParVal(2));
4409       tempAlignObj.SetRotation(    md->GetParVal(3),md->GetParVal(4),md->GetParVal(5));
4410       tempAlignObj.GetMatrix(tempMat1);
4411       tempMatX.MultiplyLeft(&tempMat1);  // acount correction to varied module
4412     }
4413     //
4414     // the corrections to all non-alignable modules from current on
4415     // till first alignable should add up to its matrix
4416     while (parent && !parent->IsAlignable()) {
4417       tempAlignObj.SetTranslation( parent->GetParVal(0),parent->GetParVal(1),parent->GetParVal(2));
4418       tempAlignObj.SetRotation(    parent->GetParVal(3),parent->GetParVal(4),parent->GetParVal(5));
4419       tempAlignObj.GetMatrix(tempMat1);
4420       tempMatX.MultiplyLeft(&tempMat1);  // add matrix of non-alignable module
4421       parent = parent->GetParent();
4422     }
4423     // create matrix X (see comment) ------------------------------------------------<<<
4424     //
4425     // create matrix Y (see comment) ------------------------------------------------>>>
4426     // start from unity matrix
4427     tempMatY.Clear();
4428     while ( parent ) {
4429       tempAlignObj.SetTranslation( parent->GetParVal(0),parent->GetParVal(1),parent->GetParVal(2));
4430       tempAlignObj.SetRotation(    parent->GetParVal(3),parent->GetParVal(4),parent->GetParVal(5));
4431       tempAlignObj.GetMatrix(tempMat1);
4432       tempMatY.MultiplyLeft(&tempMat1);
4433       parent = parent->GetParent();
4434     }
4435     // create matrix Y (see comment) ------------------------------------------------<<<
4436     //
4437     tempMatX.MultiplyLeft(&tempMatY);
4438     tempMatX.Multiply(&tempMatY.Inverse());
4439     //
4440     if (tempMatX.IsIdentity()) continue; // do not store dummy matrices
4441     UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(algname);
4442     new(alobj[idx++]) AliAlignObjParams(algname,vid,tempMatX,kTRUE);
4443     //
4444   }
4445   //
4446   if (convLoc) ConvertParamsToLocal();
4447   //
4448   return array;
4449   //
4450 }
4451
4452 //_______________________________________________________________________________________
4453 AliITSresponseSDD* AliITSAlignMille2::CreateSDDResponse()
4454 {
4455   // create object with SDD repsonse (t0 and vdrift corrections) accounting for
4456   // eventual precalibration
4457   //
4458   // if there was a precalibration provided, copy it to new arrray
4459   AliITSresponseSDD *precal = GetSDDPrecalResp();
4460   if (!precal && fIniVDriftSDD)       precal = GetSDDInitResp();    // InitResp is used only when IniVDrift is provided
4461   Bool_t isPreCalMult = precal&&precal->IsVDCorrMult() ? kTRUE : kFALSE;
4462   AliITSresponseSDD *calibSDD = new AliITSresponseSDD();
4463   calibSDD->SetVDCorrMult(fIsSDDVDriftMult);
4464   //
4465   // copy initial values to the new object
4466   if (precal) {
4467     calibSDD->SetTimeOffset(precal->GetTimeOffset());
4468     calibSDD->SetADC2keV(precal->GetADC2keV());
4469     calibSDD->SetChargevsTime(precal->GetChargevsTime());
4470     for (int ind=kSDDoffsID;ind<kSDDoffsID+kNSDDmod;ind++) {
4471       calibSDD->SetModuleTimeZero(ind, precal->GetTimeZero(ind));
4472       calibSDD->SetDeltaVDrift(ind, precal->GetDeltaVDrift(ind,kFALSE),kFALSE); // left
4473       calibSDD->SetDeltaVDrift(ind, precal->GetDeltaVDrift(ind,kTRUE ),kTRUE);  // right
4474       calibSDD->SetADCtokeV(ind,precal->GetADCtokeV(ind));
4475     }
4476   }
4477   else for (int ind=kSDDoffsID;ind<kSDDoffsID+kNSDDmod;ind++) calibSDD->SetModuleTimeZero(ind,0);
4478   //
4479   Bool_t save = kFALSE;
4480   for (int imd=GetNModules();imd--;) {
4481     AliITSAlignMille2Module* md = GetMilleModule(imd);
4482     if (!md->IsSDD()) continue;
4483     if (md->IsFreeDOF(AliITSAlignMille2Module::kDOFT0)  ||
4484         md->IsFreeDOF(AliITSAlignMille2Module::kDOFDVL) ||
4485         md->IsFreeDOF(AliITSAlignMille2Module::kDOFDVR)) save = kTRUE;
4486         //
4487     for (int is=0;is<md->GetNSensitiveVolumes();is++) {
4488       int ind  = md->GetSensVolIndex(is);
4489       float t0  = calibSDD->GetTimeZero(ind)      + md->GetParVal(AliITSAlignMille2Module::kDOFT0);
4490       double dvL = md->GetParVal(AliITSAlignMille2Module::kDOFDVL);
4491       double dvR = md->GetParVal(AliITSAlignMille2Module::kDOFDVR);
4492       if (!calibSDD->IsVDCorrMult()) { // save as additive correction
4493         dvL *= 1e4;
4494         dvR *= 1e4;
4495         //
4496         double conv = 1;
4497         if (isPreCalMult) conv = 6.4; // convert multiplicative precal correction to additive
4498         dvL += calibSDD->GetDeltaVDrift(ind,kFALSE)*conv;
4499         dvR += calibSDD->GetDeltaVDrift(ind,kTRUE)*conv;
4500       }
4501       else { // save as multipicative correction
4502         double conv = 1;
4503         if (!isPreCalMult) conv = 1./6.4; // convert additive precal correction to multiplicative
4504         dvL += calibSDD->GetDeltaVDrift(ind,kFALSE)*conv;
4505         dvR += calibSDD->GetDeltaVDrift(ind,kTRUE)*conv;
4506       }
4507       //
4508       calibSDD->SetModuleTimeZero(ind, t0);
4509       calibSDD->SetDeltaVDrift(ind, dvL, kFALSE); // left  side correction
4510       calibSDD->SetDeltaVDrift(ind, dvR, kTRUE); // right side correction
4511     }
4512   }
4513   //
4514   if (!save) {
4515     AliInfo("No free parameters for SDD calibration, nothing to save");
4516     delete calibSDD;
4517     calibSDD = 0;
4518   }
4519   //
4520   return calibSDD;
4521 }
4522
4523 //_______________________________________________________________________________________
4524 Int_t AliITSAlignMille2::ReloadInitCalib(TList *userInfo)
4525 {
4526   // Use provided UserInfo to
4527   // load the initial calib parameters (geometry, SDD response...)
4528   // Can be used if set of data was processed with different calibration
4529   //
4530   if (!userInfo) {
4531     AliInfo("Reloading of the Calibration parameters was called with empty userInfo");
4532     return 1;
4533   }
4534   if (ProcessUserInfo(userInfo)) {
4535     AliInfo("Error in processing user info");
4536     userInfo->Print();
4537     exit(1);
4538   }
4539   return ReloadInitCalib();
4540 }
4541
4542 //_______________________________________________________________________________________
4543 Int_t AliITSAlignMille2::ReloadInitCalib()
4544 {
4545   // Load the initial calib parameters (geometry, SDD response...)
4546   // Can be used if set of data was processed with different calibration
4547   //
4548   // 1st cache original matrices
4549   if (!(TestBit(kSameInitDeltasBit) && TestBit(kSameInitGeomBit))) { // need to reload geometry
4550     //
4551     if (CacheMatricesOrig()) {
4552       AliInfo("Failed to cache new initial geometry");
4553       exit(1);
4554     }
4555     // RS : commented because we don't need to reload prealignment deltas, they are already loaded
4556     // then reload the prealignment geometry
4557     //    if (LoadDeltas(fPreDeltaPath,fPrealignment)) {
4558     //      AliInfo(Form("Failed to reload the prealigned geometry %s",fPreDeltaPath.Data()));
4559     //      exit(1);
4560     //    }
4561     //
4562     if (fPrealignment && ApplyToGeometry()) {
4563       AliInfo(Form("Failed re-apply prealigned geometry %s",fPreDeltaPath.Data()));
4564       exit(1);
4565     }
4566     //
4567     // usually no need to re-cache the prealignment geometry, it was not changed
4568     if (fCacheMatrixCurr.GetEntriesFast() != fCacheMatrixOrig.GetEntriesFast()) {
4569       //      CacheMatricesCurr();
4570       AliInfo(Form("Failed to cache the prealigned geometry %s",fPreDeltaPath.Data()));
4571       exit(1);
4572     }
4573   }
4574   else ResetBit(kSameInitDeltasBit);
4575   //
4576   // reload initial SDD response
4577   if (!TestBit(kSameInitSDDRespBit)) {
4578     if (LoadSDDResponse(fIniSDDRespPath, fIniRespSDD) ) {
4579       AliInfo(Form("Failed to load new SDD response %s",fIniSDDRespPath.Data()));
4580       exit(1);
4581     }
4582   }
4583   else ResetBit(kSameInitSDDRespBit);
4584   //
4585   // reload initial SDD vdrift
4586   if (!TestBit(kSameInitSDDVDriftBit)) {
4587     if (LoadSDDVDrift(fIniSDDVDriftPath, fIniVDriftSDD) ) {
4588       AliInfo(Form("Failed to load new SDD VDrift %s",fIniSDDVDriftPath.Data()));
4589       exit(1);
4590     }
4591   }
4592   else ResetBit(kSameInitSDDRespBit);
4593   //
4594   // reload SDD corr.map
4595   if (!TestBit(kSameInitSDDCorrMapBit)) {
4596     if (LoadSDDCorrMap(fIniSDDCorrMapPath, fIniCorrMapSDD) ) {
4597       AliInfo(Form("Failed to load new SDD Correction Map %s",fIniSDDCorrMapPath.Data()));
4598       exit(1);
4599     }
4600   }
4601   else ResetBit(kSameInitSDDRespBit);
4602   //
4603   // reload diamond info
4604   if (!TestBit(kSameDiamondBit)) {
4605     if (LoadDiamond(fDiamondPath) ) {
4606       AliInfo(Form("Failed to load new Diamond constraint %s",fDiamondPath.Data()));
4607       exit(1);
4608     }
4609   }
4610   else ResetBit(kSameInitSDDRespBit);
4611   //
4612   return 0;
4613 }
4614
4615 //_______________________________________________________________________________________
4616 void AliITSAlignMille2::JacobianPosGloLoc(int locid,double* jacobian)
4617 {
4618   // calculate the locid row of the jacobian for transformation of the local coordinate to global at current point
4619   TGeoHMatrix* mat = GetSensorCurrMatrixSID(fCurrentSensID);
4620   const Double_t dpar = 1e-2;
4621   double sav = fMeasLoc[locid];
4622   fMeasLoc[locid] += dpar;
4623   mat->LocalToMaster(fMeasLoc,jacobian);
4624   fMeasLoc[locid] = sav; // recover original value
4625   for (int i=3;i--;) jacobian[i] = (jacobian[i]-fMeasGlo[i])/dpar; // the transformation is linear!!!
4626 }
4627
4628 //_______________________________________________________________________________________
4629 void AliITSAlignMille2::TieSDDVDriftsLR(AliITSAlignMille2Module* mod)
4630 {
4631   // impose equality of Left/Right sides VDrift correction for SDD
4632   ResetLocalEquation();
4633   if ( (mod->IsFreeDOF(AliITSAlignMille2Module::kDOFDVL) + mod->IsFreeDOF(AliITSAlignMille2Module::kDOFDVR))==1) {
4634     AliError("Left/Right VDrift equality is requested for SDD module with only one side VDrift free");
4635     mod->Print();
4636     return;
4637   }
4638   if (mod->GetParOffset(AliITSAlignMille2Module::kDOFDVL)>=0) SetGlobalDerivative(mod->GetParOffset(AliITSAlignMille2Module::kDOFDVL),  1.);
4639   if (mod->GetParOffset(AliITSAlignMille2Module::kDOFDVR)>=0) SetGlobalDerivative(mod->GetParOffset(AliITSAlignMille2Module::kDOFDVR), -1.);
4640   AddConstraint(fGlobalDerivatives, 0, 1e-12);
4641   //
4642 }
4643
4644 //_______________________________________________________________________________________
4645 void AliITSAlignMille2::ProcessSDDPointInfo(const AliTrackPoint* pnt,Int_t sID, Int_t pntID)
4646 {
4647   // extract the drift information from SDD track point
4648   //
4649   fDriftTime0[pntID] = fIniRespSDD ? fIniRespSDD->GetTimeZero(sID) : 0.;
4650   double tdif = pnt->GetDriftTime() - fDriftTime0[pntID];
4651   if (tdif<0) tdif = 1;
4652   //
4653   // VDrift extraction
4654   double vdrift=0,vdrift0=0;
4655   Bool_t sddSide = kFALSE;
4656   int sID0 = 2*(sID-kSDDoffsID);
4657   double zanode = -999;
4658   //
4659   if (fIniVDriftSDD) { // SDD VDrift object is provided, use the vdrift from it
4660     AliITSDriftSpeedArraySDD* drarr;
4661     double vdR,vdL,xlR,xlL;
4662     // sometimes xlocal on right side is negative due to the wrong calibration, need to test both hypothesis
4663     double xlabs = TMath::Abs(fMeasLoc[kX]);
4664     drarr  = (AliITSDriftSpeedArraySDD*)fIniVDriftSDD->At(sID0); // left side, xloc>0
4665     zanode = fSegmentationSDD->GetAnodeFromLocal(xlabs,fMeasLoc[kZ]);
4666     vdL    = drarr->GetDriftSpeed(0, zanode);
4667     if (fIniRespSDD) {
4668       double corr = fIniRespSDD->GetDeltaVDrift(sID, kFALSE);
4669       if (fIniRespSDD->IsVDCorrMult()) vdL *= (1+corr);
4670       else vdL += corr;
4671     }
4672     xlL    = (fSegmentationSDD->Dx() - vdL*tdif)*1e-4;
4673     //
4674     drarr  = (AliITSDriftSpeedArraySDD*)fIniVDriftSDD->At(sID0+1); // right side, xloc<0
4675     zanode = fSegmentationSDD->GetAnodeFromLocal(-xlabs,fMeasLoc[kZ]) - 256;
4676     vdR    = drarr->GetDriftSpeed(0, zanode);
4677     if (fIniRespSDD) {
4678       double corr = fIniRespSDD->GetDeltaVDrift(sID, kTRUE);
4679       if (fIniRespSDD->IsVDCorrMult()) vdR *= (1+corr);
4680       else vdR += corr;
4681     }
4682     xlR    = -(fSegmentationSDD->Dx() - vdR*tdif)*1e-4;
4683     //
4684     if (TMath::Abs(xlL-fMeasLoc[kX])<TMath::Abs(xlR-fMeasLoc[kX])) {
4685       sddSide = 0; // left side
4686       vdrift  = vdL*1e-4;
4687     }
4688     else {         // right side
4689       sddSide = 1;
4690       vdrift  = vdR*1e-4;
4691     }
4692     //
4693   }
4694   else { // try to determine the vdrift from the xloc
4695     vdrift = (fSegmentationSDD->Dx()*1e-4 - TMath::Abs(fMeasLoc[kX]))/tdif;
4696     sddSide = fMeasLoc[kX]<0; // 0 = left (xloc>0) ; 1 = right (xloc<1)
4697   }
4698   //
4699   if (fPreVDriftSDD) { // use imposed vdrift as a starting point
4700     zanode = fSegmentationSDD->GetAnodeFromLocal(0.5-sddSide,fMeasLoc[kZ]);
4701     if (sddSide) zanode -= 256;
4702     vdrift = ((AliITSDriftSpeedArraySDD*)fPreVDriftSDD->At(sID0+sddSide))->GetDriftSpeed(0, zanode)*1e-4;
4703   }
4704   //
4705   if (vdrift<0) vdrift = 0;
4706   vdrift0 = vdrift;
4707   // at this point we have vdrift and t0 used to create the original point.
4708   // see if precalibration was provided
4709   if (fPreRespSDD) {
4710     float t0Upd = fPreRespSDD->GetTimeZero(sID);
4711     double corr = fPreRespSDD->GetDeltaVDrift(sID, sddSide);
4712     if (fPreRespSDD->IsVDCorrMult()) vdrift *= 1+corr; // right side (xloc<0) may have different correction
4713     else                             vdrift += corr*1e-4;
4714     //
4715     // if IniRespSDD was used, it should be subtracted back, since it is accounted in the PreResp
4716     if (fIniVDriftSDD&&fIniRespSDD) {
4717       double corr1 = fIniRespSDD->GetDeltaVDrift(sID, sddSide);
4718       if (fIniRespSDD->IsVDCorrMult()) vdrift *= (1-corr1);
4719       else vdrift -= corr1*1e-4;
4720     }
4721     tdif    = pnt->GetDriftTime() - t0Upd;
4722     // correct Xlocal
4723     fMeasLoc[0] = fSegmentationSDD->Dx()*1e-4 - vdrift*tdif;
4724     if (sddSide) fMeasLoc[0] = -fMeasLoc[0];
4725     fDriftTime0[pntID] =  t0Upd;
4726   }
4727   //
4728   if (fPreCorrMapSDD) { // apply correction map
4729     fMeasLoc[0] += fPreCorrMapSDD->GetCorrection(sID,fMeasLoc[2],fMeasLoc[0]);
4730   }
4731
4732   // TEMPORARY CORRECTION (if provided) --------------<<<
4733   fDriftSpeed[pntID]  = sddSide ? -vdrift  : vdrift;
4734   fDriftSpeed0[pntID] = sddSide ? -vdrift0 : vdrift0;
4735   //
4736   //  printf("#%d: t:%+e x:%+e v:%+e: side:%d\n",pntID,fDriftTime0[pntID],fMeasLoc[0],fDriftSpeed[pntID],sddSide);
4737 }
4738
4739 //_______________________________________________________________________________________
4740 AliITSAlignMille2Module* AliITSAlignMille2::CreateVertexModule()
4741 {
4742   // creates dummy module for vertex constraint
4743   TGeoHMatrix mt;
4744   AliITSAlignMille2Module* mod = new AliITSAlignMille2Module(kVtxSensID,kVtxSensVID,"VTX",&mt,0,0);
4745   fMilleModule.AddAtAndExpand(mod,fNModules);
4746   mod->SetGeomParamsGlobal(fUseGlobalDelta);
4747   fDiamondModID = fNModules;
4748   mod->SetUniqueID(fNModules++);
4749   mod->SetNotInConf(kTRUE);
4750   return mod;
4751   //
4752 }
4753
4754 //_______________________________________________________________________________________
4755 AliCDBEntry* AliITSAlignMille2::GetCDBEntry(const char* path)
4756 {
4757   // return object from the OCDB
4758   AliCDBEntry *entry = 0;
4759   AliInfo(Form("Loading object %s",path));
4760   AliCDBManager* man = AliCDBManager::Instance();
4761   AliCDBId* cdbId = AliCDBId::MakeFromString(path);
4762   if (!cdbId) {
4763     AliError("Failed to create cdbId");
4764     return 0;
4765   }
4766   //
4767   AliCDBStorage* stor = man->GetDefaultStorage();
4768   if (!stor && !man->GetRaw()) man->SetDefaultStorage("raw://");
4769   if (man->GetRaw()) man->SetRun(fRunID>0 ? fRunID : cdbId->GetFirstRun());
4770   if (stor) {
4771     TString tp = stor->GetType();
4772     if (tp.Contains("alien",TString::kIgnoreCase) && !gGrid) TGrid::Connect("alien:");
4773   }
4774   entry = man->Get(cdbId->GetPath(),cdbId->GetFirstRun(),cdbId->GetVersion(),cdbId->GetSubVersion());
4775   //  entry = man->Get( *cdbId );
4776   man->ClearCache();
4777   //
4778   delete cdbId;
4779   return entry;
4780   //
4781 }
4782
4783 //_______________________________________________________________________________________
4784 void AliITSAlignMille2::SetVertexConstraint(const AliESDVertex* vtx)
4785 {
4786   // set vertex for constraint
4787   if (!vtx) return;
4788   //
4789   double cmat[6];
4790   float cmatF[6];
4791   vtx->GetCovMatrix(cmat);
4792   AliITSAlignMille2Module* diamMod = GetMilleModuleByVID(kVtxSensVID);
4793   if (diamMod) {
4794     cmat[0] *= diamMod->GetSigmaXFactor()*diamMod->GetSigmaXFactor();
4795     cmat[2] *= diamMod->GetSigmaYFactor()*diamMod->GetSigmaYFactor();
4796     cmat[5] *= diamMod->GetSigmaZFactor()*diamMod->GetSigmaZFactor();
4797     cmat[1] *= diamMod->GetSigmaXFactor()*diamMod->GetSigmaYFactor();
4798     cmat[3] *= diamMod->GetSigmaXFactor()*diamMod->GetSigmaZFactor();
4799     cmat[4] *= diamMod->GetSigmaYFactor()*diamMod->GetSigmaZFactor();
4800   }
4801   cmatF[0] = cmat[0]; // xx
4802   cmatF[1] = cmat[1]; // xy
4803   cmatF[2] = cmat[3]; // xz
4804   cmatF[3] = cmat[2]; // yy
4805   cmatF[4] = cmat[4]; // yz
4806   cmatF[5] = cmat[5]; // zz
4807
4808   fDiamond.SetXYZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), cmatF);
4809   //
4810   Double_t t0 = cmat[2]*cmat[5] - cmat[4]*cmat[4];
4811   Double_t t1 = cmat[1]*cmat[5] - cmat[3]*cmat[4];
4812   Double_t t2 = cmat[1]*cmat[4] - cmat[2]*cmat[3];
4813   Double_t det = cmat[0]*t0 - cmat[1]*t1 + cmat[3]*t2;
4814   if (TMath::Abs(det)<1e-36) {
4815     vtx->Print();
4816     AliFatal("Vertex constraint cov.matrix is singular");
4817   }
4818   cmatF[0] =  t0/det;
4819   cmatF[1] = -t1/det;
4820   cmatF[2] =  t2/det;
4821   cmatF[3] =  (cmat[0]*cmat[5] - cmat[3]*cmat[3])/det;
4822   cmatF[4] =  (cmat[1]*cmat[3] - cmat[0]*cmat[4])/det;
4823   cmatF[5] =  (cmat[0]*cmat[2] - cmat[1]*cmat[1])/det;
4824   fDiamondI.SetXYZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), cmatF);
4825   fVertexSet = kTRUE;
4826   //
4827 }
4828
4829 //_______________________________________________________________________________________
4830 void AliITSAlignMille2::ConvertDeltas()
4831 {
4832   // convert prealignment deltas from old geometry to new one
4833   // NOTE: the target geometry must be loaded at time this method is called
4834   //
4835   // NOTE: This method can be ONLY used when as a prealignment deltas those used for the production
4836   // of trackpoints (e.g. extracted from the UserInfo).
4837   // The prealignment deltas provided by user via config file must be already converted to target geometry:
4838   // this can be done externally using the macro ConvertDeltas.C
4839   //
4840   // delta_j_new = delta_j_old * Xj_old * Xj_new^-1
4841   // where X = Prod{delta_i,i=j-1:0} M_j
4842   // with j - the level of the alignable volume in the hierarchy, M - corresponding ideal matrix
4843   // Note that delta_j * Xj is equal to final (misaligned) matrix of corresponding geometry, G_j.
4844   // Since this method is used ONLY in the case where the prealignment deltas are equal to production deltas,
4845   // we have already loaded G_j_old in the fConvAlgMatOld (filled in the CacheMatricesOrig)
4846   // Hence, delta_j_new = G_j_old * Xj_new^-1
4847   //
4848   AliInfo("Converting deltas from initial to target geometry");
4849   int nMatOld = fConvAlgMatOld.GetEntriesFast(); // number of alignable matrices
4850   TClonesArray* deltArrNew = new TClonesArray("AliAlignObjParams",10);
4851   //
4852   TGeoHMatrix dmPar;
4853   int nDelNew = 0;
4854   //
4855   for (int im=0;im<nMatOld;im++) {
4856     TGeoHMatrix* mtGjold = (TGeoHMatrix*)fConvAlgMatOld[im];
4857     TString algname = mtGjold->GetTitle();
4858     UShort_t vID = AliITSAlignMille2Module::GetVolumeIDFromSymname(algname.Data());
4859     //
4860     // build X_new >>>
4861     TGeoHMatrix* parent = mtGjold;
4862     TGeoHMatrix xNew;
4863     int parID;
4864     while ( (parID=parent->GetUniqueID()-1)>=0 ) {
4865       parent = (TGeoHMatrix*)fConvAlgMatOld[parID];
4866       AliAlignObjParams* deltaPar = ConvFindDelta(deltArrNew,parent->GetTitle());
4867       if (deltaPar) deltaPar->GetMatrix(dmPar); xNew *= dmPar;
4868     }
4869     AliGeomManager::GetOrigGlobalMatrix(algname,dmPar); // ideal matrix of new geometry
4870     xNew *= dmPar;
4871     // build X_new <<<
4872     //
4873     dmPar  = *mtGjold;
4874     dmPar *= xNew.Inverse();
4875     new((*deltArrNew)[nDelNew++]) AliAlignObjParams(algname.Data(),vID,dmPar,kTRUE);
4876     //
4877   }
4878   delete fPrealignment;
4879   fPrealignment = deltArrNew;
4880   //
4881   // we don't need anymore old matrices
4882   fConvAlgMatOld.Delete();
4883   //
4884 }
4885
4886 //_______________________________________________________________________________________
4887 void AliITSAlignMille2::ConvSortHierarchically(TObjArray& matArr)
4888 {
4889   // Used only for the deltas conversion from one geometry to another
4890   // Sort the matrices according to hiearachy (coarse -> fine)
4891   //
4892   int nmat = matArr.GetEntriesFast();
4893   //
4894   for (int i=0;i<nmat;i++) {
4895     for (int j=i+1;j<nmat;j++) {
4896       TGeoHMatrix* matI = (TGeoHMatrix*) matArr[i];
4897       TGeoHMatrix* matJ = (TGeoHMatrix*) matArr[j];
4898       if (ConvIsJParentOfI(matI,matJ)) { // swap
4899         matArr[i] = matJ;
4900         matArr[j] = matI;
4901       }
4902     }
4903   }
4904   //
4905   // set direct parent id's in the UniqueID's
4906   for (int i=nmat;i--;) {
4907     TGeoHMatrix* matI = (TGeoHMatrix*) matArr[i];
4908     matI->SetUniqueID(0);
4909     for (int j=i;j--;) {
4910       TGeoHMatrix* matJ = (TGeoHMatrix*) matArr[j];
4911       if (ConvIsJParentOfI(matI,matJ)) { matI->SetUniqueID(j+1); break; }
4912     }
4913   }
4914 }
4915
4916 //_______________________________________________________________________________________
4917 Bool_t AliITSAlignMille2::ConvIsJParentOfI(const TGeoHMatrix* matI,const TGeoHMatrix* matJ) const
4918 {
4919   // Used only for the deltas conversion from one geometry to another
4920   // True if matJ is higher in hierarchy than
4921   //
4922   TString nmI = matI->GetTitle();
4923   TString nmJ = matJ->GetTitle();
4924   //
4925   int nlrI = nmI.CountChar('/');
4926   int nlrJ = nmJ.CountChar('/');
4927   if (nlrJ>=nlrI) return kFALSE;
4928   //
4929   // special case of SPD sectors
4930   if (nmI.BeginsWith("ITS/SPD1") && nmJ.BeginsWith("ITS/SPD0") && nlrJ==2) nmJ.ReplaceAll("SPD0","SPD1");
4931   return (nmI.BeginsWith(nmJ)) ? kTRUE:kFALSE;
4932   //
4933 }
4934
4935 //_______________________________________________________________________________________
4936 AliAlignObjParams* AliITSAlignMille2::ConvFindDelta(const TClonesArray* arrDelta,const TString& algname) const
4937 {
4938   // find the delta for given module
4939   if (!arrDelta) return 0;
4940   AliAlignObjParams* delta = 0;
4941   int nDeltas = arrDelta->GetEntries();
4942   for (int id=0;id<nDeltas;id++) {
4943     delta = (AliAlignObjParams*)arrDelta->At(id);
4944     if (algname==delta->GetSymName()) break;
4945     delta = 0;
4946   }
4947   return delta;
4948 }