ITS/AliITSAlignMille2.cxx

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