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