ITS/AliITSAlignMille.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 /// \class AliITSAlignMille
  20 /// Alignment class fro 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)
  28 //-----------------------------------------------------------------------------
  29
  30 #include <TF1.h>
  31 #include <TGraph.h>
  32 #include <TGeoMatrix.h>
  33 #include <TMath.h>
  34
  35 #include "AliITSAlignMille.h"
  36 #include "AliITSgeomTGeo.h"
  37 #include "AliGeomManager.h"
  38 #include "AliMillepede.h"
  39 #include "AliTrackPointArray.h"
  40 #include "AliAlignObjParams.h"
  41 #include "AliLog.h"
  42 #include "TSystem.h"  // come si fa?
  43
  44 /// \cond CLASSIMP
  45 ClassImp(AliITSAlignMille)
  46 /// \endcond
  47
  48 Int_t AliITSAlignMille::fgNDetElem = ITSMILLE_NDETELEM;
  49 Int_t AliITSAlignMille::fgNParCh = ITSMILLE_NPARCH;
  50
  51 AliITSAlignMille::AliITSAlignMille(const Char_t *configFilename, Bool_t initmille)
  52   : TObject(),
  53     fMillepede(0),
  54     fStartFac(16.),
  55     fResCutInitial(100.),
  56     fResCut(100.),
  57     fNGlobal(ITSMILLE_NDETELEM*ITSMILLE_NPARCH),
  58     fNLocal(ITSMILLE_NLOCAL),
  59     fNStdDev(ITSMILLE_NSTDEV),
  60     fIsMilleInit(kFALSE),
  61     fParSigTranslations(0.0100),
  62     fParSigRotations(0.1),
  63     fTrack(NULL),
  64     fCluster(),
  65     fGlobalDerivatives(NULL),
  66     fTempHMat(NULL),
  67     fTempAlignObj(NULL),
  68     fDerivativeXLoc(0),
  69     fDerivativeZLoc(0),
  70     fDeltaPar(0),
  71     fMinNPtsPerTrack(3),
  72     fGeometryFileName("geometry.root"),
  73     fGeoManager(0),
  74     fCurrentModuleIndex(0),
  75     fCurrentModuleInternalIndex(0),
  76     fNModules(0),
  77     fUseLocalShifts(kTRUE),
  78     fCurrentModuleHMatrix(NULL)
  79 {
  80   /// main constructor that takes input from configuration file
  81
  82   fMillepede = new AliMillepede();
  83   fGlobalDerivatives = new Double_t[fNGlobal];
  84   fTempHMat = new TGeoHMatrix;
  85   fCurrentModuleHMatrix = new TGeoHMatrix;
  86
  87   Int_t lc=LoadConfig(configFilename);
  88   if (lc) {
  89     AliInfo(Form("Error %d loading configuration from %s",lc,configFilename));
  90   }
  91   else {
  92     if (initmille && fNGlobal<20000) {
  93       AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, fNLocal, fNStdDev));
  94       Init(fNGlobal, fNLocal, fNStdDev);
  95       ResetLocalEquation();
  96       AliInfo("Parameters initialized to zero");
  97     }
  98     else {
  99       AliInfo("Millepede has not been initialized ... ");
 100     }
 101   }
 102
 103   fDeltaPar=0.0; // not used at the moment - to be checked later...
 104
 105 }
 106
 107 AliITSAlignMille::~AliITSAlignMille() {
 108   /// Destructor
 109   if (fMillepede) delete fMillepede;
 110   delete [] fGlobalDerivatives;
 111   delete fCurrentModuleHMatrix;
 112   delete fTempHMat;
 113 }
 114
 115 ///////////////////////////////////////////////////////////////////////
 116
 117 Int_t AliITSAlignMille::LoadConfig(const Char_t *cfile) {
 118   /// return 0 if success
 119   ///        1 if error in module index or voluid
 120
 121   FILE *pfc=fopen(cfile,"r");
 122   if (!pfc) return -1;
 123
 124   Char_t st[200],st2[200];
 125   Char_t tmp[100];
 126   Int_t idx,itx,ity,itz,ith,ips,iph;
 127   Float_t f1;
 128   UShort_t voluid;
 129   Int_t nmod=0;
 130
 131   while (fgets(st,200,pfc)) {
 132
 133     // skip comments
 134     for (int i=0; i<int(strlen(st)); i++) {
 135       if (st[i]=='#') st[i]=0;
 136     }
 137
 138     if (strstr(st,"GEOMETRY_FILE")) {
 139       sscanf(st,"%s %s",tmp,st2);
 140       if (gSystem->AccessPathName(st2)) {
 141         AliInfo("*** WARNING! *** geometry file not found! ");
 142         return -1;
 143       }
 144       fGeometryFileName=st2;
 145       InitGeometry();
 146     }
 147
 148     if (strstr(st,"SET_PARSIG_TRA")) {
 149       sscanf(st,"%s %f",tmp,&f1);
 150       fParSigTranslations=f1;
 151     }
 152
 153     if (strstr(st,"SET_PARSIG_ROT")) {
 154       sscanf(st,"%s %f",tmp,&f1);
 155       fParSigRotations=f1;
 156     }
 157
 158     if (strstr(st,"SET_NSTDDEV")) {
 159       sscanf(st,"%s %d",tmp,&idx);
 160       fNStdDev=idx;
 161     }
 162
 163     if (strstr(st,"SET_RESCUT_INIT")) {
 164       sscanf(st,"%s %f",tmp,&f1);
 165       fResCutInitial=f1;
 166     }
 167
 168     if (strstr(st,"SET_RESCUT_OTHER")) {
 169       sscanf(st,"%s %f",tmp,&f1);
 170       fResCut=f1;
 171     }
 172
 173     if (strstr(st,"SET_STARTFAC")) {
 174       sscanf(st,"%s %f",tmp,&f1);
 175       fStartFac=f1;
 176     }
 177
 178     if (strstr(st,"SET_LOCAL_SHIFTS")) { // only enabled mode...
 179       fUseLocalShifts = kTRUE;
 180     }
 181
 182     if (strstr(st,"MODULE_INDEX")) {
 183       sscanf(st,"%s %d %d %d %d %d %d %d",tmp,&idx,&itx,&ity,&itz,&iph,&ith,&ips);
 184       voluid=GetModuleVolumeID(idx);
 185       if (!voluid) return 1; // bad index
 186       fModuleIndex[nmod]=idx;
 187       fModuleVolumeID[nmod]=voluid;
 188       fFreeParam[nmod][0]=itx;
 189       fFreeParam[nmod][1]=ity;
 190       fFreeParam[nmod][2]=itz;
 191       fFreeParam[nmod][3]=iph;
 192       fFreeParam[nmod][4]=ith;
 193       fFreeParam[nmod][5]=ips;
 194       nmod++;
 195     }
 196
 197     if (strstr(st,"MODULE_VOLUID")) {
 198       // to be implemented
 199     }
 200
 201   } // end while
 202
 203   fNModules = nmod;
 204   fNGlobal = fNModules*fgNParCh;
 205
 206   fclose(pfc);
 207   return 0;
 208 }
 209
 210 Int_t AliITSAlignMille::GetModuleIndex(const Char_t *symname) {
 211   /// index from symname
 212   if (!symname) return -1;
 213   for (Int_t i=0; i<2198; i++) {
 214     if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i;
 215   }
 216   return -1;
 217 }
 218
 219 Int_t AliITSAlignMille::GetModuleIndex(UShort_t voluid) {
 220   /// index from volume ID
 221   AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
 222   if (lay<1|| lay>6) return -1;
 223   Int_t idx=Int_t(voluid)-2048*lay;
 224   if (idx>=AliGeomManager::LayerSize(lay)) return -1;
 225   for (Int_t ilay=1; ilay<lay; ilay++)
 226     idx += AliGeomManager::LayerSize(ilay);
 227   return idx;
 228 }
 229
 230 UShort_t AliITSAlignMille::GetModuleVolumeID(const Char_t *symname) {
 231   /// volume ID from symname
 232   if (!symname) return 0;
 233
 234   for (UShort_t voluid=2000; voluid<13300; voluid++) {
 235     Int_t modId;
 236     AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
 237     if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
 238       if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
 239     }
 240   }
 241
 242   return 0;
 243 }
 244
 245 UShort_t AliITSAlignMille::GetModuleVolumeID(Int_t index) {
 246   /// volume ID from index
 247   if (index<0 || index>2197) return 0;
 248   return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index));
 249 }
 250
 251 void AliITSAlignMille::InitGeometry() {
 252   /// initialize geometry
 253   AliGeomManager::LoadGeometry(fGeometryFileName.Data());
 254   fGeoManager = AliGeomManager::GetGeometry();
 255   if (!fGeoManager) {
 256     AliInfo("Couldn't initialize geometry");
 257     return;
 258   }
 259   // temporary align object, just use the rotation...
 260   fTempAlignObj=new AliAlignObjParams(AliITSgeomTGeo::GetSymName(7),2055,0,0,0,0,0,0,kFALSE);
 261 }
 262
 263 void AliITSAlignMille::Init(Int_t nGlobal,  /* number of global paramers */
 264                            Int_t nLocal,   /* number of local parameters */
 265                            Int_t nStdDev   /* std dev cut */ )
 266 {
 267   /// Initialization of AliMillepede. Fix parameters, define constraints ...
 268   fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
 269   fIsMilleInit = kTRUE;
 270
 271   /// Fix non free parameters
 272   for (Int_t i=0; i<fNModules; i++) {
 273     for (Int_t j=0; j<ITSMILLE_NPARCH; j++) {
 274       if (!fFreeParam[i][j]) FixParameter(i*ITSMILLE_NPARCH+j,0.0);
 275       else {
 276         // pepopepo: da sistemare il settaggio delle sigma...
 277         Double_t parsig=0;
 278         if (j<3) parsig=fParSigTranslations; // translations (0.0100 cm)
 279         else parsig=fParSigRotations; // rotations (1/10 deg)
 280         FixParameter(i*ITSMILLE_NPARCH+j,parsig);
 281       }
 282     }
 283   }
 284
 285
 286 //   // Set iterations
 287   if (fStartFac>1) fMillepede->SetIterations(fStartFac);
 288 }
 289
 290
 291 void AliITSAlignMille::AddConstraint(Double_t *par, Double_t value) {
 292   /// Constrain equation defined by par to value
 293   if (!fIsMilleInit) {
 294     AliInfo("Millepede has not been initialized!");
 295     return;
 296   }
 297   fMillepede->SetGlobalConstraint(par, value);
 298   AliInfo("Adding constraint");
 299 }
 300
 301 void AliITSAlignMille::InitGlobalParameters(Double_t *par) {
 302   /// Initialize global parameters with par array
 303   if (!fIsMilleInit) {
 304     AliInfo("Millepede has not been initialized!");
 305     return;
 306   }
 307   fMillepede->SetGlobalParameters(par);
 308   AliInfo("Init Global Parameters");
 309 }
 310
 311 void AliITSAlignMille::FixParameter(Int_t iPar, Double_t value) {
 312   /// Parameter iPar is encourage to vary in [-value;value].
 313   /// If value == 0, parameter is fixed
 314   if (!fIsMilleInit) {
 315     AliInfo("Millepede has not been initialized!");
 316     return;
 317   }
 318   fMillepede->SetParSigma(iPar, value);
 319   if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
 320 }
 321
 322 void AliITSAlignMille::ResetLocalEquation()
 323 {
 324   /// Reset the derivative vectors
 325   for(int i=0; i<fNLocal; i++) {
 326     fLocalDerivatives[i] = 0.0;
 327   }
 328   for(int i=0; i<fNGlobal; i++) {
 329     fGlobalDerivatives[i] = 0.0;
 330   }
 331 }
 332
 333 Int_t AliITSAlignMille::InitModuleParams() {
 334   /// initialize geometry parameters for a given detector
 335   /// for current cluster (fCluster)
 336   /// fGlobalInitParam[] is set as:
 337   ///    [tx,ty,tz,psi,theta,phi]
 338   ///    (old was [tx,ty,tz,theta,psi,phi] ROOT's angles...)
 339   /// *** At the moment: using Raffalele's angles definition ***
 340   ///
 341   /// Num of Dets: ITSMILLE_NDETELEM = fgNDetElem
 342   /// Num of Par : ITSMILLE_NPARCH = fgNParCh
 343   /// return 0 if success
 344
 345   if (!fGeoManager) {
 346     AliInfo("ITS geometry not initialized!");
 347     return -1;
 348   }
 349
 350   // set the internal index (index in module list)
 351   UShort_t voluid=fCluster.GetVolumeID();
 352   Int_t k=fNModules-1;
 353   while (k>=0 && !(voluid==fModuleVolumeID[k]) ) k--;
 354   if (k<0) return -3;
 355   fCurrentModuleInternalIndex=k;
 356
 357   // set the index
 358   Int_t index = GetModuleIndex(AliGeomManager::SymName(voluid));
 359   if (index<0) return -2;
 360   fCurrentModuleIndex = index;
 361
 362   fModuleInitParam[0] = 0.0;
 363   fModuleInitParam[1] = 0.0;
 364   fModuleInitParam[2] = 0.0;
 365   fModuleInitParam[3] = 0.0; // psi   (X)
 366   fModuleInitParam[4] = 0.0; // theta (Y)
 367   fModuleInitParam[5] = 0.0; // phi   (Z)
 368
 369   /// get global (corrected) matrix
 370   //  if (!AliITSgeomTGeo::GetOrigMatrix(index,*fCurrentModuleHMatrix)) return -3;
 371   Double_t rott[9];
 372   if (!AliITSgeomTGeo::GetRotation(index,rott)) return -3;
 373   fCurrentModuleHMatrix->SetRotation(rott);
 374   Double_t oLoc[3]={0,0,0};
 375   if (!AliITSgeomTGeo::LocalToGlobal(index,oLoc,fCurrentModuleTranslation)) return -4;
 376   fCurrentModuleHMatrix->SetTranslation(fCurrentModuleTranslation);
 377
 378   /// get back local coordinates
 379   fMeasGlo[0] = fCluster.GetX();
 380   fMeasGlo[1] = fCluster.GetY();
 381   fMeasGlo[2] = fCluster.GetZ();
 382   fCurrentModuleHMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
 383
 384   // set stdev from cluster
 385   TGeoHMatrix hcov;
 386   Double_t hcovel[9];
 387   hcovel[0]=double(fCluster.GetCov()[0]);
 388   hcovel[1]=double(fCluster.GetCov()[1]);
 389   hcovel[2]=double(fCluster.GetCov()[3]);
 390   hcovel[3]=double(fCluster.GetCov()[1]);
 391   hcovel[4]=double(fCluster.GetCov()[2]);
 392   hcovel[5]=double(fCluster.GetCov()[4]);
 393   hcovel[6]=double(fCluster.GetCov()[3]);
 394   hcovel[7]=double(fCluster.GetCov()[4]);
 395   hcovel[8]=double(fCluster.GetCov()[5]);
 396   hcov.SetRotation(hcovel);
 397   // now rotate in local system
 398   hcov.MultiplyLeft(&fCurrentModuleHMatrix->Inverse());
 399   hcov.Multiply(fCurrentModuleHMatrix);
 400
 401   // per i ruotati c'e' delle sigmaY che compaiono... prob
 402   // e' un problema di troncamento
 403   fSigmaLoc[0] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[0]));
 404   fSigmaLoc[1] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[4]));
 405   fSigmaLoc[2] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[8]));
 406
 407     AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%f  fSigmaLocY=%f fSigmaLocZ=%f \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
 408
 409   return 0;
 410 }
 411
 412 void AliITSAlignMille::SetCurrentModule(Int_t index) {
 413   ///
 414   UShort_t voluid=GetModuleVolumeID(index);
 415   if (voluid) {
 416     fCluster.SetVolumeID(voluid);
 417     fCluster.SetXYZ(0,0,0);
 418     InitModuleParams();
 419   }
 420 }
 421
 422 void AliITSAlignMille::Print() {
 423   ///
 424   printf("*** AliMillepede for ITS ***\n");
 425   printf("    number of defined modules: %d\n",fNModules);
 426   if (fGeoManager)
 427     printf("    geometry loaded from %s\n",fGeometryFileName.Data());
 428   else
 429     printf("    geometry not loaded\n");
 430   if (fUseLocalShifts)
 431     printf("    Alignment shifts will be computed in LOCAL RS\n");
 432   else
 433     printf("    Alignment shifts will be computed in GLOBAL RS\n");
 434   printf("    Millepede configuration parameters:\n");
 435   printf("       init parsig for translations  : %.4f\n",fParSigTranslations);
 436   printf("       init parsig for rotations     : %.4f\n",fParSigRotations);
 437   printf("       init value for chi2 cut       : %.4f\n",fStartFac);
 438   printf("       first iteration cut value     : %.4f\n",fResCutInitial);
 439   printf("       other iterations cut value    : %.4f\n",fResCut);
 440   printf("       number of stddev for chi2 cut : %d\n",fNStdDev);
 441
 442   printf("List of defined modules:\n");
 443   printf("  intidx\tindex\tvoluid\tsymname\n");
 444   for (int i=0; i<fNModules; i++)
 445     printf("  %d\t%d\t%d\t%s\n",i,fModuleIndex[i],fModuleVolumeID[i],AliITSgeomTGeo::GetSymName(fModuleIndex[i]));
 446 }
 447
 448 void AliITSAlignMille::PrintCurrentModuleInfo() {
 449   ///
 450   if (fCurrentModuleIndex<0 || fCurrentModuleIndex>2197) return;
 451   UShort_t voluid=fModuleVolumeID[fCurrentModuleInternalIndex];
 452   printf("Current module: index=%d   voluid=%d\n",fCurrentModuleIndex,voluid);
 453   printf("                symname:%s\n",AliGeomManager::SymName(voluid));
 454   printf("  TGeoHMatrix: \n");
 455   fCurrentModuleHMatrix->Print();
 456 }
 457
 458
 459 void AliITSAlignMille::InitTrackParams(int meth) {
 460   /// initialize local parameters with different methods
 461   /// for current track (fTrack)
 462
 463   switch (meth) {
 464   case 1:   // simple linear interpolation
 465     // get local starting parameters (to be substituted by ESD track parms)
 466     // local parms (fLocalInitParam[]) are:
 467     //      [0] = global x coord. of straight line intersection at y=0 plane
 468     //      [1] = global z coord. of straight line intersection at y=0 plane
 469     //      [2] = px/py
 470     //      [3] = pz/py
 471
 472     // test #1: linear fit in x(y) and z(y)
 473     Int_t npts = fTrack->GetNPoints();
 474
 475     TF1 *f1=new TF1("f1","[0]+x*[1]",-50,50);
 476
 477     TGraph *g=new TGraph(npts,fTrack->GetY(),fTrack->GetX());
 478     g->Fit(f1,"RNQ");
 479     fLocalInitParam[0] = f1->GetParameter(0);
 480     fLocalInitParam[2] = f1->GetParameter(1);
 481     AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f +- %f    ugx = %f +- %f\n",fLocalInitParam[0],f1->GetParError(0),fLocalInitParam[2],f1->GetParError(1)));
 482     delete g; g=NULL;
 483
 484     g=new TGraph(npts,fTrack->GetY(),fTrack->GetZ());
 485     g->Fit(f1,"RNQ");
 486     fLocalInitParam[1] = f1->GetParameter(0);
 487     fLocalInitParam[3] = f1->GetParameter(1);
 488     AliDebug(2,Form("Z = p0gz + ugz*Y : p0gz=%f  ugz=%f\n",fLocalInitParam[1],fLocalInitParam[3]));
 489     delete g;
 490     delete f1;
 491
 492     break;
 493   }
 494
 495 }
 496 Bool_t AliITSAlignMille::CheckVolumeID(UShort_t voluid) const
 497 {
 498   ///
 499   Int_t k=fNModules-1;
 500   while (k>=0 && !(voluid==fModuleVolumeID[k]) ) k--;
 501   //printf("selected element with voluid=%d : %d\n",voluid,k);
 502   if (k>=0) return kTRUE;
 503   return kFALSE;
 504 }
 505
 506 Int_t AliITSAlignMille::CheckCurrentTrack() {
 507   /// checks if AliTrackPoints belongs to defined modules
 508   /// return number of good poins
 509   /// return 0 if not enough points
 510
 511   Int_t npts = fTrack->GetNPoints();
 512   Int_t ngoodpts=0;
 513   // debug points
 514   for (int j=0; j<npts; j++) {
 515     UShort_t voluid = fTrack->GetVolumeID()[j];
 516     if (CheckVolumeID(voluid)) {
 517       ngoodpts++;
 518     }
 519   }
 520   // pepo da controllare...
 521   if (ngoodpts<fMinNPtsPerTrack) return 0;
 522
 523   return ngoodpts;
 524 }
 525
 526 Int_t AliITSAlignMille::ProcessTrack(AliTrackPointArray *track) {
 527   /// Process track; Loop over hits and set local equations
 528   /// here 'track' is a AliTrackPointArray
 529   /// return 0 if success;
 530
 531   if (!fIsMilleInit) {
 532     AliInfo("Millepede has not been initialized!");
 533     return -1;
 534   }
 535
 536   Int_t npts = track->GetNPoints();
 537   AliDebug(2,Form("\n*** Processing track with %d points ***\n",npts));
 538   fTrack = track;
 539
 540   // checks if there are enough good points
 541   if (!CheckCurrentTrack()) {
 542     AliInfo("Track with not enough good points - will not be used...");
 543     return -1;
 544   }
 545
 546   // set local starting parameters (to be substituted by ESD track parms)
 547   // local parms (fLocalInitParam[]) are:
 548   //      [0] = global x coord. of straight line intersection at y=0 plane
 549   //      [1] = global z coord. of straight line intersection at y=0 plane
 550   //      [2] = px/py
 551   //      [3] = pz/py
 552   InitTrackParams(1);
 553
 554   for (Int_t ipt=0; ipt<npts; ipt++) {
 555     fTrack->GetPoint(fCluster,ipt);
 556     if (!CheckVolumeID(fCluster.GetVolumeID())) continue;
 557
 558     // set geometry parameters for the the current module
 559     AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
 560     if (InitModuleParams()) continue;
 561
 562     AliDebug(2,Form("    VolID=%d  Index=%d  InternalIdx=%d  symname=%s\n", track->GetVolumeID()[ipt], fCurrentModuleIndex ,fCurrentModuleInternalIndex, AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
 563     AliDebug(2,Form("    Point = ( %f , %f , %f ) \n",track->GetX()[ipt],track->GetY()[ipt],track->GetZ()[ipt]));
 564
 565     if (SetLocalEquations()) return -1;
 566
 567   } // end loop over points
 568
 569   return 0;
 570 }
 571
 572 Int_t AliITSAlignMille::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar) {
 573   /// calculate track intersection point in local coordinates
 574   /// according with a given set of parameters (local(4) and global(6))
 575   /// and fill fPintLoc/Glo
 576   ///    local are:   pgx0, pgz0, ugx0, ugz0
 577   ///    global are:  tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
 578   /// return 0 if success
 579
 580   AliDebug(3,Form("lpar = %g %g %g %g \ngpar= %g %g %g %g %g %g\n",lpar[0],lpar[1],lpar[2],lpar[3],gpar[0],gpar[1],gpar[2],gpar[3],gpar[4],gpar[5]));
 581
 582   // vector of initial straight line direction in glob. coord
 583   // ATTENZIONE: forse va rinormalizzato tutto...
 584   Double_t v0g[3];
 585   //Double_t
 586   v0g[0]=lpar[2];
 587   v0g[1]=1.0;
 588   v0g[2]=lpar[3];
 589
 590   // intercept in yg=0 plane in glob coord
 591   Double_t p0g[3];
 592   p0g[0]=lpar[0];
 593   p0g[1]=0.0;
 594   p0g[2]=lpar[1];
 595
 596
 597   // prepare the TGeoHMatrix
 598   Double_t tr[3],ang[3];
 599   //Double_t rad2deg=180./TMath::Pi();
 600   if (fUseLocalShifts) { // just Delta matrix
 601     tr[0]=gpar[0];
 602     tr[1]=gpar[1];
 603     tr[2]=gpar[2];
 604     ang[0]=gpar[3]; // psi   (X)
 605     ang[1]=gpar[4]; // theta (Y)
 606     ang[2]=gpar[5]; // phi   (Z)
 607   }
 608   else { // total matrix with shifted parameter
 609     AliInfo("global shifts not implemented yet!");
 610     return -1;
 611   }
 612
 613   //printf("fTempRot = 0x%x  - ang = %g %g %g \n",fTempRot,gpar[5]*rad2deg,gpar[3]*rad2deg,gpar[4]*rad2deg);
 614
 615   fTempAlignObj->SetRotation(ang[0],ang[1],ang[2]);
 616   AliDebug(3,Form("Delta angles: psi=%f  theta=%f   phi=%f",ang[0],ang[1],ang[2]));
 617   TGeoHMatrix hm;
 618   fTempAlignObj->GetMatrix(hm);
 619   fTempHMat->SetRotation(hm.GetRotationMatrix());
 620   fTempHMat->SetTranslation(tr);
 621
 622   // in this case the gpar[] array contains only shifts
 623   // and fInitModuleParam[] are set to 0
 624   // fCurrentModuleHMatrix is then modified as fCurrentHM*fTempHM
 625   if (fUseLocalShifts)
 626     fTempHMat->MultiplyLeft(fCurrentModuleHMatrix);
 627
 628   // same in local coord.
 629   Double_t p0l[3],v0l[3];
 630   fTempHMat->MasterToLocalVect(v0g,v0l);
 631   fTempHMat->MasterToLocal(p0g,p0l);
 632
 633   if (TMath::Abs(v0l[1])<1e-15) {
 634     AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
 635     return -1;
 636   }
 637
 638   // local intersection point
 639   fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
 640   fPintLoc[1] = 0;
 641   fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
 642
 643   // global intersection point
 644   fTempHMat->LocalToMaster(fPintLoc,fPintGlo);
 645   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]));
 646
 647   return 0;
 648 }
 649
 650 Int_t AliITSAlignMille::CalcDerivatives(Int_t paridx, Bool_t islpar) {
 651   /// calculate numerically (ROOT's style) the derivatives for
 652   /// local X intersection  and local Z intersection
 653   /// parlist: local  (islpar=kTRUE)  pgx0, pgz0, ugx0, ugz0
 654   ///          global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
 655   /// return 0 if success
 656
 657   // copy initial parameters
 658   Double_t lpar[ITSMILLE_NLOCAL];
 659   Double_t gpar[ITSMILLE_NPARCH];
 660   for (Int_t i=0; i<ITSMILLE_NLOCAL; i++) lpar[i]=fLocalInitParam[i];
 661   for (Int_t i=0; i<ITSMILLE_NPARCH; i++) gpar[i]=fModuleInitParam[i];
 662
 663   // pepopepo
 664   // trial with fixed dpar...
 665   Double_t dpar=0.0;
 666   if (islpar) {
 667     //dpar=fLocalInitParam[paridx]*0.001;
 668     // set minimum dpar
 669     if (paridx<2) dpar=1.0e-4; // translations
 670     else dpar=1.0e-6; // direction
 671   }
 672   else {
 673     //dpar=fModuleInitParam[paridx]*0.001;
 674     if (paridx<3) dpar=1.0e-4; // translations
 675     else dpar=1.0e-2; // angles
 676   }
 677   AliDebug(3,Form("\n+++ automatic dpar=%g\n",dpar));
 678   if (fDeltaPar) dpar=fDeltaPar;
 679   AliDebug(3,Form("+++ using dpar=%g\n\n",dpar));
 680
 681   // calculate derivative ROOT's like:
 682   //  using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
 683   Double_t pintl1[3]; // f(x-h)
 684   Double_t pintl2[3]; // f(x-h/2)
 685   Double_t pintl3[3]; // f(x+h/2)
 686   Double_t pintl4[3]; // f(x+h)
 687
 688   // first values
 689   if (islpar) lpar[paridx] -= dpar;
 690   else gpar[paridx] -= dpar;
 691   if (CalcIntersectionPoint(lpar, gpar)) return -2;
 692   for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
 693
 694   // second values
 695   if (islpar) lpar[paridx] += dpar/2;
 696   else gpar[paridx] += dpar/2;
 697   if (CalcIntersectionPoint(lpar, gpar)) return -2;
 698   for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
 699
 700   // third values
 701   if (islpar) lpar[paridx] += dpar;
 702   else gpar[paridx] += dpar;
 703   if (CalcIntersectionPoint(lpar, gpar)) return -2;
 704   for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
 705
 706   // fourth values
 707   if (islpar) lpar[paridx] += dpar/2;
 708   else gpar[paridx] += dpar/2;
 709   if (CalcIntersectionPoint(lpar, gpar)) return -2;
 710   for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
 711
 712   Double_t h2 = 1./(2.*dpar);
 713   Double_t d0 = pintl4[0]-pintl1[0];
 714   Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
 715   fDerivativeXLoc = h2*(4*d2 - d0)/3.;
 716   if (TMath::Abs(fDerivativeXLoc) < 1.0e-9) fDerivativeXLoc=0.0;
 717
 718   d0 = pintl4[2]-pintl1[2];
 719   d2 = 2.*(pintl3[2]-pintl2[2]);
 720   fDerivativeZLoc = h2*(4*d2 - d0)/3.;
 721   if (TMath::Abs(fDerivativeZLoc) < 1.0e-9) fDerivativeZLoc=0.0;
 722
 723   AliDebug(3,Form("\n+++ derivatives +++ \n"));
 724   AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",fDerivativeXLoc));
 725   AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",fDerivativeZLoc));
 726
 727   return 0;
 728 }
 729
 730 Int_t AliITSAlignMille::SetLocalEquations() {
 731   /// Define local equation for current track and cluster in x coor.
 732   /// return 0 if success
 733
 734   // store first interaction point
 735   CalcIntersectionPoint(fLocalInitParam, fModuleInitParam);
 736   for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
 737
 738   // calculate local derivatives numerically
 739   Double_t dXdL[ITSMILLE_NLOCAL],dZdL[ITSMILLE_NLOCAL];
 740   for (Int_t i=0; i<ITSMILLE_NLOCAL; i++) {
 741     if (CalcDerivatives(i,kTRUE)) return -1;
 742     dXdL[i]=fDerivativeXLoc;
 743     dZdL[i]=fDerivativeZLoc;
 744   }
 745
 746   Double_t dXdG[ITSMILLE_NPARCH],dZdG[ITSMILLE_NPARCH];
 747   for (Int_t i=0; i<ITSMILLE_NPARCH; i++) {
 748     if (CalcDerivatives(i,kFALSE)) return -1;
 749     dXdG[i]=fDerivativeXLoc;
 750     dZdG[i]=fDerivativeZLoc;
 751   }
 752
 753   AliDebug(2,Form("\n***************\n"));
 754   for (Int_t i=0; i<ITSMILLE_NLOCAL; i++)
 755     AliDebug(2,Form("Local parameter %d - dXdpar = %g  - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
 756   for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
 757     AliDebug(2,Form("Global parameter %d - dXdpar = %g  - dZdpar = %g\n",i,dXdG[i],dZdG[i]));
 758   AliDebug(2,Form("\n***************\n"));
 759
 760
 761   AliDebug(2,Form("setting local equation X with fMeas=%.6f  and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
 762   // set equation for Xloc coordinate
 763   for (Int_t i=0; i<ITSMILLE_NLOCAL; i++)
 764     SetLocalDerivative(i,dXdL[i]);
 765   for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
 766     SetGlobalDerivative(fCurrentModuleInternalIndex*ITSMILLE_NPARCH+i,dXdG[i]);
 767   fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, (fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]);
 768
 769
 770   AliDebug(2,Form("setting local equation Z with fMeas=%.6f  and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
 771   // set equation for Zloc coordinate
 772   for (Int_t i=0; i<ITSMILLE_NLOCAL; i++)
 773     SetLocalDerivative(i,dZdL[i]);
 774   for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
 775     SetGlobalDerivative(fCurrentModuleInternalIndex*ITSMILLE_NPARCH+i,dZdG[i]);
 776   fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, (fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]);
 777
 778   return 0;
 779 }
 780
 781
 782 void AliITSAlignMille::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) {
 783   /// Call local fit for this track
 784   if (!fIsMilleInit) {
 785     AliInfo("Millepede has not been initialized!");
 786     return;
 787   }
 788   Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
 789   AliDebug(2,Form("iRes = %d",iRes));
 790   if (iRes && !lSingleFit) {
 791     fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1);
 792   }
 793 }
 794
 795 void AliITSAlignMille::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) {
 796   /// Call global fit; Global parameters are stored in parameters
 797   if (!fIsMilleInit) {
 798     AliInfo("Millepede has not been initialized!");
 799     return;
 800   }
 801   fMillepede->GlobalFit(parameters,errors,pulls);
 802   AliInfo("Done fitting global parameters!");
 803 }
 804
 805 Double_t AliITSAlignMille::GetParError(Int_t iPar) {
 806   /// Get error of parameter iPar
 807   if (!fIsMilleInit) {
 808     AliInfo("Millepede has not been initialized!");
 809     return 0;
 810   }
 811   Double_t lErr = fMillepede->GetParError(iPar);
 812   return lErr;
 813 }
 814
 815 void AliITSAlignMille::PrintGlobalParameters() {
 816   /// Print global parameters
 817   if (!fIsMilleInit) {
 818     AliInfo("Millepede has not been initialized!");
 819     return;
 820   }
 821   fMillepede->PrintGlobalParameters();
 822 }
 823
 824 // //_________________________________________________________________________
 825