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