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   if (md) return int(md->GetUniqueID());
2371   else return -1;
2372 }
2373
2374 //________________________________________________________________________________________________________
2375 Int_t AliITSAlignMille2::CheckCurrentTrack()
2376 {
2377   /// checks if AliTrackPoints belongs to defined modules
2378   /// return number of good poins
2379   /// return 0 if not enough points
2380
2381   Int_t npts = fTrack->GetNPoints();
2382   Int_t ngoodpts=0;
2383   // debug points
2384   for (int j=0; j<npts; j++) if (IsVIDContained(fTrack->GetVolumeID()[j])>=0) ngoodpts++;
2385   //
2386   if (ngoodpts<fMinNPtsPerTrack) return 0;
2387
2388   return ngoodpts;
2389 }
2390
2391 //________________________________________________________________________________________________________
2392 Int_t AliITSAlignMille2::ProcessTrack(const AliTrackPointArray *track, Double_t wgh)
2393 {
2394   /// Process track; Loop over hits and set local equations
2395   /// here 'track' is a AliTrackPointArray
2396   /// return 0 if success;
2397   //
2398   if (!fIsMilleInit) Init();
2399   //
2400   Int_t npts = track->GetNPoints();
2401   AliDebug(2,Form("*** Input track with %d points ***",npts));
2402
2403   // preprocessing of the input track: keep only points in defined volumes,
2404   // move points if prealignment is set, sort by Yglo if required
2405   fTrackWeight = wgh;
2406   fTrack=PrepareTrack(track);
2407   if (!fTrack) {
2408     RemoveHelixFitConstraint();
2409     RemoveVertexConstraint();
2410     return -1;
2411   }
2412   npts = fTrack->GetNPoints();
2413   if (npts>kMaxPoints) {
2414     AliError(Form("Compiled with kMaxPoints=%d, current track has %d points",kMaxPoints,npts));
2415   }
2416   AliDebug(2,Form("*** Processing prepared track with %d points ***",npts));
2417   //
2418   npts = FitTrack();
2419   if (npts<0) return npts;
2420   //
2421   //  printf("Params: "); for (int i=0;i<fNLocal;i++) printf("%+.2e ",fLocalInitParam[i]); printf("\n");//RRR
2422   Int_t nloceq=0;
2423   Int_t ngloeq=0;
2424   static Mille2Data md[kMaxPoints];
2425   //
2426   for (Int_t ipt=0; ipt<npts; ipt++) {
2427     fTrack->GetPoint(fCluster,ipt);
2428     fCluster.SetUniqueID(ipt+1);
2429     AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
2430
2431     // set geometry parameters for the the current module
2432     if (InitModuleParams()) continue;
2433     AliDebug(2,Form("    VolID=%d  Index=%d  InternalIdx=%d  symname=%s\n",
2434                     track->GetVolumeID()[ipt], fCurrentModule->GetIndex(),
2435                     fCurrentModule->GetUniqueID(), AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
2436     AliDebug(2,Form("    Preprocessed Point = ( %+f , %+f , %+f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
2437     int res = fTPAFitter ? AddLocalEquationTPA(md[nloceq]) : AddLocalEquation(md[nloceq]);
2438     if (res<0) {fTotBadLocEqPoints++; nloceq = 0; break;}
2439     else if (res==0) nloceq++;
2440     else {nloceq++; ngloeq++;}
2441   } // end loop over points
2442   //
2443   fTrack=NULL;
2444   // not enough good points?
2445   if (nloceq<fMinNPtsPerTrack || ngloeq<1) return -1;
2446   //
2447   // finally send local equations to millepede
2448   SetLocalEquations(md,nloceq);
2449   fMillepede->SaveRecordData(); // RRR
2450   //
2451   return 0;
2452 }
2453
2454 //________________________________________________________________________________________________________
2455 Int_t AliITSAlignMille2::FitTrack()
2456 {
2457   // Fit the track with selected constraints
2458   //
2459   const Double_t kfDiamondTolerance = 0.1;  //diamond tolerance on top of the MS error
2460   if (!fTrack) return -1;
2461   int npts = fTrack->GetNPoints();
2462   //
2463   if (fTPAFitter) {  // use dediacted fitter
2464     //
2465     // if the diamond point is attached, for the moment don't include it in the fit
2466     fTPAFitter->AttachPoints(fTrack,0, npts-1);
2467     fTPAFitter->SetBz(fBField);
2468     fTPAFitter->SetTypeCosmics(IsTypeCosmics());
2469     if (fIniTrackParamsMeth==1) fTPAFitter->SetIgnoreCov();
2470     //
2471     double chi2;
2472     double chi2f = 0;
2473     double dca2err;
2474     double dca2 = 0.;
2475     Bool_t fitIsDone = kFALSE;
2476     if (fUseDiamond && fDiamondPointID>0 && fCheckDiamondPoint==kDiamondCheckIfPrompt) { // diamond constraint was added, check if the track looks like prompt
2477       fTPAFitter->SetFirstLast(0,fDiamondPointID-1);
2478       if (IsCovIScaleTouched()) for (int i=npts;i--;) fTPAFitter->SetCovIScale(i,GetCovIScale(i));
2479       //
2480       chi2f = fTPAFitter->Fit(fConstrCharge,fConstrPT,fConstrPTErr);
2481       if ( chi2f<0 || (chi2f>fNStdDev*fStartFac && fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations()) ) { //RRR
2482         AliInfo(Form("Track fit failed on checking if it is prompt! skipping this track... Chi2:%+e",chi2f));
2483         fTPAFitter->Reset();
2484         //      fTrack = NULL;
2485         return -5;
2486       }
2487       double xyzRes[3];
2488       fTPAFitter->GetResiduals(xyzRes,&fDiamondI,kTRUE);
2489       dca2 = xyzRes[0]*xyzRes[0] + xyzRes[1]*xyzRes[1];
2490       double pT = IsFieldON() ? fTPAFitter->GetPt() : 0.45;
2491       if (pT<0.1) pT = 0.1;
2492       dca2err = kfDiamondTolerance + 0.02/pT;
2493       if (dca2>dca2err*dca2err) { // this is secondary
2494         int* clst = (int*) fTrack->GetClusterType();
2495         clst[fDiamondPointID] = -1;;
2496         fDiamondPointID = -1;
2497         fitIsDone = kTRUE;
2498         npts--;
2499       }
2500       else fTPAFitter->SetFirstLast(0,fDiamondPointID); // fit with diamond
2501     }
2502     //    fTPAFitter->SetParAxis(1);
2503     if (!fitIsDone) {
2504       if (IsCovIScaleTouched()) for (int i=npts;i--;) fTPAFitter->SetCovIScale(i,GetCovIScale(i));
2505       chi2 = fTPAFitter->Fit(fConstrCharge,fConstrPT,fConstrPTErr);
2506     }
2507     //
2508     RemoveHelixFitConstraint();  // suppress eventual constraints to not affect fit of the next track
2509     RemoveVertexConstraint(); // same ...
2510     //
2511     if ( !fitIsDone && (chi2<0 || (chi2>fNStdDev*fStartFac && fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations())) ) { //RRR
2512       AliInfo(Form("Track fit failed! skipping this track... Chi2:%+e",chi2));
2513       if (fUseDiamond && fDiamondPointID>0 && fCheckDiamondPoint==kDiamondCheckIfPrompt) AliInfo(Form("VertexFree fit gave Chi2:%+e with residual %+e",chi2f,TMath::Sqrt(dca2)));
2514       /*
2515         fTrack->Print("");
2516         fTPAFitter->FitHelixCrude();
2517         fTPAFitter->SetFitDone();
2518         fTPAFitter->Print();
2519       */
2520       fTPAFitter->Reset();
2521       //      fTrack = NULL;
2522       return -5;
2523     }
2524     fNLocal = fTPAFitter->IsFieldON() ? 5:4; // Attention: the fitter might have decided to work in line mode
2525     npts  = fTPAFitter->GetLast() - fTPAFitter->GetFirst() + 1; // actual number of points
2526     /*
2527       double *pr = fTPAFitter->GetParams();
2528       printf("FtPar: %+.5e  %+.5e  %+.5e  %+.5e | chi2:%.3e\n",pr[2],pr[0],pr[3],pr[1],chi2); // RRR
2529     */
2530   }
2531   else {
2532     //
2533     if (!fBOn) { // straight lines
2534       // set local starting parameters (to be substituted by ESD track parms)
2535       // local parms (fLocalInitParam[]) are:
2536       //      [0] = global x coord. of straight line intersection at y=0 plane
2537       //      [1] = global z coord. of straight line intersection at y=0 plane
2538       //      [2] = px/py
2539       //      [3] = pz/py
2540       InitTrackParams(fIniTrackParamsMeth);
2541       /*
2542       double *pr = fLocalInitParam;
2543       printf("FtPar: %+.5e  %+.5e  %+.5e  %+.5e |\n",pr[0],pr[1],pr[2],pr[3]); // RRR
2544       */
2545     }
2546     else {
2547       // local parms (fLocalInitParam[]) are the Riemann Fitter params
2548       if (!InitRiemanFit()) {
2549         AliInfo("Riemann fit failed! skipping this track...");
2550         fTrack=NULL;
2551         return -5;
2552       }
2553     }
2554   }
2555   return npts;
2556   //
2557 }
2558
2559 //________________________________________________________________________________________________________
2560 Int_t AliITSAlignMille2::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar)
2561 {
2562   /// calculate track intersection point in local coordinates
2563   /// according with a given set of parameters (local(4) and global(6))
2564   /// and fill fPintLoc/Glo
2565   ///    local are:   pgx0, pgz0, ugx, ugz   OR   riemann fitters pars
2566   ///    global are:  tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
2567   /// return 0 if success
2568
2569   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]));
2570   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]));
2571
2572
2573   // prepare the TGeoHMatrix
2574   TGeoHMatrix *tempHMat = fCurrentModule->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar,
2575                                                                            !fUseGlobalDelta);
2576   if (!tempHMat) return -1;
2577
2578   Double_t v0g[3]; // vector with straight line direction in global coord.
2579   Double_t p0g[3]; // point of the straight line (glo)
2580
2581   if (fBOn) { // B FIELD!
2582     AliTrackPoint prf;
2583     for (int ip=0; ip<5; ip++)
2584       fRieman->SetParam(ip,lpar[ip]);
2585
2586     if (!fRieman->GetPCA(fCluster,prf))  {
2587       AliInfo(Form("error in GetPCA for point %d",fCluster.GetVolumeID()));
2588       return -3;
2589     }
2590     // now determine straight line passing tangent to fit curve at prf
2591     // ugx = dX/dY_glo = DeltaX/DeltaY_glo
2592     // mo' P1=(X,Y,Z)_glo_prf
2593     //       => (x,y,Z)_trk_prf ruotando di alpha...
2594     Double_t alpha=fRieman->GetAlpha();
2595     Double_t x1g = prf.GetX();
2596     Double_t y1g = prf.GetY();
2597     Double_t z1g = prf.GetZ();
2598     Double_t x1t =  x1g*TMath::Cos(alpha) + y1g*TMath::Sin(alpha);
2599     Double_t y1t = -x1g*TMath::Sin(alpha) + y1g*TMath::Cos(alpha);
2600     Double_t z1t =  z1g;
2601
2602     Double_t x2t = x1t+1.0;
2603     Double_t y2t = y1t+fRieman->GetDYat(x1t);
2604     Double_t z2t = z1t+fRieman->GetDZat(x1t);
2605     Double_t x2g =  x2t*TMath::Cos(alpha) - y2t*TMath::Sin(alpha);
2606     Double_t y2g =  x2t*TMath::Sin(alpha) + y2t*TMath::Cos(alpha);
2607     Double_t z2g =  z2t;
2608
2609     AliDebug(3,Form("Riemann frame:  fAlpha = %+f  =  %+f  ",alpha,alpha*180./TMath::Pi()));
2610     AliDebug(3,Form("   prf_glo=( %+f , %+f , %+f )  prf_rf=( %+f , %+f , %+f )\n", x1g,y1g,z1g, x1t,y1t,z1t));
2611     AliDebug(3,Form("   mov_glo=( %+f , %+f , %+f )      rf=( %+f , %+f , %+f )\n",x2g,y2g,z2g, x2t,y2t,z2t));
2612
2613     if (TMath::Abs(y2g-y1g)<1e-15) {
2614       AliInfo("DeltaY=0! Cannot proceed...");
2615       return -1;
2616     }
2617     // ugx,1,ugz
2618     v0g[0] = (x2g-x1g)/(y2g-y1g);
2619     v0g[1]=1.0;
2620     v0g[2] = (z2g-z1g)/(y2g-y1g);
2621
2622     // point: just keep prf
2623     p0g[0]=x1g;
2624     p0g[1]=y1g;
2625     p0g[2]=z1g;
2626   }
2627   else { // staight line
2628     // vector of initial straight line direction in glob. coord
2629     v0g[0]=lpar[2];
2630     v0g[1]=1.0;
2631     v0g[2]=lpar[3];
2632
2633     // intercept in yg=0 plane in glob coord
2634     p0g[0]=lpar[0];
2635     p0g[1]=0.0;
2636     p0g[2]=lpar[1];
2637   }
2638   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]));
2639
2640   // same in local coord.
2641   Double_t p0l[3],v0l[3];
2642   tempHMat->MasterToLocalVect(v0g,v0l);
2643   tempHMat->MasterToLocal(p0g,p0l);
2644
2645   if (TMath::Abs(v0l[1])<1e-15) {
2646     AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
2647     return -1;
2648   }
2649
2650   // local intersection point
2651   fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
2652   fPintLoc[1] = 0;
2653   fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
2654
2655   // global intersection point
2656   tempHMat->LocalToMaster(fPintLoc,fPintGlo);
2657   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]));
2658
2659   return 0;
2660 }
2661
2662 //________________________________________________________________________________________________________
2663 Int_t AliITSAlignMille2::CalcDerivatives(Int_t paridx, Bool_t islpar)
2664 {
2665   /// calculate numerically (ROOT's style) the derivatives for
2666   /// local X intersection  and local Z intersection
2667   /// parlist: local  (islpar=kTRUE)  pgx0, pgz0, ugx0, ugz0  OR riemann's params
2668   ///          global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
2669   /// return 0 if success
2670
2671   // copy initial parameters
2672   Double_t lpar[kNLocal];
2673   Double_t gpar[kNParCh];
2674   Double_t *derivative;
2675   for (Int_t i=0; i<kNLocal; i++) lpar[i]=fLocalInitParam[i];
2676   for (Int_t i=0; i<kNParCh; i++) gpar[i]=fModuleInitParam[i];
2677
2678   // trial with fixed dpar...
2679   Double_t dpar = 0.;
2680
2681   if (islpar) { // track parameters
2682     //dpar=fLocalInitParam[paridx]*0.001;
2683     // set minimum dpar
2684     derivative = fDerivativeLoc[paridx];
2685     if (!fBOn) {
2686       if (paridx<3) dpar=1.0e-4; // translations
2687       else dpar=1.0e-6; // direction
2688     }
2689     else { // B Field
2690       // pepo: proviamo con 1/1000, poi evenctually 1/100...
2691       Double_t dfrac=0.01;
2692       switch(paridx) {
2693       case 0:
2694         // RMS cosmics: 1e-4
2695         dpar = TMath::Max(1.0e-6,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2696         break;
2697       case 1:
2698         // RMS cosmics: 0.2
2699         dpar = TMath::Max(0.002,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2700         break;
2701       case 2:
2702         // RMS cosmics: 9
2703         dpar = TMath::Max(0.09,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2704         break;
2705       case 3:
2706         // RMS cosmics: 7
2707         dpar = TMath::Max(0.07,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2708         break;
2709       case 4:
2710         // RMS cosmics: 0.3
2711         dpar = TMath::Max(0.003,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2712         break;
2713       }
2714     }
2715   }
2716   else { // alignment global parameters
2717     derivative = fDerivativeGlo[paridx];
2718     //dpar=fModuleInitParam[paridx]*0.001;
2719     if (paridx<3) dpar=1.0e-4; // translations
2720     else dpar=1.0e-2; // angles
2721   }
2722
2723   AliDebug(3,Form("+++ using dpar=%g",dpar));
2724
2725   // calculate derivative ROOT's like:
2726   //  using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
2727   Double_t pintl1[3]; // f(x-h)
2728   Double_t pintl2[3]; // f(x-h/2)
2729   Double_t pintl3[3]; // f(x+h/2)
2730   Double_t pintl4[3]; // f(x+h)
2731
2732   // first values
2733   if (islpar) lpar[paridx] -= dpar;
2734   else gpar[paridx] -= dpar;
2735   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2736   for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
2737
2738   // second values
2739   if (islpar) lpar[paridx] += dpar/2;
2740   else gpar[paridx] += dpar/2;
2741   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2742   for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
2743
2744   // third values
2745   if (islpar) lpar[paridx] += dpar;
2746   else gpar[paridx] += dpar;
2747   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2748   for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
2749
2750   // fourth values
2751   if (islpar) lpar[paridx] += dpar/2;
2752   else gpar[paridx] += dpar/2;
2753   if (CalcIntersectionPoint(lpar, gpar)) return -2;
2754   for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
2755
2756   Double_t h2 = 1./(2.*dpar);
2757   Double_t d0 = pintl4[0]-pintl1[0];
2758   Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
2759   derivative[0] = h2*(4*d2 - d0)/3.;
2760   if (TMath::Abs(derivative[0]) < 1.0e-9) derivative[0] = 0.0;
2761
2762   d0 = pintl4[2]-pintl1[2];
2763   d2 = 2.*(pintl3[2]-pintl2[2]);
2764   derivative[2] = h2*(4*d2 - d0)/3.;
2765   if (TMath::Abs(derivative[2]) < 1.0e-9) derivative[2]=0.0;
2766
2767   AliDebug(3,Form("\n+++ derivatives +++ \n"));
2768   AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",derivative[0]));
2769   AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",derivative[2]));
2770
2771   return 0;
2772 }
2773
2774 //________________________________________________________________________________________________________
2775 Int_t AliITSAlignMille2::AddLocalEquation(Mille2Data &m)
2776 {
2777   /// Define local equation for current cluster in X and Z coor.
2778   /// and store them to memory
2779   /// return -1 in case of failure to build some equation
2780   ///         0 if no free global parameters were found but local eq is built
2781   ///         1 if both local and global eqs are built
2782   //
2783   // store first intersection point
2784   if (CalcIntersectionPoint(fLocalInitParam, fModuleInitParam)) return -1;
2785   for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
2786
2787   AliDebug(2,Form("Intersect. point: L( %+f , %+f , %+f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
2788
2789   // calculate local derivatives numerically
2790   Bool_t zeroX = kTRUE;
2791   Bool_t zeroZ = kTRUE;
2792   //
2793   for (Int_t i=0; i<fNLocal; i++) {
2794     if (CalcDerivatives(i,kTRUE)) return -1;
2795     m.fDerLoc[i][kX] = fDerivativeLoc[i][0];
2796     m.fDerLoc[i][kZ] = fDerivativeLoc[i][2];
2797     if (zeroX) zeroX = IsZero(fDerivativeLoc[i][0]);
2798     if (zeroZ) zeroZ = IsZero(fDerivativeLoc[i][2]);
2799   }
2800   //  for (Int_t i=0; i<fNLocal; i++) AliDebug(2,Form("Local parameter %d - dXdpar = %g  - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
2801   //
2802   if (zeroX) {AliInfo("Skipping: zero local X derivatives!"); return -1;}
2803   if (zeroZ) {AliInfo("Skipping: zero local Z derivatives!"); return -1;}
2804   //
2805   int status = 0;
2806   int ifill = 0;
2807   //
2808   AliITSAlignMille2Module* endModule = fCurrentModule;
2809   //
2810   zeroX = zeroZ = kTRUE;
2811   Bool_t dfDone[kNParCh];
2812   for (int i=kNParCh;i--;) dfDone[i] = kFALSE;
2813   m.fNModFilled = 0;
2814   //
2815   // special block for SDD derivatives
2816   Double_t jacobian[kNParChGeom];
2817   Int_t nmodTested = 0;
2818   //
2819   do {
2820     if (fCurrentModule->GetNParFree()==0) continue;
2821     nmodTested++;
2822     for (Int_t i=0; i<kNParChGeom; i++) {   // common for all sensors: derivatives over geom params
2823       //
2824       if (!fUseGlobalDelta) dfDone[i] = kFALSE; // for global deltas the derivatives at diff. levels are different
2825       if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2826       if (!dfDone[i]) {
2827         if (CalcDerivatives(i,kFALSE)) return -1;
2828         else {
2829           dfDone[i] = kTRUE;
2830           if (zeroX) zeroX = IsZero(fDerivativeGlo[i][0]);
2831           if (zeroZ) zeroZ = IsZero(fDerivativeGlo[i][2]);
2832         }
2833       }
2834       //
2835       m.fDerGlo[ifill][kX] = fDerivativeGlo[i][0];
2836       m.fDerGlo[ifill][kZ] = fDerivativeGlo[i][2];
2837       m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
2838     }
2839     //
2840     // specific for special sensors
2841     Int_t sddLR = -1;
2842     if ( fCurrentModule->IsSDD() &&
2843          (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0  ||
2844           //      fCurrentModule->GetParOffset(sddLR = fMeasLoc[kX]>0 ?
2845           fCurrentModule->GetParOffset(sddLR = GetVDriftSDD()>0 ?
2846                                        AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR)>=0)
2847          ) {
2848       //
2849       // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
2850       // where V0 and T are the nominal drift velocity, time and time0
2851       // and the dT0 and dV are the corrections:
2852       // dX/dT0 = dX/dxloc * dxloc/dT0 = dX/dxloc * V0
2853       // dX/dV  = dX/dxloc * dxloc/dV =  dX/dxloc * (T-T0)
2854       // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
2855       //
2856       if (!dfDone[AliITSAlignMille2Module::kDOFT0] ||  !dfDone[sddLR]) {
2857         //
2858         double dXdxlocsens=0., dZdxlocsens=0.;
2859         //
2860         // if the current module is the sensor itself and we work with local params, then
2861         // we can directly take dX/dxloc_sens dZ/dxloc_sens
2862         if (!fUseGlobalDelta && fCurrentModule->GetVolumeID()==fCluster.GetVolumeID()) {
2863           if (!dfDone[AliITSAlignMille2Module::kDOFTX]) {
2864             CalcDerivatives(AliITSAlignMille2Module::kDOFTX,kFALSE);
2865             dfDone[AliITSAlignMille2Module::kDOFTX] = kTRUE;
2866           }
2867           dXdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][0];
2868           dZdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][2];
2869         }
2870         //
2871         else { // need to perform some transformations
2872           // fetch the jacobian of the transformation from the sensors local frame to the frame
2873           // where the parameters are defined:
2874           // Global: dX/dxloc_sens = dX/dxgl*dxgl/dxloc_sens + ...dX/dphigl*dphigl/dxloc_sens
2875           if (fUseGlobalDelta) fCurrentModule->CalcDerivGloLoc(fCluster.GetVolumeID(),
2876                                                                AliITSAlignMille2Module::kDOFTX, jacobian);
2877           // Local:  dX/dxloc_sens = dX/dxcurr*dxcurr/dxloc_sens +..+dX/dphicurr * dphicurr/dxloc_sens
2878           else                 fCurrentModule->CalcDerivCurLoc(fCluster.GetVolumeID(),
2879                                                                AliITSAlignMille2Module::kDOFTX, jacobian);
2880           //
2881           for (int j=0;j<kNParChGeom;j++) {
2882             // need global derivative even if the j-th param is locked
2883             if (!dfDone[j]) {CalcDerivatives(j,kFALSE); dfDone[j] = kTRUE;}
2884             dXdxlocsens += fDerivativeGlo[j][0] * jacobian[j];
2885             dZdxlocsens += fDerivativeGlo[j][2] * jacobian[j];
2886           }
2887         }
2888         //
2889         if (zeroX) zeroX = IsZero(dXdxlocsens);
2890         if (zeroZ) zeroZ = IsZero(dZdxlocsens);
2891         //
2892         double vdrift = GetVDriftSDD();
2893         double tdrift = GetTDriftSDD();
2894         //
2895         fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0] = dXdxlocsens*vdrift;
2896         fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2] = dZdxlocsens*vdrift;
2897         dfDone[AliITSAlignMille2Module::kDOFT0] = kTRUE;
2898         //
2899         double mltCorr = fIsSDDVDriftMult ? TMath::Abs(vdrift) : 1;
2900         fDerivativeGlo[sddLR][0] = -dXdxlocsens*mltCorr*TMath::Sign(tdrift,vdrift);
2901         fDerivativeGlo[sddLR][2] = -dZdxlocsens*mltCorr*TMath::Sign(tdrift,vdrift);
2902         dfDone[sddLR] = kTRUE;
2903         //
2904       }
2905       //
2906       if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
2907         m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0];
2908         m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2];
2909         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);
2910       }
2911       //
2912       if (fCurrentModule->GetParOffset(sddLR)>=0) {
2913         m.fDerGlo[ifill][kX] = fDerivativeGlo[sddLR][0];
2914         m.fDerGlo[ifill][kZ] = fDerivativeGlo[sddLR][2];
2915         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(sddLR);
2916       }
2917     }
2918     //
2919     m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
2920   } while( (fCurrentModule=fCurrentModule->GetParent()) );
2921   //
2922   if (nmodTested>0 && zeroX) {AliInfo("Skipping: zero global X derivatives!");return -1;}
2923   if (nmodTested>0 && zeroZ) {AliInfo("Skipping: zero global Z derivatives!");return -1;}
2924   //
2925   // ok, can copy to m
2926   AliDebug(2,Form("Adding local equation X with fMeas=%.6f  and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
2927   m.fMeas[kX] = fMeasLoc[0]-fPintLoc0[0];
2928   m.fSigma[kX] = fSigmaLoc[0];
2929   //
2930   AliDebug(2,Form("Adding local equation Z with fMeas=%.6f  and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
2931   m.fMeas[kZ] = fMeasLoc[2]-fPintLoc0[2];
2932   m.fSigma[kZ] = fSigmaLoc[2];
2933   //
2934   m.fNGlobFilled = ifill;
2935   fCurrentModule = endModule;
2936   //
2937   status += Int_t(!zeroX && !zeroZ); // 0 - only locals, 1 locals + globals
2938   return status;
2939 }
2940
2941 //________________________________________________________________________________________________________
2942 Int_t AliITSAlignMille2::AddLocalEquationTPA(Mille2Data &m)
2943 {
2944   /// Define local equation for current cluster in X Y and Z coor.
2945   /// and store them to memory
2946   /// return -1 in case of failure to build some equation
2947   ///         0 if no free global parameters were found but local eq is built
2948   ///         1 if both local and global eqs are built
2949   //
2950   static int cnt = 0;
2951   Bool_t dbg = kFALSE;//kTRUE;
2952   if (++cnt>100000) dbg = kFALSE;
2953
2954   int curpoint = fCluster.GetUniqueID()-1;
2955   TGeoHMatrix *tempHMat = GetSensorCurrMatrixSID(fCurrentSensID);// fCurrentModule->GetSensitiveVolumeMatrix(fCluster.GetVolumeID());
2956   //
2957   fTPAFitter->GetDResDParams(&fDerivativeLoc[0][0], curpoint);    // resid. derivatives over the track parameters
2958   for (Int_t i=fNLocal; i--;) tempHMat->MasterToLocalVect(fDerivativeLoc[i],m.fDerLoc[i]);
2959   //
2960   int status = 0;
2961   // derivatives over the global parameters ---------------------------------------->>>
2962   Double_t dGL[3];     // derivative of global position vs local X (for SDD)
2963   Double_t dRdP[3][3]; // derivative of local residuals vs local position
2964   Double_t dPdG[AliITSAlignMille2Module::kMaxParGeom][3]; // derivatives of local position vs global params
2965   fTPAFitter->GetDResDPos(&fDerivativeGlo[0][0], curpoint);
2966   if (fCurrentSensID!=kVtxSensID) for (int i=3;i--;) tempHMat->MasterToLocalVect(fDerivativeGlo[i],dRdP[i]);
2967   else for (int i=3;i--;) for (int j=3;j--;) dRdP[i][j] = fDerivativeGlo[i][j];  // vertex constraint is in lab
2968   //
2969   if (dbg) {
2970     printf("\nCurrentMod: %s Sens:%d\n",fCurrentModule->GetName(),fCurrentSensID); //RRR
2971     printf("Module Matrix: ");
2972     fCurrentModule->GetMatrix()->Print(); //RRR
2973     for (int i=0;i<3;i++) {
2974       printf("dRdP[M%d][resI] ",i); for (int j=0;j<3;j++) printf(":[%d] %+.3e ",j,dRdP[i][j]); printf("\n");
2975     }//RRR
2976     printf("Sensor Matrix: "); tempHMat->Print();
2977   }
2978   UInt_t ifill=0, dfDone = 0;
2979   m.fNModFilled = 0;
2980   //
2981   AliITSAlignMille2Module* endModule = fCurrentModule;
2982   //
2983   m.fModuleID[0] = fCurrentModule->GetUniqueID(); // always register id of the base module, even if it has no DOF
2984   //
2985   do {
2986     if (fCurrentModule->GetNParFree()==0) continue;
2987     status = 1;
2988     if (!fUseGlobalDelta) dfDone = 0; // for local deltas the derivatives at diff. levels are different
2989     Bool_t jacobOK = kFALSE;
2990     //
2991     for (Int_t i=0; i<kNParChGeom; i++) {              // common for all sensors: derivatives over geom params
2992       if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2993       //
2994       if (!TestWordBit(dfDone,i)) {                    // need to calculate new derivative
2995         if (!jacobOK) {
2996           if (fCurrentSensID!=kVtxSensID) {
2997             fCurrentModule->CalcDerivDPosDPar(fCluster.GetVolumeID(),fMeasLoc,&dPdG[0][0]);
2998             if (dbg) {
2999               for (int i1=0;i1<3;i1++) {
3000                 printf("Jacob:dPdG[gpar%d][Mj]",i1); for (int j1=0;j1<3;j1++) printf(":[%d] %+.3e ",j1,dPdG[i1][j1]); printf("\n");//RRR
3001               }
3002             }
3003           }
3004           else {
3005             // this is a vertex constraint: only lateral shifts are allowed, no rotations
3006             for (int ip=AliITSAlignMille2Module::kMaxParGeom;ip--;) for (int jp=3;jp--;) dPdG[ip][jp] = (ip==jp) ? 1:0;
3007           }
3008           jacobOK = kTRUE;
3009         }
3010         // dRes_j/dGlo_i = \sum_{k=1:3}  dRes_j/dPos_k * dPos_k/dGlo_i
3011         fDerivativeGlo[i][kX] = dRdP[kX][kX]*dPdG[i][kX] + dRdP[kY][kX]*dPdG[i][kY] + dRdP[kZ][kX]*dPdG[i][kZ];
3012         fDerivativeGlo[i][kY] = dRdP[kX][kY]*dPdG[i][kX] + dRdP[kY][kY]*dPdG[i][kY] + dRdP[kZ][kY]*dPdG[i][kZ];
3013         fDerivativeGlo[i][kZ] = dRdP[kX][kZ]*dPdG[i][kX] + dRdP[kY][kZ]*dPdG[i][kY] + dRdP[kZ][kZ]*dPdG[i][kZ];
3014         SetWordBit(dfDone,i);
3015       }
3016       //
3017       if (dbg) {
3018         printf("Level %s DGlob[par%d][resJ] ",fCurrentModule->GetName(),i); //RRR
3019         for (int k=0;k<3;k++) printf(":[%d] %+.3e ",k, fDerivativeGlo[i][k]); printf("\n");//RRR
3020       }
3021       m.fDerGlo[ifill][kX] = fDerivativeGlo[i][kX];
3022       m.fDerGlo[ifill][kY] = fDerivativeGlo[i][kY];
3023       m.fDerGlo[ifill][kZ] = fDerivativeGlo[i][kZ];
3024       m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
3025       //
3026     }
3027     //
3028     if ( fCurrentModule->IsSDD() ) {     // specific for SDD
3029       //
3030       // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
3031       // where V0 and T are the nominal drift velocity, time and time0
3032       // and the dT0 and dV are the corrections:
3033       // drloc_i/dT0 = sum_j drloc_i/dMeasGlo_j * dMeasGlo_j/dT0 =
3034       //             = sum_j drloc_i/dMeasGlo_j sum_k dMeasGlo_j/dMeasLoc_k * dMeasLoc_k/dT0
3035       //             = sum_j drloc_i/dMeasGlo_j dMeasGlo_j/dMeasLoc_X * V0
3036       //
3037       // drloc_i/dV0 = sum_j drloc_i/dMeasGlo_j * dMeasGlo_j/dV0 =
3038       //             = sum_j drloc_i/dMeasGlo_j sum_k dMeasGlo_j/dMeasLoc_k * dMeasLoc_k/dV0
3039       //             = sum_j drloc_i/dMeasGlo_j dMeasGlo_j/dMeasLoc_X * T0
3040
3041       // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
3042       //
3043       Bool_t jacOK = kFALSE;
3044       //Int_t sddLR = fMeasLoc[kX]>0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR;
3045       Int_t sddLR = GetVDriftSDD()>0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR;
3046       if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
3047         if (!TestWordBit(dfDone, AliITSAlignMille2Module::kDOFT0)) {
3048           double vdrift = GetVDriftSDD();
3049           JacobianPosGloLoc(kX,dGL);
3050           jacOK = kTRUE;
3051           fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kX] =
3052             vdrift*(dRdP[kX][kX]*dGL[kX] + dRdP[kY][kX]*dGL[kY] + dRdP[kZ][kX]*dGL[kZ]);
3053           fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kY] =
3054             vdrift*(dRdP[kX][kY]*dGL[kX] + dRdP[kY][kY]*dGL[kY] + dRdP[kZ][kY]*dGL[kZ]);
3055           fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kZ] =
3056             vdrift*(dRdP[kX][kZ]*dGL[kX] + dRdP[kY][kZ]*dGL[kY] + dRdP[kZ][kZ]*dGL[kZ]);
3057           //
3058           SetWordBit(dfDone, AliITSAlignMille2Module::kDOFT0);
3059         }
3060         m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kX];
3061         m.fDerGlo[ifill][kY] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kY];
3062         m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kZ];
3063         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);
3064       }
3065       //
3066       if (fCurrentModule->GetParOffset(sddLR)>=0) {
3067         if (!TestWordBit(dfDone, sddLR)) {
3068           double tdrift = TMath::Sign(GetTDriftSDD(), GetVDriftSDD());
3069           double vdrift = fIsSDDVDriftMult ? TMath::Abs(GetVDriftSDD()) : 1;
3070           if (!jacOK) JacobianPosGloLoc(kX,dGL);
3071           fDerivativeGlo[sddLR][kX] =
3072             -tdrift*vdrift*(dRdP[kX][kX]*dGL[kX] + dRdP[kY][kX]*dGL[kY] + dRdP[kZ][kX]*dGL[kZ]);
3073           fDerivativeGlo[sddLR][kY] =
3074             -tdrift*vdrift*(dRdP[kX][kY]*dGL[kX] + dRdP[kY][kY]*dGL[kY] + dRdP[kZ][kY]*dGL[kZ]);
3075           fDerivativeGlo[sddLR][kZ] =
3076             -tdrift*vdrift*(dRdP[kX][kZ]*dGL[kX] + dRdP[kY][kZ]*dGL[kY] + dRdP[kZ][kZ]*dGL[kZ]);
3077           SetWordBit(dfDone, sddLR);
3078         }
3079         m.fDerGlo[ifill][kX] = fDerivativeGlo[sddLR][kX];
3080         m.fDerGlo[ifill][kY] = fDerivativeGlo[sddLR][kY];
3081         m.fDerGlo[ifill][kZ] = fDerivativeGlo[sddLR][kZ];
3082         m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(sddLR);
3083       }
3084     }
3085     //
3086     m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
3087   } while( (fCurrentModule=fCurrentModule->GetParent()) );
3088   //
3089   // store first local residuals
3090   fTPAFitter->GetResiduals(fPintLoc , curpoint);       // lab residuals
3091   for (int i=3;i--;) fPintLoc[i] = -fPintLoc[i];
3092   if   (fCurrentSensID!=kVtxSensID) tempHMat->MasterToLocalVect(fPintLoc,m.fMeas);       // local residuals
3093   else for (int i=3;i--;) m.fMeas[i] = fPintLoc[i];
3094   if (dbg) {
3095     printf("res(meas-loc) "); for (int k=0;k<3;k++) printf(":[%d] %+.3e ",k,m.fMeas[k]); printf("\n");
3096     printf("Fin:%s %+e %+e\n",endModule->GetName(), fDerivativeGlo[kZ][kZ], fPintLoc[kZ]);
3097   }//RRR
3098   m.fSigma[kX] = fSigmaLoc[kX];
3099   m.fSigma[kY] = fSigmaLoc[kY];
3100   m.fSigma[kZ] = fSigmaLoc[kZ];
3101   //
3102   m.fNGlobFilled = ifill;
3103   fCurrentModule = endModule;
3104   //
3105   return status;
3106 }
3107
3108 //________________________________________________________________________________________________________
3109 void AliITSAlignMille2::SetLocalEquations(const Mille2Data *marr, Int_t neq)
3110 {
3111   /// Set local equations with data stored in m
3112   /// return 0 if success
3113   //
3114   for (Int_t j=0; j<neq; j++) {
3115     //
3116     const Mille2Data &m = marr[j];
3117     //
3118     Bool_t filled = kFALSE;
3119     for (int ic=3;ic--;) {
3120       // for the diamond point (if any) the Y residual is accounted
3121       if (ic==kY && !fUseLocalYErr && !(m.fModuleID[0]==fDiamondModID)) continue;
3122       AliDebug(2,Form("setting local equation %c with fMeas=%.6f  and fSigma=%.6f",fgkXYZ[ic],m.fMeas[ic], m.fSigma[ic]));
3123       Int_t nzero = 0;
3124       for (int i=fNLocal; i--;) nzero += SetLocalDerivative(i,m.fDerLoc[i][ic] );
3125       if (nzero==fNLocal) {
3126         AliInfo(Form("Skipping %c residual due to the zero derivatives!",fgkXYZ[ic]));
3127         continue;
3128       }
3129       for (int i=m.fNGlobFilled;i--;) SetGlobalDerivative( m.fParMilleID[i] , m.fDerGlo[i][ic] );
3130       fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.fMeas[ic], m.fSigma[ic]);
3131       filled = kTRUE;
3132       //
3133     }
3134     //
3135     if (filled) for (int i=m.fNModFilled;i--;) GetMilleModule(m.fModuleID[i])->IncNProcessedPoints();
3136   }
3137   //
3138   double wgh = 1;
3139   if (GetWeightPt() && fTPAFitter) {
3140     wgh = fTPAFitter->GetPt();
3141     if (wgh>10) wgh = 10.;
3142     if (wgh<0) wgh = fTPAFitter->IsTypeCosmics() ? 7 : 0.5;
3143     if (GetWeightPt()>0) wgh = TMath::Power(wgh,GetWeightPt());
3144   }
3145   fMillepede->SetRecordWeight(wgh*fTrackWeight);
3146   fMillepede->SetRecordRun(fRunID);
3147   //
3148 }
3149
3150 //________________________________________________________________________________________________________
3151 Int_t AliITSAlignMille2::GlobalFit()
3152 {
3153   /// Call global fit; Global parameters are stored in parameters
3154   if (!fIsMilleInit) Init();
3155   //
3156   ApplyPreConstraints();
3157   int res = fMillepede->GlobalFit();
3158   AliInfo(Form("%s fitting global parameters!",res ? "Done":"Failed"));
3159   if (res) {
3160     // fetch the parameters
3161     for (int imd=fNModules;imd--;) {
3162       AliITSAlignMille2Module* mod = GetMilleModule(imd);
3163       int nprocp = 0;
3164       for (int ip=mod->GetNParTot();ip--;) {
3165         int idp = mod->GetParOffset(ip);
3166         if (idp<0) continue;    // was not in the explicit fit
3167         mod->SetParVal(ip,fMillepede->GetFinalParam(idp));
3168         mod->SetParErr(ip,fMillepede->GetFinalError(idp));
3169         int np = fMillepede->GetProcessedPoints(idp);
3170         if (TMath::Abs(np)>TMath::Abs(nprocp)) nprocp = np;
3171       }
3172       if (!mod->GetNProcessedPoints()) mod->SetNProcessedPoints(nprocp);
3173     }
3174
3175   }
3176   ApplyPostConstraints();
3177   return res;
3178 }
3179
3180 //________________________________________________________________________________________________________
3181 void AliITSAlignMille2::PrintGlobalParameters()
3182 {
3183   /// Print global parameters
3184   if (!fIsMilleInit) {
3185     AliInfo("Millepede has not been initialized!");
3186     return;
3187   }
3188   fMillepede->PrintGlobalParameters();
3189 }
3190
3191 //________________________________________________________________________________________________________
3192 Int_t AliITSAlignMille2::LoadSuperModuleFile(const Char_t *sfile)
3193 {
3194   // load definitions of supermodules from a root file
3195   // return 0 if success
3196   AliInfo(Form("Loading SuperModule definitions from %s",sfile));
3197   TFile *smf=TFile::Open(sfile);
3198   if (!smf->IsOpen()) {
3199     AliInfo(Form("Cannot open supermodule file %s",sfile));
3200     return -1;
3201   }
3202
3203   TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules");
3204   if (!sma) {
3205     AliInfo(Form("Cannot find ITSMilleSuperModules array in file"));
3206     return -2;
3207   }
3208   Int_t nsma=sma->GetEntriesFast();
3209   AliInfo(Form("Array of SuperModules with %d entries\n",nsma));
3210   //
3211   // pepo200709
3212   Char_t st[2048];
3213   char symname[250];
3214   // end pepo200709
3215
3216   UShort_t volid;
3217   TGeoHMatrix m;
3218   //
3219   for (Int_t i=0; i<nsma; i++) {
3220     AliAlignObjParams *a = (AliAlignObjParams*)sma->UncheckedAt(i);
3221     volid=a->GetVolUID();
3222     strncpy(st,a->GetSymName(),TMath::Min(sizeof(st),strlen(a->GetSymName())+1));
3223     a->GetMatrix(m);
3224     //
3225     symname[0] = '\0';
3226     sscanf(st,"%249s",symname);
3227     //
3228     // decode module list
3229     char *stp=strstr(st,"ModuleList:");
3230     if (!stp) return -3;
3231     stp += 11;
3232     int idx[2200];
3233     char spp[200]; int jp=0;
3234     char cl[20];
3235     strncpy(st,stp,TMath::Min(sizeof(st),strlen(stp)+1));
3236     int l=strlen(st);
3237     int j=0;
3238     int n=0;
3239     //
3240     while (j<=l) {
3241       if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) {
3242         spp[jp]=0;
3243         jp=0;
3244         if (strlen(spp)) {
3245           int k=strcspn(spp,"-");
3246           if (k<int(strlen(spp))) { // c'e' il -
3247             strncpy(cl,&(spp[k+1]), TMath::Min(sizeof(cl),strlen(&spp[k+1])+1));
3248             spp[k]=0;
3249             int ifrom=atoi(spp); int ito=atoi(cl);
3250             for (int b=ifrom; b<=ito; b++) {
3251               idx[n]=b;
3252               n++;
3253             }
3254           }
3255           else { // numerillo singolo
3256             idx[n]=atoi(spp);
3257             n++;
3258           }
3259         }
3260       }
3261       else {
3262         spp[jp]=st[j];
3263         jp++;
3264       }
3265       j++;
3266     }
3267     UShort_t volidsv[2198];
3268     for (j=0;j<n;j++) {
3269       volidsv[j]=AliITSAlignMille2Module::GetVolumeIDFromIndex(idx[j]);
3270       if (!volidsv[j]) {
3271         AliInfo(Form("Index %d not valid (range 0->%d)",idx[j],kMaxITSSensID));
3272         return -5;
3273       }
3274     }
3275     Int_t smindex=int(2198+volid-14336); // virtual index
3276     //
3277     fSuperModule.AddAtAndExpand(new AliITSAlignMille2Module(smindex,volid,symname,&m,n,volidsv),fNSuperModules);
3278     //
3279     fNSuperModules++;
3280   }
3281   //
3282   smf->Close();
3283   //
3284   return 0;
3285 }
3286
3287 //________________________________________________________________________________________________________
3288 void AliITSAlignMille2::ConstrainModuleSubUnitsMean(Int_t idm, Double_t val, UInt_t pattern)
3289 {
3290   // require that sum of modifications for the childs of this module is = val, i.e.
3291   // the internal corrections moves the module as a whole by fixed value (0 by default).
3292   // pattern is the bit pattern for the parameters to constrain
3293   //
3294   if (fIsMilleInit) {
3295     AliInfo("Millepede has been already initialized: no constrain may be added!");
3296     return;
3297   }
3298   if (!GetMilleModule(idm)->GetNChildren()) return;
3299   TString nm = "cstrSUMean";
3300   nm += GetNConstraints();
3301   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
3302                                                                       idm,val,pattern);
3303   cstr->SetConstraintID(GetNConstraints());
3304   fConstraints.Add(cstr);
3305 }
3306
3307 //________________________________________________________________________________________________________
3308 void AliITSAlignMille2::ConstrainModuleSubUnitsMedian(Int_t idm, Double_t val, UInt_t pattern)
3309 {
3310   // require that median of the modifications for the childs of this module is = val, i.e.
3311   // the internal corrections moves the module as a whole by fixed value (0 by default)
3312   // module the outliers.
3313   // pattern is the bit pattern for the parameters to constrain
3314   // The difference between the mean and the median will be transfered to the parent
3315   if (fIsMilleInit) {
3316     AliInfo("Millepede has been already initialized: no constrain may be added!");
3317     return;
3318   }
3319   if (!GetMilleModule(idm)->GetNChildren()) return;
3320   TString nm = "cstrSUMed";
3321   nm += GetNConstraints();
3322   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
3323                                                                       idm,val,pattern);
3324   cstr->SetConstraintID(GetNConstraints());
3325   fConstraints.Add(cstr);
3326 }
3327
3328 //________________________________________________________________________________________________________
3329 void AliITSAlignMille2::ConstrainOrphansMean(Double_t val, UInt_t pattern)
3330 {
3331   // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3332   // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3333   // pattern is the bit pattern for the parameters to constrain
3334   //
3335   if (fIsMilleInit) {
3336     AliInfo("Millepede has been already initialized: no constrain may be added!");
3337     return;
3338   }
3339   TString nm = "cstrOMean";
3340   nm += GetNConstraints();
3341   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
3342                                                                       -1,val,pattern);
3343   cstr->SetConstraintID(GetNConstraints());
3344   fConstraints.Add(cstr);
3345 }
3346
3347 //________________________________________________________________________________________________________
3348 void AliITSAlignMille2::ConstrainOrphansMedian(Double_t val, UInt_t pattern)
3349 {
3350   // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3351   // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3352   // pattern is the bit pattern for the parameters to constrain
3353   //
3354   if (fIsMilleInit) {
3355     AliInfo("Millepede has been already initialized: no constrain may be added!");
3356     return;
3357   }
3358   TString nm = "cstrOMed";
3359   nm += GetNConstraints();
3360   AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
3361                                                                       -1,val,pattern);
3362   cstr->SetConstraintID(GetNConstraints());
3363   fConstraints.Add(cstr);
3364 }
3365
3366 //________________________________________________________________________________________________________
3367 void AliITSAlignMille2::ConstrainLocal(const Char_t* name,Double_t *parcf,Int_t npar,Double_t val,Double_t err)
3368 {
3369   // apply constraint on parameters in the local frame
3370   if (fIsMilleInit) {
3371     AliInfo("Millepede has been already initialized: no constrain may be added!");
3372     return;
3373   }
3374   AliITSAlignMille2ConstrArray *cstr = new AliITSAlignMille2ConstrArray(name,parcf,npar,val,err);
3375   cstr->SetConstraintID(GetNConstraints());
3376   fConstraints.Add(cstr);
3377 }
3378
3379 //________________________________________________________________________________________________________
3380 void AliITSAlignMille2::ApplyGaussianConstraint(const AliITSAlignMille2ConstrArray* cstr)
3381 {
3382   // apply the constraint on the local corrections of a list of modules
3383   int nmod = cstr->GetNModules();
3384   double jacobian[AliITSAlignMille2Module::kMaxParGeom][AliITSAlignMille2Module::kMaxParGeom];
3385   //
3386   // check if this not special SDDT0 constraint
3387   if (cstr->GetPattern()==BIT(AliITSAlignMille2Module::kDOFT0)) {
3388     for (int i=0;i<cstr->GetNModules()-1;i++) {
3389       AliITSAlignMille2Module *mdI = GetMilleModule(cstr->GetModuleID(i));
3390       if (!mdI->IsFreeDOF(AliITSAlignMille2Module::kDOFT0)) continue;
3391       for (int j=i+1;j<cstr->GetNModules();j++) {
3392         AliITSAlignMille2Module *mdJ = GetMilleModule(cstr->GetModuleID(j));
3393         if (!mdJ->IsFreeDOF(AliITSAlignMille2Module::kDOFT0)) continue;
3394         //
3395         ResetLocalEquation();
3396         fGlobalDerivatives[mdI->GetParOffset(AliITSAlignMille2Module::kDOFT0)] = 1;
3397         fGlobalDerivatives[mdJ->GetParOffset(AliITSAlignMille2Module::kDOFT0)] =-1;
3398         AddConstraint(fGlobalDerivatives, 0, 1.E-6);
3399       }
3400     }
3401     return;
3402   }
3403
3404   for (int imd=nmod;imd--;) {
3405     int modID = cstr->GetModuleID(imd);
3406     AliITSAlignMille2Module* mod = GetMilleModule(modID);
3407     ResetLocalEquation();
3408     int nadded = 0;
3409     double value = cstr->GetValue();
3410     double sigma = cstr->GetError();
3411     //
3412     // in case the reference (survey) deltas were imposed for Gaussian constraints
3413     // already accumulated corrections: they must be subtracted from the constraint value.
3414     if (IsConstraintWrtRef()) {
3415       //
3416       Double_t precal[AliITSAlignMille2Module::kMaxParTot];
3417       Double_t refcal[AliITSAlignMille2Module::kMaxParTot];
3418       for (int ip=AliITSAlignMille2Module::kMaxParTot;ip--;) {precal[ip]=0; refcal[ip] = 0.;}
3419       //
3420       // check if there was a reference delta provided for this module
3421       AliAlignObjParams* parref = GetConstrRefObject(mod->GetName());
3422       if (parref) parref->GetPars(refcal, refcal+3);    // found reference delta
3423       //
3424       // extract already applied local corrections for this module
3425       if (fPrealignment) {
3426         //
3427         AliAlignObjParams *preo = GetPrealignedObject(mod->GetName());
3428         if (preo) {
3429           TGeoHMatrix preMat,tmpMat = *mod->GetMatrix(); //  Delta_Glob * Delta_Glob_Par * M
3430           preo->GetMatrix(preMat);                       //  Delta_Glob
3431           preMat.MultiplyLeft( &tmpMat.Inverse() );      //  M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
3432           tmpMat.MultiplyLeft( &preMat );                //  (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
3433           AliAlignObjParams algob;
3434           algob.SetMatrix(tmpMat);
3435           algob.GetPars(precal,precal+3); // local corrections for geometry
3436         }
3437       }
3438       //
3439       // subtract the contribution to constraint from precalibration
3440       for (int ipar=cstr->GetNCoeffs();ipar--;) value += (refcal[ipar]-precal[ipar])*cstr->GetCoeff(ipar);
3441       //
3442     }
3443     //
3444     if (fUseGlobalDelta) mod->CalcDerivLocGlo(&jacobian[0][0]);
3445     //
3446     for (int ipar=cstr->GetNCoeffs();ipar--;) {
3447       double coef = cstr->GetCoeff(ipar);
3448       if (IsZero(coef)) continue;
3449       //
3450       if (!fUseGlobalDelta || ipar>= AliITSAlignMille2Module::kMaxParGeom) { //
3451         // we are working with local params or if the given param is not related to geometry,
3452         // apply the constraint directly
3453         int parPos = mod->GetParOffset(ipar);
3454         if (parPos<0) continue; // not in the fit
3455         fGlobalDerivatives[parPos] += coef;
3456         nadded++;
3457       }
3458       else { // we are working with global params, while the constraint is on local ones -> jacobian
3459         for (int jpar=AliITSAlignMille2Module::kMaxParGeom;jpar--;) {
3460           int parPos = mod->GetParOffset(jpar);
3461           if (parPos<0) continue;
3462           fGlobalDerivatives[parPos] += coef*jacobian[ipar][jpar];
3463           nadded++;
3464         }
3465       }
3466     }
3467     if (nadded) AddConstraint(fGlobalDerivatives, value, sigma);
3468   }
3469   //
3470 }
3471
3472 //________________________________________________________________________________________________________
3473 void AliITSAlignMille2::ApplyPreConstraints()
3474 {
3475   // apply constriants which cannot be imposed after the fit
3476   int nconstr = GetNConstraints();
3477   for (int i=0;i<nconstr;i++) {
3478     AliITSAlignMille2Constraint* cstr = GetConstraint(i);
3479     //
3480     if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) {
3481       ApplyGaussianConstraint( (AliITSAlignMille2ConstrArray*)cstr);
3482       continue;
3483     }
3484     //
3485     if (cstr->GetType() == AliITSAlignMille2Constraint::kTypeMedian) continue; // post type constraint
3486     //
3487     if (!fUseGlobalDelta) continue; // mean/med constraints must be applied to global deltas
3488     // apply constraint on the mean's before the fit
3489     int imd = cstr->GetModuleID();
3490     if (imd>=0) {
3491       AliITSAlignMille2Module* mod = GetMilleModule(imd);
3492       UInt_t pattern = 0;
3493       for (int ipar=mod->GetNParTot();ipar--;) {
3494         if (!cstr->IncludesParam(ipar)) continue;
3495         if (mod->GetParOffset(ipar)<0) continue; // parameter is not in the explicit fit -> post constraint
3496         pattern |= 0x1<<ipar;
3497         cstr->SetApplied(ipar);
3498       }
3499       ConstrainModuleSubUnits(imd,cstr->GetValue(),pattern);
3500       //
3501     }
3502     else if (!PseudoParentsAllowed()) {
3503       ConstrainOrphans(cstr->GetValue(),(UInt_t)cstr->GetPattern());
3504       cstr->SetApplied(-1);
3505     }
3506   }
3507   //
3508   // do we need to tie the SDD left/right VDrift corrections
3509   for (int imd=0;imd<fNModules;imd++) {
3510     AliITSAlignMille2Module* mod = GetMilleModule(imd);
3511     if (mod->IsSDD() && mod->IsVDriftLRSame()) TieSDDVDriftsLR(mod);
3512   }
3513   //
3514 }
3515
3516 //________________________________________________________________________________________________________
3517 void AliITSAlignMille2::ApplyPostConstraints()
3518 {
3519   // apply constraints which can be imposed after the fit
3520   int nconstr = GetNConstraints();
3521   Bool_t convGlo      = kFALSE;
3522   // check if there is something to do
3523   int ntodo = 0;
3524   for (int i=0;i<nconstr;i++) {
3525     AliITSAlignMille2Constraint* cstr = GetConstraint(i);
3526     if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
3527     if (cstr->GetRemainingPattern() == 0) continue;
3528     ntodo++;
3529   }
3530   if (!ntodo) return;
3531   //
3532   if (!fUseGlobalDelta) { // need to convert to global params
3533     ConvertParamsToGlobal();
3534     convGlo = kTRUE;
3535   }
3536   //
3537   for (int i=0;i<nconstr;i++) {
3538     AliITSAlignMille2Constraint* cstr = GetConstraint(i);
3539     if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
3540     //
3541     int imd = cstr->GetModuleID();
3542     //
3543     if (imd>=0) {
3544       AliITSAlignMille2Module* mod = GetMilleModule(imd);
3545       if (mod->IsNotInConf()) continue;
3546       UInt_t pattern = 0;
3547       for (int ipar=mod->GetNParTot();ipar--;) {
3548         if (cstr->IsApplied(ipar))      continue;
3549         if (!cstr->IncludesParam(ipar)) continue;
3550         if (!mod->IsFreeDOF(ipar))      continue; // parameter is fixed, will not apply constraint
3551         pattern |= 0x1<<ipar;
3552         cstr->SetApplied(ipar);
3553       }
3554       if (pattern) PostConstrainModuleSubUnits(cstr->GetType(),cstr->GetModuleID(),cstr->GetValue(),pattern);
3555       //
3556     }
3557     else if (PseudoParentsAllowed()) {
3558       UInt_t pattern = (UInt_t)cstr->GetRemainingPattern();
3559       PostConstrainOrphans(cstr->GetType(),cstr->GetValue(),pattern);
3560       cstr->SetApplied(-1);
3561     }
3562   }
3563   // if there was a conversion, rewind it
3564   if (convGlo) ConvertParamsToLocal();
3565   //
3566 }
3567
3568 //________________________________________________________________________________________________________
3569 void AliITSAlignMille2::ConstrainModuleSubUnits(Int_t idm, Double_t val, UInt_t pattern)
3570 {
3571   // require that sum of modifications for the childs of this module is = val, i.e.
3572   // the internal corrections moves the module as a whole by fixed value (0 by default).
3573   // pattern is the bit pattern for the parameters to constrain
3574   //
3575   //
3576   AliITSAlignMille2Module* mod = GetMilleModule(idm);
3577   //
3578   for (int ip=0;ip<kNParCh;ip++) {
3579     if ( !((pattern>>ip)&0x1) /*|| !parent->IsFreeDOF(ip)*/) continue;
3580     ResetLocalEquation();
3581     int nadd = 0;
3582     for (int ich=mod->GetNChildren();ich--;) {
3583       int idpar = ((AliITSAlignMille2Module*)mod->GetChild(ich))->GetParOffset(ip);
3584       if (idpar<0) continue;
3585       fGlobalDerivatives[idpar] = 1.0;
3586       nadd++;
3587     }
3588     //
3589     if (nadd>0) {
3590       AddConstraint(fGlobalDerivatives,val);
3591       AliInfo(Form("Constrained param %d for %d submodules of module #%d: %s",ip,nadd,idm,mod->GetName()));
3592     }
3593   }
3594   //
3595 }
3596
3597 //________________________________________________________________________________________________________
3598 void AliITSAlignMille2::ConstrainOrphans(Double_t val, UInt_t pattern)
3599 {
3600   // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3601   // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3602   // pattern is the bit pattern for the parameters to constrain
3603   //
3604   for (int ip=0;ip<kNParCh;ip++) {
3605     //
3606     if ( !((pattern>>ip)&0x1) ) continue;
3607     ResetLocalEquation();
3608     int nadd = 0;
3609     for (int imd=fNModules;imd--;) {
3610       AliITSAlignMille2Module* mod = GetMilleModule(imd);
3611       if (mod->IsNotInConf()) continue; // dummy module
3612       AliITSAlignMille2Module* par = mod->GetParent();
3613       while (par && par->IsNotInConf() ) par = par->GetParent(); // use only decalred parents
3614       if (par) continue; // this is not an orphan
3615       int idpar = mod->GetParOffset(ip);
3616       if (idpar<0) continue;
3617       fGlobalDerivatives[idpar] = 1.0;
3618       nadd++;
3619     }
3620     if (nadd>0) {
3621       AddConstraint(fGlobalDerivatives,val);
3622       AliInfo(Form("Constrained param %d for %d orphan modules",ip,nadd));
3623     }
3624   }
3625   //
3626   //
3627 }
3628
3629 //________________________________________________________________________________________________________
3630 void AliITSAlignMille2::PostConstrainModuleSubUnits(Int_t type,Int_t idm, Double_t val, UInt_t pattern)
3631 {
3632   // require that median or mean of the modifications for the childs of this module is = val, i.e.
3633   // the internal corrections moves the module as a whole by fixed value (0 by default)
3634   // module the outliers.
3635   // pattern is the bit pattern for the parameters to constrain
3636   // The difference between the mean and the median will be transfered to the parent
3637   //
3638   AliITSAlignMille2Module* parent = GetMilleModule(idm);
3639   int nc = parent->GetNChildren();
3640   //
3641   double *tmpArr = new double[nc];
3642   //
3643   for (int ip=0;ip<kNParCh;ip++) {
3644     int npc = 0;
3645     if ( !((pattern>>ip)&0x1) || !parent->IsFreeDOF(ip)) continue;
3646     // compute the mean and median of the deltas
3647     int nfree = 0;
3648     for (int ich=nc;ich--;) {
3649       AliITSAlignMille2Module* child = parent->GetChild(ich);
3650       //      if (!child->IsFreeDOF(ip)) continue;
3651       tmpArr[nfree++] = child->GetParVal(ip);
3652     }
3653     double median=0,mean=0;
3654     for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
3655       mean += tmpArr[ic0];
3656       for (int ic1=ic0+1;ic1<nfree;ic1++)
3657         if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
3658     }
3659     //
3660     int kmed = nfree/2;
3661     median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
3662     if (nfree>0) mean /= nfree;
3663     //
3664     double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
3665     //
3666     for (int ich=nc;ich--;) {
3667       AliITSAlignMille2Module* child = parent->GetChild(ich);
3668       //    if (!child->IsFreeDOF(ip)) continue;
3669       child->SetParVal(ip, child->GetParVal(ip) + shift);
3670       npc++;
3671     }
3672     //
3673     parent->SetParVal(ip, parent->GetParVal(ip) - shift);
3674     AliInfo(Form("%s constraint: added %+f shift to param[%d] of %d children of module %d: %s",
3675                  type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
3676                  ip,npc,idm,parent->GetName()));
3677   }
3678   delete[] tmpArr;
3679   //
3680   //
3681 }
3682
3683 //________________________________________________________________________________________________________
3684 void AliITSAlignMille2::PostConstrainOrphans(Int_t type,Double_t val, UInt_t pattern)
3685 {
3686   // require that median or mean of modifications for the supermodules which have no parents is = val, i.e.
3687   // the corrections moves the whole setup by fixed value (0 by default).
3688   // pattern is the bit pattern for the parameters to constrain
3689   //
3690   int nc = fNModules;
3691   //
3692   int norph = 0;
3693   for (int ich=nc;ich--;) {
3694     AliITSAlignMille2Module *par= GetMilleModule(ich)->GetParent();
3695     while (par && par->IsNotInConf()) par = par->GetParent(); // use only decalred parents
3696     if (!par) norph ++;
3697   }
3698   //
3699   if (!norph) return;
3700   double *tmpArr = new double[norph];
3701   for (int i=norph;i--;) tmpArr[i] = 0;
3702   //
3703   for (int ip=0;ip<kNParCh;ip++) {
3704     int npc = 0;
3705     if ( !((pattern>>ip)&0x1)) continue;
3706     // compute the mean and median of the deltas
3707     int nfree = 0;
3708     for (int ich=nc;ich--;) {
3709       AliITSAlignMille2Module* child = GetMilleModule(ich);
3710       if (child->IsNotInConf()) continue; // dummy module
3711       //      if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
3712       AliITSAlignMille2Module* par = child->GetParent();
3713       while (par && par->IsNotInConf()) par = par->GetParent(); // count only declared parents
3714       if (par) continue;
3715       tmpArr[nfree++] = child->GetParVal(ip);
3716     }
3717     double median=0,mean=0;
3718     for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
3719       mean += tmpArr[ic0];
3720       for (int ic1=ic0+1;ic1<nfree;ic1++)
3721         if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
3722     }
3723     //
3724     int kmed = nfree/2;
3725     median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
3726     if (nfree>0) mean /= nfree;
3727     //
3728     double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
3729     //
3730     for (int ich=nc;ich--;) {
3731       AliITSAlignMille2Module* child = GetMilleModule(ich);
3732       if (child->IsNotInConf()) continue; // dummy module
3733       //      if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
3734       AliITSAlignMille2Module* par = child->GetParent();
3735       while (par && par->IsNotInConf()) par = par->GetParent(); // count only declared parents
3736       if (par) continue;
3737       child->SetParVal(ip, child->GetParVal(ip) + shift);
3738       npc++;
3739     }
3740     //
3741     AliInfo(Form("%s constraint: added %+f shift to param[%d] of %d orphan modules",
3742                  type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
3743                  ip,npc));
3744   }
3745   delete[] tmpArr;
3746   //
3747 }
3748
3749 //________________________________________________________________________________________________________
3750 Bool_t AliITSAlignMille2::IsParModConstrained(const AliITSAlignMille2Module* mod,Int_t par, Bool_t &meanmed, Bool_t &gaussian) const
3751 {
3752   // check if par of the module participates in some constraint, and set the flag for their types
3753   meanmed = gaussian = kFALSE;
3754   //
3755   if ( mod->IsParConstrained(par) ) gaussian = kTRUE;     // direct constraint on this param
3756   //
3757   for (int icstr=GetNConstraints();icstr--;) {
3758     AliITSAlignMille2Constraint* cstr = GetConstraint(icstr);
3759     //
3760     if (!cstr->IncludesModPar(mod,par)) continue;
3761     if (cstr->GetType()==AliITSAlignMille2ConstrArray::kTypeGaussian) gaussian = kTRUE;
3762     else meanmed = kTRUE;
3763     //
3764     if (meanmed && gaussian) break; // no sense to check further
3765   }
3766   //
3767   return meanmed||gaussian;
3768 }
3769
3770 //________________________________________________________________________________________________________
3771 Bool_t AliITSAlignMille2::IsParModFamilyVaried(const AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3772 {
3773   // check if parameter par is varied for this module or its children up to the level depth
3774   if (depth<0) return kFALSE;
3775   if (mod->GetParOffset(par)>=0) return kTRUE;
3776   for (int icld=mod->GetNChildren();icld--;) {
3777     AliITSAlignMille2Module* child = mod->GetChild(icld);
3778     if (IsParModFamilyVaried(child, par, depth-1)) return kTRUE;
3779   }
3780   return kFALSE;
3781   //
3782 }
3783
3784 /*
3785 //________________________________________________________________________________________________________
3786 Bool_t AliITSAlignMille2::IsParFamilyFree(AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3787 {
3788   // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
3789   if (depth<0) return kTRUE;
3790   for (int icld=mod->GetNChildren();icld--;) {
3791     AliITSAlignMille2Module* child = mod->GetChild(icld);
3792     //if (child->GetParOffset(par)<0) continue;                  // fixed
3793     Bool_t cstMM=kFALSE,cstGS=kFALSE;
3794     // does this child have gaussian constraint ?
3795     if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
3796     // check its children
3797     if (!IsParFamilyFree(child,par,depth-1)) return kTRUE;
3798   }
3799   return kFALSE;
3800   //
3801 }
3802 */
3803
3804 //________________________________________________________________________________________________________
3805 Bool_t AliITSAlignMille2::IsParFamilyFree(const AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3806 {
3807   // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
3808   if (depth<0) return kFALSE;
3809   for (int icld=mod->GetNChildren();icld--;) {
3810     AliITSAlignMille2Module* child = mod->GetChild(icld);
3811     //if (child->GetParOffset(par)<0) continue;                  // fixed
3812     Bool_t cstMM=kFALSE,cstGS=kFALSE;
3813     // does this child have gaussian constraint ?
3814     if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
3815     // check its children
3816     if (IsParFamilyFree(child,par,depth-1)) return kTRUE;
3817   }
3818   return kFALSE;
3819   //
3820 }
3821
3822 //________________________________________________________________________________________________________
3823 Double_t AliITSAlignMille2::GetTDriftSDD() const
3824 {
3825   // obtain drift time corrected for t0
3826   double t = fCluster.GetDriftTime();
3827   return t - fDriftTime0[ fCluster.GetUniqueID()-1 ];
3828 }
3829
3830 //________________________________________________________________________________________________________
3831 Double_t AliITSAlignMille2::GetVDriftSDD() const
3832 {
3833   // obtain corrected drift speed
3834   return fDriftSpeed[ fCluster.GetUniqueID()-1 ];
3835 }
3836
3837 //________________________________________________________________________________________________________
3838 Bool_t AliITSAlignMille2::FixedOrphans() const
3839 {
3840   // are there fixed modules with no parent (normally in such a case
3841   // the constraints on the orphans should not be applied
3842   if (!IsConfigured()) {
3843     AliInfo("Still not configured");
3844     return kFALSE;
3845   }
3846   for (int i=0;i<fNModules;i++) {
3847     AliITSAlignMille2Module* md = GetMilleModule(i);
3848     if (md->IsNotInConf()) continue;
3849     if (md->GetParent()==0 && md->GetNParFree()==0) return kTRUE;
3850   }
3851   return kFALSE;
3852 }
3853
3854 //________________________________________________________________________________________________________
3855 void AliITSAlignMille2::ConvertParamsToGlobal()
3856 {
3857   // convert params in local frame to global one
3858   double pars[AliITSAlignMille2Module::kMaxParGeom];
3859   for (int imd=fNModules;imd--;) {
3860     AliITSAlignMille2Module* mod = GetMilleModule(imd);
3861     if (mod->GeomParamsGlobal()) continue;
3862     mod->GetGeomParamsGlo(pars);
3863     mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
3864     mod->SetGeomParamsGlobal(kTRUE);
3865   }
3866 }
3867
3868 //________________________________________________________________________________________________________
3869 void AliITSAlignMille2::ConvertParamsToLocal()
3870 {
3871   // convert params in global frame to local one
3872   double pars[AliITSAlignMille2Module::kMaxParGeom];
3873   for (int imd=fNModules;imd--;) {
3874     AliITSAlignMille2Module* mod = GetMilleModule(imd);
3875     if (!mod->GeomParamsGlobal()) continue;
3876     mod->GetGeomParamsLoc(pars);
3877     mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
3878     mod->SetGeomParamsGlobal(kFALSE);
3879   }
3880 }
3881
3882 //________________________________________________________________________________________________________
3883 void AliITSAlignMille2::SetBField(Double_t b)
3884 {
3885   // set Bz value
3886   if (IsZero(b,1e-5)) {
3887     fBField = 0.0;
3888     fBOn = kFALSE;
3889     fNLocal = 4;
3890   }
3891   else {
3892     fBField = b;
3893     fBOn = kTRUE;
3894     fNLocal = 5; // helices
3895   }
3896 }
3897
3898 //________________________________________________________________________________________________________
3899 Int_t AliITSAlignMille2::ProcessUserInfo(TList* userInfo)
3900 {
3901   // extract calibration information used for TrackPointArray creation from run info
3902   //
3903   if (!userInfo) { AliInfo("No UserInfo is provided"); return 0;}
3904   //
3905   TMap *cdbMap=0;
3906   TList* cdbList=0;
3907   TObjString *objStr,*objStr1,*keyStr;
3908   TString cdbStr;
3909   AliCDBManager* man = AliCDBManager::Instance();
3910   man->SetCacheFlag(kFALSE);
3911   //
3912   int run = userInfo->GetUniqueID();
3913   if (run>0) SetRunID(run);
3914   AliInfo(Form("UserInfo corresponds to run#%d",run));
3915   cdbMap  = (TMap*)userInfo->FindObject("cdbMap");
3916   const TMap *curMap = man->GetStorageMap();
3917   if (!cdbMap) {AliInfo("No CDB Map found in UserInfo");}
3918   else {
3919     if ((objStr=(TObjString*)cdbMap->GetValue("default"))) { // first set default CDB path
3920       if ((objStr1=(TObjString*)curMap->GetValue("default")) && objStr1->GetUniqueID()) {
3921         AliInfo(Form("OCDB default path from UserInfo: %s is overriden by user setting %s",objStr->GetName(),objStr1->GetName()));
3922       }
3923       else {
3924         cdbStr = objStr->GetString();
3925         man->UnsetDefaultStorage();
3926         if (man->GetRaw()) man->SetRaw(kFALSE);
3927         if (cdbStr.BeginsWith("raw://")) cdbStr = "raw://";
3928         AliInfo(Form("Default CDB Storage from UserInfo: %s",cdbStr.Data()));
3929         man->SetDefaultStorage( cdbStr.Data() ); // this may be overriden later by configuration file
3930       }
3931     }
3932     if (man->GetRaw() && run>0) man->SetRun(run);
3933     //
3934     // set specific paths relevant for alignment
3935     TIter itMap(cdbMap);
3936     while( (keyStr=(TObjString*)itMap.Next()) ) {
3937       TString keyS = keyStr->GetString();
3938       if ( keyS == "default" ) continue;
3939       //
3940       TObjString* curPath = (TObjString*)curMap->GetValue(keyStr->GetName());
3941       if (curPath && curPath->GetUniqueID()) {
3942         AliInfo(Form("Storage for %s from UserInfo\n is overriden by user setting %s",keyS.Data(),curPath->GetName()));
3943         continue;
3944       }
3945       man->SetSpecificStorage( keyS.Data(), cdbMap->GetValue(keyS)->GetName() );
3946     }
3947   }
3948   //
3949   cdbList = (TList*)userInfo->FindObject("cdbList");
3950   if (!cdbList) {AliInfo("No CDB List found in UserInfo");}
3951   else {
3952     // Objects used for TrackPointArray production
3953     GetPathFromUserInfo(cdbList,"GRP/Geometry/Data",fIniGeomPath ,kSameInitGeomBit);
3954     GetPathFromUserInfo(cdbList,"ITS/Align/Data"   ,fIniDeltaPath,kSameInitDeltasBit);
3955     GetPathFromUserInfo(cdbList,"ITS/Calib/RespSDD",fIniSDDRespPath,kSameInitSDDRespBit);
3956     GetPathFromUserInfo(cdbList,"ITS/Calib/DriftSpeedSDD",fIniSDDVDriftPath,kSameInitSDDVDriftBit);
3957     GetPathFromUserInfo(cdbList,"ITS/Calib/MapsTimeSDD",fIniSDDCorrMapPath,kSameInitSDDCorrMapBit);
3958     GetPathFromUserInfo(cdbList,"GRP/Calib/MeanVertexSPD",fDiamondPath,kSameDiamondBit);
3959   }
3960   //
3961   TList *bzlst = (TList*)userInfo->FindObject("BzkGauss");
3962   if (bzlst && bzlst->At(0)) {
3963     objStr = (TObjString*)bzlst->At(0);
3964     SetBField( objStr->GetString().Atof() );
3965     AliInfo(Form("Magnetic field from UserInfo: %+.2e",GetBField()));
3966   }
3967   return 0;
3968 }
3969
3970 //________________________________________________________________________________________________________
3971 Int_t AliITSAlignMille2::GetPathFromUserInfo(TList* cdbList,const char* calib,TString& path, Int_t useBit)
3972 {
3973   // extract the path for specific CDB path from user info. If it is the same as already loaded, set corresponing bit
3974   TIter itList(cdbList);
3975   if (useBit>=0) ResetBit(useBit);
3976   TObjString* objStr;
3977   while( (objStr=(TObjString*)itList.Next()) )
3978     if (objStr->GetString().Contains(calib)) {
3979       TString newpath = objStr->GetString();
3980       AliInfo(Form("Found path in UserInfo: %s",newpath.Data()));
3981       if ( useBit>=0 && (fUserProvided&useBit) ) {
3982         AliInfo(Form("Will use the one provided in config: %s",path.Data()));
3983         SetBit(useBit);
3984       }
3985       else if ( useBit>=0 && (newpath == path) ) {
3986         AliInfo(Form("Path %s is the same as already loaded",path.Data()));
3987         SetBit(useBit);
3988       }
3989       else path = newpath;
3990       //
3991       return 0;
3992     }
3993   AliInfo(Form("Did not find path for %s in UserInfo",calib));
3994   path = "";
3995   return -1;
3996 }
3997
3998 //________________________________________________________________________________________________________
3999 Int_t AliITSAlignMille2::LoadSDDResponse(TString& path, AliITSresponseSDD *&resp)
4000 {
4001   // load SDD response
4002   if (path.IsNull()) return 0;
4003   AliInfo(Form("Loading SDD response from %s",path.Data()));
4004   //
4005   AliCDBEntry *entry = 0;
4006   delete resp;
4007   resp = 0;
4008   //
4009   while(1) {
4010     if (path.BeginsWith("path: ")) { // must load from OCDB
4011       entry = GetCDBEntry(path.Data());
4012       if (!entry) break;
4013       resp = (AliITSresponseSDD*) entry->GetObject();
4014       entry->SetObject(NULL);
4015       entry->SetOwner(kTRUE);
4016       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4017       //      delete cdbId;
4018       //      delete entry;
4019       break;
4020     }
4021     //
4022     if (gSystem->AccessPathName(path.Data())) break;
4023     TFile* precf = TFile::Open(path.Data());
4024     if (precf->FindKey("AliITSresponseSDD")) resp = (AliITSresponseSDD*)precf->Get("AliITSresponseSDD");
4025     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4026       resp = (AliITSresponseSDD*) entry->GetObject();
4027       if (resp && resp->InheritsFrom(AliITSresponseSDD::Class())) entry->SetObject(NULL);
4028       else resp = 0;
4029       entry->SetObject(NULL);
4030       entry->SetOwner(kTRUE);
4031       delete entry;
4032     }
4033     //
4034     precf->Close();
4035     delete precf;
4036     break;
4037   }
4038   //
4039   if (!resp) {AliError(Form("Failed to load SDD response from %s",path.Data())); return -1;}
4040   return 0;
4041 }
4042
4043 //________________________________________________________________________________________________________
4044 Int_t AliITSAlignMille2::LoadSDDVDrift(TString& path, TObjArray *&arr)
4045 {
4046   // load VDrift object
4047   if (path.IsNull()) return 0;
4048   AliInfo(Form("Loading SDD VDrift from %s",path.Data()));
4049   //
4050   AliCDBEntry *entry = 0;
4051   delete arr;
4052   arr = 0;
4053   while(1) {
4054     if (path.BeginsWith("path: ")) { // must load from OCDB
4055       entry = GetCDBEntry(path.Data());
4056       if (!entry) break;
4057       arr = (TObjArray*) entry->GetObject();
4058       entry->SetObject(NULL);
4059       entry->SetOwner(kTRUE);
4060       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4061       //      delete cdbId;
4062       //      delete entry;
4063       break;
4064     }
4065     //
4066     if (gSystem->AccessPathName(path.Data())) break;
4067     TFile* precf = TFile::Open(path.Data());
4068     if (precf->FindKey("TObjArray")) arr = (TObjArray*)precf->Get("TObjArray");
4069     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4070       arr = (TObjArray*) entry->GetObject();
4071       if (arr && arr->InheritsFrom(TObjArray::Class())) entry->SetObject(NULL);
4072       else arr = 0;
4073       entry->SetObject(NULL);
4074       entry->SetOwner(kTRUE);
4075       delete entry;
4076     }
4077     //
4078     precf->Close();
4079     delete precf;
4080     break;
4081   }
4082   //
4083   if (!arr) {AliError(Form("Failed to load SDD vdrift from %s",path.Data())); return -1;}
4084   arr->SetOwner(kTRUE);
4085   return 0;
4086 }
4087
4088 //________________________________________________________________________________________________________
4089 Int_t AliITSAlignMille2::LoadSDDCorrMap(TString& path, AliITSCorrectSDDPoints *&map)
4090 {
4091   // Load SDD correction map
4092   //
4093   if (path.IsNull()) return 0;
4094   AliInfo(Form("Loading SDD Correction Maps from %s",path.Data()));
4095   //
4096   AliCDBEntry *entry = 0;
4097   delete map;
4098   map = 0;
4099   TObjArray* arr = 0;
4100   while(1) {
4101     if (path.BeginsWith("path: ")) { // must load from OCDB
4102       entry = GetCDBEntry(path.Data());
4103       if (!entry) break;
4104       arr = (TObjArray*) entry->GetObject();
4105       entry->SetObject(NULL);
4106       entry->SetOwner(kTRUE);
4107       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4108       //      delete cdbId;
4109       //      delete entry;
4110       break;
4111     }
4112     //
4113     if (gSystem->AccessPathName(path.Data())) break;
4114     TFile* precf = TFile::Open(path.Data());
4115     if (precf->FindKey("TObjArray")) arr = (TObjArray*)precf->Get("TObjArray");
4116     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4117       arr = (TObjArray*) entry->GetObject();
4118       if (arr && arr->InheritsFrom(TObjArray::Class())) entry->SetObject(NULL);
4119       else arr = 0;
4120       entry->SetObject(NULL);
4121       entry->SetOwner(kTRUE);
4122       delete entry;
4123     }
4124     //
4125     precf->Close();
4126     delete precf;
4127     break;
4128   }
4129   //
4130   if (!arr) {AliError(Form("Failed to load SDD Correction Map from %s",path.Data())); return -1;}
4131   arr->SetOwner(kTRUE);
4132   map = new AliITSCorrectSDDPoints(arr);
4133
4134   return 0;
4135 }
4136
4137 //________________________________________________________________________________________________________
4138 Int_t AliITSAlignMille2::LoadDiamond(TString& path)
4139 {
4140   // load vertex constraint
4141   if (path.IsNull()) return 0;
4142   AliInfo(Form("Loading Diamond Constraint from %s",path.Data()));
4143   //
4144   AliCDBEntry *entry = 0;
4145   AliESDVertex *vtx = 0;
4146   while(1) {
4147     if (path.BeginsWith("path: ")) { // must load from OCDB
4148       entry = GetCDBEntry(path.Data());
4149       if (!entry) break;
4150       vtx = (AliESDVertex*) entry->GetObject();
4151       entry->SetObject(NULL);
4152       entry->SetOwner(kTRUE);
4153       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4154       //      delete cdbId;
4155       //      delete entry;
4156       break;
4157     }
4158     //
4159     if (gSystem->AccessPathName(path.Data())) break;
4160     TFile* precf = TFile::Open(path.Data());
4161     if (precf->FindKey("AliESDVertex")) vtx = (AliESDVertex*)precf->Get("AliESDVertex");
4162     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4163       vtx = (AliESDVertex*) entry->GetObject();
4164       if (vtx && vtx->InheritsFrom(AliESDVertex::Class())) entry->SetObject(NULL);
4165       else vtx = 0;
4166       entry->SetObject(NULL);
4167       entry->SetOwner(kTRUE);
4168       delete entry;
4169     }
4170     //
4171     precf->Close();
4172     delete precf;
4173     break;
4174   }
4175   //
4176   if (!vtx) {AliError(Form("Failed to load Diamond constraint from %s",path.Data())); return -1;}
4177   //
4178   double vtxXYZ[3];
4179   vtx->GetXYZ(vtxXYZ);
4180   for (int i=3;i--;) vtxXYZ[i] -= fCorrDiamond[i];
4181   vtx->SetXYZ(vtxXYZ);
4182   SetVertexConstraint(vtx);
4183   AliInfo("Will use following Diamond Constraint (errors inverted):");
4184   fDiamondI.Print("");
4185   delete vtx;
4186   return 0;
4187 }
4188
4189 //________________________________________________________________________________________________________
4190 Int_t AliITSAlignMille2::LoadDeltas(TString& path, TClonesArray *&arr)
4191 {
4192   // load ITS geom deltas
4193   if (path.IsNull()) return 0;
4194   AliInfo(Form("Loading Alignment Deltas from %s",path.Data()));
4195   //
4196   AliCDBEntry *entry = 0;
4197   delete arr;
4198   arr = 0;
4199   while(1) {
4200     if (path.BeginsWith("path: ")) { // must load from OCDB
4201       entry = GetCDBEntry(path.Data());
4202       if (!entry) break;
4203       arr = (TClonesArray*) entry->GetObject();
4204       entry->SetObject(NULL);
4205       entry->SetOwner(kTRUE);
4206       //      AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4207       //      delete cdbId;
4208       //      delete entry;
4209       break;
4210     }
4211     //
4212     if (gSystem->AccessPathName(path.Data())) break;
4213     TFile* precf = TFile::Open(path.Data());
4214     if (precf->FindKey("ITSAlignObjs")) arr = (TClonesArray*)precf->Get("ITSAlignObjs");
4215     else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4216       arr = (TClonesArray*) entry->GetObject();
4217       if (arr && arr->InheritsFrom(TClonesArray::Class())) entry->SetObject(NULL);
4218       else arr = 0;
4219       entry->SetObject(NULL);
4220       entry->SetOwner(kTRUE);
4221       delete entry;
4222     }
4223     precf->Close();
4224     delete precf;
4225     break;
4226   }
4227   //
4228   if (!arr) {AliError(Form("Failed to load Deltas from %s",path.Data())); return -1;}
4229   //
4230   return 0;
4231 }
4232
4233 //________________________________________________________________________________________________________
4234 Int_t AliITSAlignMille2::CacheMatricesCurr()
4235 {
4236   // build arrays for the fast access to sensor matrices from their sensor ID
4237   //
4238   TGeoHMatrix mdel;
4239   AliInfo("Building sensors current matrices cache");
4240   //
4241   fCacheMatrixCurr.Delete();
4242   for (int idx=0;idx<=kMaxITSSensID;idx++) {
4243     int volID = AliITSAlignMille2Module::GetVolumeIDFromIndex(idx);
4244     TGeoHMatrix *mcurr = new TGeoHMatrix();
4245     AliITSAlignMille2Module::SensVolMatrix(volID, mcurr);
4246     fCacheMatrixCurr.AddAtAndExpand(mcurr,idx);
4247     //
4248   }
4249   //
4250   TGeoHMatrix *mcurr = new TGeoHMatrix();
4251   fCacheMatrixCurr.AddAtAndExpand(mcurr,kVtxSensID); // special unit matrix for diamond constraint
4252   //
4253   fCacheMatrixCurr.SetOwner(kTRUE);
4254   return 0;
4255 }
4256
4257 //________________________________________________________________________________________________________
4258 Int_t AliITSAlignMille2::CacheMatricesOrig()
4259 {
4260   // build arrays for the fast access to sensor original matrices (used for production)
4261   //
4262   TGeoHMatrix mdel;
4263   AliInfo("Building sensors original matrices cache");
4264   //
4265   /*if (fIniGeomPath!=fGeometryPath)*/ if (LoadGeometry(fIniGeomPath)) {AliInfo("Failed to re-load ideal geometry");exit(1);}
4266   //
4267   fCacheMatrixOrig.Delete();
4268   if (!fIniDeltaPath.IsNull()) {
4269     TClonesArray* prealSav = fPrealignment;
4270     fPrealignment = 0;
4271     if (LoadDeltas(fIniDeltaPath,fPrealignment) || ApplyToGeometry())
4272       { AliError("Failed to load/apply initial deltas used to produce points"); return -1;}
4273     delete fPrealignment;
4274     fPrealignment = prealSav;
4275   }
4276   //
4277   for (int idx=0;idx<=kMaxITSSensID;idx++) {
4278     int volID = AliITSAlignMille2Module::GetVolumeIDFromIndex(idx);
4279     TGeoHMatrix *morig = new TGeoHMatrix();
4280     AliITSAlignMille2Module::SensVolMatrix(volID,morig);
4281     fCacheMatrixOrig.AddAtAndExpand(morig,idx);
4282   }
4283   //
4284   if (fConvertPreDeltas) {
4285     // in order to convert deltas from old to new geometry we need the final matrices for all alignable objects
4286     int nmat = fGeoManager->GetNAlignable();
4287     fConvAlgMatOld.Delete();
4288     int nmatSel = 0;
4289     for (int i=0;i<nmat;i++) {
4290       TString nm = fGeoManager->GetAlignableEntry(i)->GetName();
4291       if (!nm.BeginsWith("ITS")) continue;
4292       TGeoHMatrix *mo = new TGeoHMatrix();
4293       (*mo) = *(AliGeomManager::GetMatrix(nm));
4294       fConvAlgMatOld.AddAtAndExpand(mo,nmatSel++);
4295       mo->SetTitle(nm);
4296       mo->SetName(nm);
4297     }
4298     ConvSortHierarchically(fConvAlgMatOld);
4299   }
4300   //
4301   TGeoHMatrix *mcurr = new TGeoHMatrix();
4302   fCacheMatrixOrig.AddAtAndExpand(mcurr,kVtxSensID); // special unit matrix for diamond constraint
4303   //
4304   fCacheMatrixOrig.SetOwner(kTRUE);
4305
4306   fUsePreAlignment = 0;
4307   LoadGeometry(fGeometryPath);   // reload target geometry
4308   //
4309   return 0;
4310 }
4311
4312 //________________________________________________________________________________________________________
4313 void AliITSAlignMille2::RemoveHelixFitConstraint()
4314 {
4315   // suppress constraint
4316   fConstrCharge = 0;
4317   fConstrPT = fConstrPTErr = -1;
4318 }
4319
4320 //________________________________________________________________________________________________________
4321 void AliITSAlignMille2::ConstrainHelixFitPT(Int_t q,Double_t pt,Double_t pterr)
4322 {
4323   // constrain q and pT of the helical fit of the track (should be set before process.track)
4324   //
4325   fConstrCharge = q==0 ? q:TMath::Sign(1,q);
4326   fConstrPT = pt;
4327   fConstrPTErr = pterr;
4328 }
4329
4330 //________________________________________________________________________________________________________
4331 void AliITSAlignMille2::ConstrainHelixFitCurv(Int_t q,Double_t crv,Double_t crverr)
4332 {
4333   // constrain charge and curvature of the helical fit of the track (should be set before process.track)
4334   //
4335   const double kCQConv = 0.299792458e-3;// R = PT/Bz/fgkCQConv with GeV,kGauss,cm
4336
4337   fConstrCharge = q==0 ? q:TMath::Sign(1,q);
4338   if (crv<0 || IsZero(crv)) {
4339     fConstrPT    = -1;
4340     fConstrPTErr = -1;
4341   }
4342   else {
4343     fConstrPT    = TMath::Abs(1./crv*fBField*kCQConv);
4344     fConstrPTErr = crverr>1e-10 ? TMath::Abs(fConstrPT/crv*crverr) : 0.;
4345   }
4346 }
4347
4348 //________________________________________________________________________________________________________
4349 TClonesArray* AliITSAlignMille2::CreateDeltas()
4350 {
4351   // Create \Deltas for every explicitly or implicitly (via non-alignable volumes) varied
4352   // or prealigned module.
4353   // If the module has inded J in the hierarchy of alignable volumes (0 - the top, most
4354   // coarse level), then its Delta is expressed via MP2 \deltas (in global frame) and
4355   // prealignment \DeltaP's as:
4356   // \Delta_J = Y X Y^-1
4357   // where X = \delta_J * \DeltaP_J
4358   // Y = Prod_{K=0,J-1} \delta_K
4359   // Note that \delta_L accounts not only for its own correction but also of all non-alignable
4360   // modules in the hierarchy chain from L up to the closest alignable:
4361   // while (parent && !parent->IsAlignable()) {
4362   //   \delta_L->MultiplyLeft( \delta_parent );
4363   //   parent = parent->GetParent();
4364   // }
4365   //
4366   Bool_t convLoc = kFALSE;
4367   if (!GetUseGlobalDelta()) {
4368     ConvertParamsToGlobal();
4369     convLoc = kTRUE;
4370   }
4371   //
4372   AliAlignObjParams tempAlignObj;
4373   TGeoHMatrix tempMatX,tempMatY,tempMat1;
4374   //
4375   TClonesArray *array = new TClonesArray("AliAlignObjParams",10);
4376   TClonesArray &alobj = *array;
4377   int idx = 0;
4378   //
4379   TGeoManager* geoManager = AliGeomManager::GetGeometry();
4380   int nalgtot = geoManager->GetNAlignable();
4381   //
4382   for (int ialg=0;ialg<nalgtot;ialg++) {             // loop over all alignable entries
4383     //
4384     const char* algname = geoManager->GetAlignableEntry(ialg)->GetName();
4385     //
4386     AliITSAlignMille2Module* md     = GetMilleModuleBySymName(algname); // explicitly varied?
4387     AliITSAlignMille2Module* parent = md ? md->GetParent(): GetMilleModuleIfContained(algname);
4388     if (md && parent) {
4389       TString mdName = md->GetName();
4390       TString prName = parent->GetName();
4391       // SPD Sector -> Layer parentship is fake, need special treatment
4392       if ( mdName.CountChar('/')==2 && mdName.BeginsWith("ITS/SPD") && // SPD sector
4393            prName.CountChar('/')==1 && mdName.BeginsWith("ITS/SPD") )  // SPD Layer
4394         parent = parent->GetParent();//: GetMilleModuleIfContained(prName.Data());
4395     }
4396     //
4397     AliAlignObjParams*       preob  = GetPrealignedObject(algname);  // was it prealigned ?
4398     //
4399     if (!preob && !md && (!parent || parent->IsAlignable())) continue; // noting to do
4400     //
4401     // create matrix X (see comment) ------------------------------------------------->>>
4402     // start from unity matrix
4403     tempMatX.Clear();
4404     if (preob) {   // account prealigngment
4405       preob->GetMatrix(tempMat1);
4406       tempMatX.MultiplyLeft(&tempMat1);
4407     }
4408     //
4409     if (md) {
4410       tempAlignObj.SetTranslation( md->GetParVal(0),md->GetParVal(1),md->GetParVal(2));
4411       tempAlignObj.SetRotation(    md->GetParVal(3),md->GetParVal(4),md->GetParVal(5));
4412       tempAlignObj.GetMatrix(tempMat1);
4413       tempMatX.MultiplyLeft(&tempMat1);  // acount correction to varied module
4414     }
4415     //
4416     // the corrections to all non-alignable modules from current on
4417     // till first alignable should add up to its matrix
4418     while (parent && !parent->IsAlignable()) {
4419       tempAlignObj.SetTranslation( parent->GetParVal(0),parent->GetParVal(1),parent->GetParVal(2));
4420       tempAlignObj.SetRotation(    parent->GetParVal(3),parent->GetParVal(4),parent->GetParVal(5));
4421       tempAlignObj.GetMatrix(tempMat1);
4422       tempMatX.MultiplyLeft(&tempMat1);  // add matrix of non-alignable module
4423       parent = parent->GetParent();
4424     }
4425     // create matrix X (see comment) ------------------------------------------------<<<
4426     //
4427     // create matrix Y (see comment) ------------------------------------------------>>>
4428     // start from unity matrix
4429     tempMatY.Clear();
4430     while ( parent ) {
4431       tempAlignObj.SetTranslation( parent->GetParVal(0),parent->GetParVal(1),parent->GetParVal(2));
4432       tempAlignObj.SetRotation(    parent->GetParVal(3),parent->GetParVal(4),parent->GetParVal(5));
4433       tempAlignObj.GetMatrix(tempMat1);
4434       tempMatY.MultiplyLeft(&tempMat1);
4435       parent = parent->GetParent();
4436     }
4437     // create matrix Y (see comment) ------------------------------------------------<<<
4438     //
4439     tempMatX.MultiplyLeft(&tempMatY);
4440     tempMatX.Multiply(&tempMatY.Inverse());
4441     //
4442     if (tempMatX.IsIdentity()) continue; // do not store dummy matrices
4443     UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(algname);
4444     new(alobj[idx++]) AliAlignObjParams(algname,vid,tempMatX,kTRUE);
4445     //
4446   }
4447   //
4448   if (convLoc) ConvertParamsToLocal();
4449   //
4450   return array;
4451   //
4452 }
4453
4454 //_______________________________________________________________________________________
4455 AliITSresponseSDD* AliITSAlignMille2::CreateSDDResponse()
4456 {
4457   // create object with SDD repsonse (t0 and vdrift corrections) accounting for
4458   // eventual precalibration
4459   //
4460   // if there was a precalibration provided, copy it to new arrray
4461   AliITSresponseSDD *precal = GetSDDPrecalResp();
4462   if (!precal && fIniVDriftSDD)       precal = GetSDDInitResp();    // InitResp is used only when IniVDrift is provided
4463   Bool_t isPreCalMult = precal&&precal->IsVDCorrMult() ? kTRUE : kFALSE;
4464   AliITSresponseSDD *calibSDD = new AliITSresponseSDD();
4465   calibSDD->SetVDCorrMult(fIsSDDVDriftMult);
4466   //
4467   // copy initial values to the new object
4468   if (precal) {
4469     calibSDD->SetTimeOffset(precal->GetTimeOffset());
4470     calibSDD->SetADC2keV(precal->GetADC2keV());
4471     calibSDD->SetChargevsTime(precal->GetChargevsTime());
4472     for (int ind=kSDDoffsID;ind<kSDDoffsID+kNSDDmod;ind++) {
4473       calibSDD->SetModuleTimeZero(ind, precal->GetTimeZero(ind));
4474       calibSDD->SetDeltaVDrift(ind, precal->GetDeltaVDrift(ind,kFALSE),kFALSE); // left
4475       calibSDD->SetDeltaVDrift(ind, precal->GetDeltaVDrift(ind,kTRUE ),kTRUE);  // right
4476       calibSDD->SetADCtokeV(ind,precal->GetADCtokeV(ind));
4477     }
4478   }
4479   else for (int ind=kSDDoffsID;ind<kSDDoffsID+kNSDDmod;ind++) calibSDD->SetModuleTimeZero(ind,0);
4480   //
4481   Bool_t save = kFALSE;
4482   for (int imd=GetNModules();imd--;) {
4483     AliITSAlignMille2Module* md = GetMilleModule(imd);
4484     if (!md->IsSDD()) continue;
4485     if (md->IsFreeDOF(AliITSAlignMille2Module::kDOFT0)  ||
4486         md->IsFreeDOF(AliITSAlignMille2Module::kDOFDVL) ||
4487         md->IsFreeDOF(AliITSAlignMille2Module::kDOFDVR)) save = kTRUE;
4488         //
4489     for (int is=0;is<md->GetNSensitiveVolumes();is++) {
4490       int ind  = md->GetSensVolIndex(is);
4491       float t0  = calibSDD->GetTimeZero(ind)      + md->GetParVal(AliITSAlignMille2Module::kDOFT0);
4492       double dvL = md->GetParVal(AliITSAlignMille2Module::kDOFDVL);
4493       double dvR = md->GetParVal(AliITSAlignMille2Module::kDOFDVR);
4494       if (!calibSDD->IsVDCorrMult()) { // save as additive correction
4495         dvL *= 1e4;
4496         dvR *= 1e4;
4497         //
4498         double conv = 1;
4499         if (isPreCalMult) conv = 6.4; // convert multiplicative precal correction to additive
4500         dvL += calibSDD->GetDeltaVDrift(ind,kFALSE)*conv;
4501         dvR += calibSDD->GetDeltaVDrift(ind,kTRUE)*conv;
4502       }
4503       else { // save as multipicative correction
4504         double conv = 1;
4505         if (!isPreCalMult) conv = 1./6.4; // convert additive precal correction to multiplicative
4506         dvL += calibSDD->GetDeltaVDrift(ind,kFALSE)*conv;
4507         dvR += calibSDD->GetDeltaVDrift(ind,kTRUE)*conv;
4508       }
4509       //
4510       calibSDD->SetModuleTimeZero(ind, t0);
4511       calibSDD->SetDeltaVDrift(ind, dvL, kFALSE); // left  side correction
4512       calibSDD->SetDeltaVDrift(ind, dvR, kTRUE); // right side correction
4513     }
4514   }
4515   //
4516   if (!save) {
4517     AliInfo("No free parameters for SDD calibration, nothing to save");
4518     delete calibSDD;
4519     calibSDD = 0;
4520   }
4521   //
4522   return calibSDD;
4523 }
4524
4525 //_______________________________________________________________________________________
4526 Int_t AliITSAlignMille2::ReloadInitCalib(TList *userInfo)
4527 {
4528   // Use provided UserInfo to
4529   // load the initial calib parameters (geometry, SDD response...)
4530   // Can be used if set of data was processed with different calibration
4531   //
4532   if (!userInfo) {
4533     AliInfo("Reloading of the Calibration parameters was called with empty userInfo");
4534     return 1;
4535   }
4536   if (ProcessUserInfo(userInfo)) {
4537     AliInfo("Error in processing user info");
4538     userInfo->Print();
4539     exit(1);
4540   }
4541   return ReloadInitCalib();
4542 }
4543
4544 //_______________________________________________________________________________________
4545 Int_t AliITSAlignMille2::ReloadInitCalib()
4546 {
4547   // Load the initial calib parameters (geometry, SDD response...)
4548   // Can be used if set of data was processed with different calibration
4549   //
4550   // 1st cache original matrices
4551   if (!(TestBit(kSameInitDeltasBit) && TestBit(kSameInitGeomBit))) { // need to reload geometry
4552     //
4553     if (CacheMatricesOrig()) {
4554       AliInfo("Failed to cache new initial geometry");
4555       exit(1);
4556     }
4557     // RS : commented because we don't need to reload prealignment deltas, they are already loaded
4558     // then reload the prealignment geometry
4559     //    if (LoadDeltas(fPreDeltaPath,fPrealignment)) {
4560     //      AliInfo(Form("Failed to reload the prealigned geometry %s",fPreDeltaPath.Data()));
4561     //      exit(1);
4562     //    }
4563     //
4564     if (fPrealignment && ApplyToGeometry()) {
4565       AliInfo(Form("Failed re-apply prealigned geometry %s",fPreDeltaPath.Data()));
4566       exit(1);
4567     }
4568     //
4569     // usually no need to re-cache the prealignment geometry, it was not changed
4570     if (fCacheMatrixCurr.GetEntriesFast() != fCacheMatrixOrig.GetEntriesFast()) {
4571       //      CacheMatricesCurr();
4572       AliInfo(Form("Failed to cache the prealigned geometry %s",fPreDeltaPath.Data()));
4573       exit(1);
4574     }
4575   }
4576   else ResetBit(kSameInitDeltasBit);
4577   //
4578   // reload initial SDD response
4579   if (!TestBit(kSameInitSDDRespBit)) {
4580     if (LoadSDDResponse(fIniSDDRespPath, fIniRespSDD) ) {
4581       AliInfo(Form("Failed to load new SDD response %s",fIniSDDRespPath.Data()));
4582       exit(1);
4583     }
4584   }
4585   else ResetBit(kSameInitSDDRespBit);
4586   //
4587   // reload initial SDD vdrift
4588   if (!TestBit(kSameInitSDDVDriftBit)) {
4589     if (LoadSDDVDrift(fIniSDDVDriftPath, fIniVDriftSDD) ) {
4590       AliInfo(Form("Failed to load new SDD VDrift %s",fIniSDDVDriftPath.Data()));
4591       exit(1);
4592     }
4593   }
4594   else ResetBit(kSameInitSDDRespBit);
4595   //
4596   // reload SDD corr.map
4597   if (!TestBit(kSameInitSDDCorrMapBit)) {
4598     if (LoadSDDCorrMap(fIniSDDCorrMapPath, fIniCorrMapSDD) ) {
4599       AliInfo(Form("Failed to load new SDD Correction Map %s",fIniSDDCorrMapPath.Data()));
4600       exit(1);
4601     }
4602   }
4603   else ResetBit(kSameInitSDDRespBit);
4604   //
4605   // reload diamond info
4606   if (!TestBit(kSameDiamondBit)) {
4607     if (LoadDiamond(fDiamondPath) ) {
4608       AliInfo(Form("Failed to load new Diamond constraint %s",fDiamondPath.Data()));
4609       exit(1);
4610     }
4611   }
4612   else ResetBit(kSameInitSDDRespBit);
4613   //
4614   return 0;
4615 }
4616
4617 //_______________________________________________________________________________________
4618 void AliITSAlignMille2::JacobianPosGloLoc(int locid,double* jacobian)
4619 {
4620   // calculate the locid row of the jacobian for transformation of the local coordinate to global at current point
4621   TGeoHMatrix* mat = GetSensorCurrMatrixSID(fCurrentSensID);
4622   const Double_t dpar = 1e-2;
4623   double sav = fMeasLoc[locid];
4624   fMeasLoc[locid] += dpar;
4625   mat->LocalToMaster(fMeasLoc,jacobian);
4626   fMeasLoc[locid] = sav; // recover original value
4627   for (int i=3;i--;) jacobian[i] = (jacobian[i]-fMeasGlo[i])/dpar; // the transformation is linear!!!
4628 }
4629
4630 //_______________________________________________________________________________________
4631 void AliITSAlignMille2::TieSDDVDriftsLR(AliITSAlignMille2Module* mod)
4632 {
4633   // impose equality of Left/Right sides VDrift correction for SDD
4634   ResetLocalEquation();
4635   if ( (mod->IsFreeDOF(AliITSAlignMille2Module::kDOFDVL) + mod->IsFreeDOF(AliITSAlignMille2Module::kDOFDVR))==1) {
4636     AliError("Left/Right VDrift equality is requested for SDD module with only one side VDrift free");
4637     mod->Print();
4638     return;
4639   }
4640   if (mod->GetParOffset(AliITSAlignMille2Module::kDOFDVL)>=0) SetGlobalDerivative(mod->GetParOffset(AliITSAlignMille2Module::kDOFDVL),  1.);
4641   if (mod->GetParOffset(AliITSAlignMille2Module::kDOFDVR)>=0) SetGlobalDerivative(mod->GetParOffset(AliITSAlignMille2Module::kDOFDVR), -1.);
4642   AddConstraint(fGlobalDerivatives, 0, 1e-12);
4643   //
4644 }
4645
4646 //_______________________________________________________________________________________
4647 void AliITSAlignMille2::ProcessSDDPointInfo(const AliTrackPoint* pnt,Int_t sID, Int_t pntID)
4648 {
4649   // extract the drift information from SDD track point
4650   //
4651   fDriftTime0[pntID] = fIniRespSDD ? fIniRespSDD->GetTimeZero(sID) : 0.;
4652   double tdif = pnt->GetDriftTime() - fDriftTime0[pntID];
4653   if (tdif<0) tdif = 1;
4654   //
4655   // VDrift extraction
4656   double vdrift=0,vdrift0=0;
4657   Bool_t sddSide = kFALSE;
4658   int sID0 = 2*(sID-kSDDoffsID);
4659   double zanode = -999;
4660   //
4661   if (fIniVDriftSDD) { // SDD VDrift object is provided, use the vdrift from it
4662     AliITSDriftSpeedArraySDD* drarr;
4663     double vdR,vdL,xlR,xlL;
4664     // sometimes xlocal on right side is negative due to the wrong calibration, need to test both hypothesis
4665     double xlabs = TMath::Abs(fMeasLoc[kX]);
4666     drarr  = (AliITSDriftSpeedArraySDD*)fIniVDriftSDD->At(sID0); // left side, xloc>0
4667     zanode = fSegmentationSDD->GetAnodeFromLocal(xlabs,fMeasLoc[kZ]);
4668     vdL    = drarr->GetDriftSpeed(0, zanode);
4669     if (fIniRespSDD) {
4670       double corr = fIniRespSDD->GetDeltaVDrift(sID, kFALSE);
4671       if (fIniRespSDD->IsVDCorrMult()) vdL *= (1+corr);
4672       else vdL += corr;
4673     }
4674     xlL    = (fSegmentationSDD->Dx() - vdL*tdif)*1e-4;
4675     //
4676     drarr  = (AliITSDriftSpeedArraySDD*)fIniVDriftSDD->At(sID0+1); // right side, xloc<0
4677     zanode = fSegmentationSDD->GetAnodeFromLocal(-xlabs,fMeasLoc[kZ]) - 256;
4678     vdR    = drarr->GetDriftSpeed(0, zanode);
4679     if (fIniRespSDD) {
4680       double corr = fIniRespSDD->GetDeltaVDrift(sID, kTRUE);
4681       if (fIniRespSDD->IsVDCorrMult()) vdR *= (1+corr);
4682       else vdR += corr;
4683     }
4684     xlR    = -(fSegmentationSDD->Dx() - vdR*tdif)*1e-4;
4685     //
4686     if (TMath::Abs(xlL-fMeasLoc[kX])<TMath::Abs(xlR-fMeasLoc[kX])) {
4687       sddSide = 0; // left side
4688       vdrift  = vdL*1e-4;
4689     }
4690     else {         // right side
4691       sddSide = 1;
4692       vdrift  = vdR*1e-4;
4693     }
4694     //
4695   }
4696   else { // try to determine the vdrift from the xloc
4697     vdrift = (fSegmentationSDD->Dx()*1e-4 - TMath::Abs(fMeasLoc[kX]))/tdif;
4698     sddSide = fMeasLoc[kX]<0; // 0 = left (xloc>0) ; 1 = right (xloc<1)
4699   }
4700   //
4701   if (fPreVDriftSDD) { // use imposed vdrift as a starting point
4702     zanode = fSegmentationSDD->GetAnodeFromLocal(0.5-sddSide,fMeasLoc[kZ]);
4703     if (sddSide) zanode -= 256;
4704     vdrift = ((AliITSDriftSpeedArraySDD*)fPreVDriftSDD->At(sID0+sddSide))->GetDriftSpeed(0, zanode)*1e-4;
4705   }
4706   //
4707   if (vdrift<0) vdrift = 0;
4708   vdrift0 = vdrift;
4709   // at this point we have vdrift and t0 used to create the original point.
4710   // see if precalibration was provided
4711   if (fPreRespSDD) {
4712     float t0Upd = fPreRespSDD->GetTimeZero(sID);
4713     double corr = fPreRespSDD->GetDeltaVDrift(sID, sddSide);
4714     if (fPreRespSDD->IsVDCorrMult()) vdrift *= 1+corr; // right side (xloc<0) may have different correction
4715     else                             vdrift += corr*1e-4;
4716     //
4717     // if IniRespSDD was used, it should be subtracted back, since it is accounted in the PreResp
4718     if (fIniVDriftSDD&&fIniRespSDD) {
4719       double corr1 = fIniRespSDD->GetDeltaVDrift(sID, sddSide);
4720       if (fIniRespSDD->IsVDCorrMult()) vdrift *= (1-corr1);
4721       else vdrift -= corr1*1e-4;
4722     }
4723     tdif    = pnt->GetDriftTime() - t0Upd;
4724     // correct Xlocal
4725     fMeasLoc[0] = fSegmentationSDD->Dx()*1e-4 - vdrift*tdif;
4726     if (sddSide) fMeasLoc[0] = -fMeasLoc[0];
4727     fDriftTime0[pntID] =  t0Upd;
4728   }
4729   //
4730   if (fPreCorrMapSDD) { // apply correction map
4731     fMeasLoc[0] += fPreCorrMapSDD->GetCorrection(sID,fMeasLoc[2],fMeasLoc[0]);
4732   }
4733
4734   // TEMPORARY CORRECTION (if provided) --------------<<<
4735   fDriftSpeed[pntID]  = sddSide ? -vdrift  : vdrift;
4736   fDriftSpeed0[pntID] = sddSide ? -vdrift0 : vdrift0;
4737   //
4738   //  printf("#%d: t:%+e x:%+e v:%+e: side:%d\n",pntID,fDriftTime0[pntID],fMeasLoc[0],fDriftSpeed[pntID],sddSide);
4739 }
4740
4741 //_______________________________________________________________________________________
4742 AliITSAlignMille2Module* AliITSAlignMille2::CreateVertexModule()
4743 {
4744   // creates dummy module for vertex constraint
4745   TGeoHMatrix mt;
4746   AliITSAlignMille2Module* mod = new AliITSAlignMille2Module(kVtxSensID,kVtxSensVID,"VTX",&mt,0,0);
4747   fMilleModule.AddAtAndExpand(mod,fNModules);
4748   mod->SetGeomParamsGlobal(fUseGlobalDelta);
4749   fDiamondModID = fNModules;
4750   mod->SetUniqueID(fNModules++);
4751   mod->SetNotInConf(kTRUE);
4752   return mod;
4753   //
4754 }
4755
4756 //_______________________________________________________________________________________
4757 AliCDBEntry* AliITSAlignMille2::GetCDBEntry(const char* path)
4758 {
4759   // return object from the OCDB
4760   AliCDBEntry *entry = 0;
4761   AliInfo(Form("Loading object %s",path));
4762   AliCDBManager* man = AliCDBManager::Instance();
4763   AliCDBId* cdbId = AliCDBId::MakeFromString(path);
4764   if (!cdbId) {
4765     AliError("Failed to create cdbId");
4766     return 0;
4767   }
4768   //
4769   AliCDBStorage* stor = man->GetDefaultStorage();
4770   if (!stor && !man->GetRaw()) man->SetDefaultStorage("raw://");
4771   if (man->GetRaw()) man->SetRun(fRunID>0 ? fRunID : cdbId->GetFirstRun());
4772   if (stor) {
4773     TString tp = stor->GetType();
4774     if (tp.Contains("alien",TString::kIgnoreCase) && !gGrid) TGrid::Connect("alien:");
4775   }
4776   entry = man->Get(cdbId->GetPath(),cdbId->GetFirstRun(),cdbId->GetVersion(),cdbId->GetSubVersion());
4777   //  entry = man->Get( *cdbId );
4778   man->ClearCache();
4779   //
4780   delete cdbId;
4781   return entry;
4782   //
4783 }
4784
4785 //_______________________________________________________________________________________
4786 void AliITSAlignMille2::SetVertexConstraint(const AliESDVertex* vtx)
4787 {
4788   // set vertex for constraint
4789   if (!vtx) return;
4790   //
4791   double cmat[6];
4792   float cmatF[6];
4793   vtx->GetCovMatrix(cmat);
4794   AliITSAlignMille2Module* diamMod = GetMilleModuleByVID(kVtxSensVID);
4795   if (diamMod) {
4796     cmat[0] *= diamMod->GetSigmaXFactor()*diamMod->GetSigmaXFactor();
4797     cmat[2] *= diamMod->GetSigmaYFactor()*diamMod->GetSigmaYFactor();
4798     cmat[5] *= diamMod->GetSigmaZFactor()*diamMod->GetSigmaZFactor();
4799     cmat[1] *= diamMod->GetSigmaXFactor()*diamMod->GetSigmaYFactor();
4800     cmat[3] *= diamMod->GetSigmaXFactor()*diamMod->GetSigmaZFactor();
4801     cmat[4] *= diamMod->GetSigmaYFactor()*diamMod->GetSigmaZFactor();
4802   }
4803   cmatF[0] = cmat[0]; // xx
4804   cmatF[1] = cmat[1]; // xy
4805   cmatF[2] = cmat[3]; // xz
4806   cmatF[3] = cmat[2]; // yy
4807   cmatF[4] = cmat[4]; // yz
4808   cmatF[5] = cmat[5]; // zz
4809
4810   fDiamond.SetXYZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), cmatF);
4811   //
4812   Double_t t0 = cmat[2]*cmat[5] - cmat[4]*cmat[4];
4813   Double_t t1 = cmat[1]*cmat[5] - cmat[3]*cmat[4];
4814   Double_t t2 = cmat[1]*cmat[4] - cmat[2]*cmat[3];
4815   Double_t det = cmat[0]*t0 - cmat[1]*t1 + cmat[3]*t2;
4816   if (TMath::Abs(det)<1e-36) {
4817     vtx->Print();
4818     AliFatal("Vertex constraint cov.matrix is singular");
4819   }
4820   cmatF[0] =  t0/det;
4821   cmatF[1] = -t1/det;
4822   cmatF[2] =  t2/det;
4823   cmatF[3] =  (cmat[0]*cmat[5] - cmat[3]*cmat[3])/det;
4824   cmatF[4] =  (cmat[1]*cmat[3] - cmat[0]*cmat[4])/det;
4825   cmatF[5] =  (cmat[0]*cmat[2] - cmat[1]*cmat[1])/det;
4826   fDiamondI.SetXYZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), cmatF);
4827   fVertexSet = kTRUE;
4828   //
4829 }
4830
4831 //_______________________________________________________________________________________
4832 void AliITSAlignMille2::ConvertDeltas()
4833 {
4834   // convert prealignment deltas from old geometry to new one
4835   // NOTE: the target geometry must be loaded at time this method is called
4836   //
4837   // NOTE: This method can be ONLY used when as a prealignment deltas those used for the production
4838   // of trackpoints (e.g. extracted from the UserInfo).
4839   // The prealignment deltas provided by user via config file must be already converted to target geometry:
4840   // this can be done externally using the macro ConvertDeltas.C
4841   //
4842   // delta_j_new = delta_j_old * Xj_old * Xj_new^-1
4843   // where X = Prod{delta_i,i=j-1:0} M_j
4844   // with j - the level of the alignable volume in the hierarchy, M - corresponding ideal matrix
4845   // Note that delta_j * Xj is equal to final (misaligned) matrix of corresponding geometry, G_j.
4846   // Since this method is used ONLY in the case where the prealignment deltas are equal to production deltas,
4847   // we have already loaded G_j_old in the fConvAlgMatOld (filled in the CacheMatricesOrig)
4848   // Hence, delta_j_new = G_j_old * Xj_new^-1
4849   //
4850   AliInfo("Converting deltas from initial to target geometry");
4851   int nMatOld = fConvAlgMatOld.GetEntriesFast(); // number of alignable matrices
4852   TClonesArray* deltArrNew = new TClonesArray("AliAlignObjParams",10);
4853   //
4854   TGeoHMatrix dmPar;
4855   int nDelNew = 0;
4856   //
4857   for (int im=0;im<nMatOld;im++) {
4858     TGeoHMatrix* mtGjold = (TGeoHMatrix*)fConvAlgMatOld[im];
4859     TString algname = mtGjold->GetTitle();
4860     UShort_t vID = AliITSAlignMille2Module::GetVolumeIDFromSymname(algname.Data());
4861     //
4862     // build X_new >>>
4863     TGeoHMatrix* parent = mtGjold;
4864     TGeoHMatrix xNew;
4865     int parID;
4866     while ( (parID=parent->GetUniqueID()-1)>=0 ) {
4867       parent = (TGeoHMatrix*)fConvAlgMatOld[parID];
4868       AliAlignObjParams* deltaPar = ConvFindDelta(deltArrNew,parent->GetTitle());
4869       if (deltaPar) deltaPar->GetMatrix(dmPar); xNew *= dmPar;
4870     }
4871     AliGeomManager::GetOrigGlobalMatrix(algname,dmPar); // ideal matrix of new geometry
4872     xNew *= dmPar;
4873     // build X_new <<<
4874     //
4875     dmPar  = *mtGjold;
4876     dmPar *= xNew.Inverse();
4877     new((*deltArrNew)[nDelNew++]) AliAlignObjParams(algname.Data(),vID,dmPar,kTRUE);
4878     //
4879   }
4880   delete fPrealignment;
4881   fPrealignment = deltArrNew;
4882   //
4883   // we don't need anymore old matrices
4884   fConvAlgMatOld.Delete();
4885   //
4886 }
4887
4888 //_______________________________________________________________________________________
4889 void AliITSAlignMille2::ConvSortHierarchically(TObjArray& matArr)
4890 {
4891   // Used only for the deltas conversion from one geometry to another
4892   // Sort the matrices according to hiearachy (coarse -> fine)
4893   //
4894   int nmat = matArr.GetEntriesFast();
4895   //
4896   for (int i=0;i<nmat;i++) {
4897     for (int j=i+1;j<nmat;j++) {
4898       TGeoHMatrix* matI = (TGeoHMatrix*) matArr[i];
4899       TGeoHMatrix* matJ = (TGeoHMatrix*) matArr[j];
4900       if (ConvIsJParentOfI(matI,matJ)) { // swap
4901         matArr[i] = matJ;
4902         matArr[j] = matI;
4903       }
4904     }
4905   }
4906   //
4907   // set direct parent id's in the UniqueID's
4908   for (int i=nmat;i--;) {
4909     TGeoHMatrix* matI = (TGeoHMatrix*) matArr[i];
4910     matI->SetUniqueID(0);
4911     for (int j=i;j--;) {
4912       TGeoHMatrix* matJ = (TGeoHMatrix*) matArr[j];
4913       if (ConvIsJParentOfI(matI,matJ)) { matI->SetUniqueID(j+1); break; }
4914     }
4915   }
4916 }
4917
4918 //_______________________________________________________________________________________
4919 Bool_t AliITSAlignMille2::ConvIsJParentOfI(const TGeoHMatrix* matI,const TGeoHMatrix* matJ) const
4920 {
4921   // Used only for the deltas conversion from one geometry to another
4922   // True if matJ is higher in hierarchy than
4923   //
4924   TString nmI = matI->GetTitle();
4925   TString nmJ = matJ->GetTitle();
4926   //
4927   int nlrI = nmI.CountChar('/');
4928   int nlrJ = nmJ.CountChar('/');
4929   if (nlrJ>=nlrI) return kFALSE;
4930   //
4931   // special case of SPD sectors
4932   if (nmI.BeginsWith("ITS/SPD1") && nmJ.BeginsWith("ITS/SPD0") && nlrJ==2) nmJ.ReplaceAll("SPD0","SPD1");
4933   return (nmI.BeginsWith(nmJ)) ? kTRUE:kFALSE;
4934   //
4935 }
4936
4937 //_______________________________________________________________________________________
4938 AliAlignObjParams* AliITSAlignMille2::ConvFindDelta(const TClonesArray* arrDelta,const TString& algname) const
4939 {
4940   // find the delta for given module
4941   if (!arrDelta) return 0;
4942   AliAlignObjParams* delta = 0;
4943   int nDeltas = arrDelta->GetEntries();
4944   for (int id=0;id<nDeltas;id++) {
4945     delta = (AliAlignObjParams*)arrDelta->At(id);
4946     if (algname==delta->GetSymName()) break;
4947     delta = 0;
4948   }
4949   return delta;
4950 }