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