ITS/AliITSAlignMille2.cxx

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