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