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