STEER/AliAlignmentTracks.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //-----------------------------------------------------------------
  17 //   Implementation of the alignment steering class
  18 //   It provides an access to the track space points
  19 //   written along the esd tracks. The class enables
  20 //   the user to plug any track fitter (deriving from
  21 //   AliTrackFitter class) and minimization fo the
  22 //   track residual sums (deriving from the AliTrackResiduals).
  23 //-----------------------------------------------------------------
  24
  25 #include <TChain.h>
  26 #include <TFile.h>
  27
  28 #include "AliAlignmentTracks.h"
  29 #include "AliTrackPointArray.h"
  30 #include "AliAlignObjParams.h"
  31 #include "AliTrackFitterRieman.h"
  32 #include "AliTrackResidualsChi2.h"
  33 #include "AliESDEvent.h"
  34 #include "AliLog.h"
  35
  36 ClassImp(AliAlignmentTracks)
  37
  38 //______________________________________________________________________________
  39 AliAlignmentTracks::AliAlignmentTracks():
  40   fESDChain(0),
  41   fPointsFilename("AliTrackPoints.root"),
  42   fPointsFile(0),
  43   fPointsTree(0),
  44   fLastIndex(0),
  45   fArrayIndex(0),
  46   fIsIndexBuilt(kFALSE),
  47   fAlignObjs(0),
  48   fMisalignObjs(0),
  49   fTrackFitter(0),
  50   fMinimizer(0),
  51   fDoUpdate(kTRUE)
  52 {
  53   // Default constructor
  54   InitIndex();
  55   InitAlignObjs();
  56 }
  57
  58 //______________________________________________________________________________
  59 AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
  60   fESDChain(esdchain),
  61   fPointsFilename("AliTrackPoints.root"),
  62   fPointsFile(0),
  63   fPointsTree(0),
  64   fLastIndex(0),
  65   fArrayIndex(0),
  66   fIsIndexBuilt(kFALSE),
  67   fAlignObjs(0),
  68   fMisalignObjs(0),
  69   fTrackFitter(0),
  70   fMinimizer(0),
  71   fDoUpdate(kTRUE)
  72 {
  73   // Constructor in the case
  74   // the user provides an already
  75   // built TChain with ESD trees
  76   InitIndex();
  77   InitAlignObjs();
  78 }
  79
  80
  81 //______________________________________________________________________________
  82 AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
  83   fESDChain(new TChain(esdtreename)),
  84   fPointsFilename("AliTrackPoints.root"),
  85   fPointsFile(0),
  86   fPointsTree(0),
  87   fLastIndex(0),
  88   fArrayIndex(0),
  89   fIsIndexBuilt(kFALSE),
  90   fAlignObjs(0),
  91   fMisalignObjs(0),
  92   fTrackFitter(0),
  93   fMinimizer(0),
  94   fDoUpdate(kTRUE)
  95 {
  96   // Constructor in the case
  97   // the user provides a single ESD file
  98   // or a directory containing ESD files
  99   fESDChain->Add(esdfilename);
 100
 101   InitIndex();
 102   InitAlignObjs();
 103 }
 104
 105
 106 //______________________________________________________________________________
 107 AliAlignmentTracks::~AliAlignmentTracks()
 108 {
 109   // Destructor
 110   if (fESDChain) delete fESDChain;
 111
 112   DeleteIndex();
 113   DeleteAlignObjs();
 114
 115   delete fTrackFitter;
 116   delete fMinimizer;
 117
 118   if (fPointsFile) fPointsFile->Close();
 119 }
 120
 121 //______________________________________________________________________________
 122 void AliAlignmentTracks::AddESD(TChain *esdchain)
 123 {
 124   // Add a chain with ESD files
 125   if (fESDChain)
 126     fESDChain->Add(esdchain);
 127   else
 128     fESDChain = esdchain;
 129 }
 130
 131 //______________________________________________________________________________
 132 void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
 133 {
 134   // Add a single file or
 135   // a directory to the chain
 136   // with the ESD files
 137   if (fESDChain)
 138     fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
 139   else {
 140     fESDChain = new TChain(esdtreename);
 141     fESDChain->Add(esdfilename);
 142   }
 143 }
 144
 145 //________________________________________________________________________
 146
 147 void AliAlignmentTracks::ProcessESD(Bool_t onlyITS,
 148                                     Int_t minITSpts,
 149                                     Bool_t cuts,
 150                                     Float_t minMom,Float_t maxMom,
 151                                     Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
 152                                     Float_t minSinTheta,Float_t maxSinTheta)
 153 {
 154   // Analyzes and filters ESD tracks
 155   // Stores the selected track space points
 156   // into the output file
 157
 158   if (!fESDChain) return;
 159
 160   AliESDEvent *esd = new AliESDEvent();
 161   esd->ReadFromTree(fESDChain);
 162   AliESDfriend *esdf = 0;
 163   fESDChain->SetBranchStatus("ESDfriend*",1);
 164   fESDChain->SetBranchAddress("ESDfriend.",&esdf);
 165
 166   // Open the output file
 167   if (fPointsFilename.Data() == "") {
 168     AliWarning("Incorrect output filename!");
 169     return;
 170   }
 171
 172   TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
 173   if (!pointsFile || !pointsFile->IsOpen()) {
 174     AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
 175     return;
 176   }
 177
 178   TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
 179   const AliTrackPointArray *array = 0;
 180   AliTrackPointArray *array2 = 0;
 181   if(onlyITS) {   // only ITS AliTrackPoints
 182     pointsTree->Branch("SP","AliTrackPointArray", &array2);
 183   } else {
 184     pointsTree->Branch("SP","AliTrackPointArray", &array);
 185   }
 186
 187   Int_t ievent = 0;
 188   while (fESDChain->GetEntry(ievent++)) {
 189     if (!esd) break;
 190
 191     esd->SetESDfriend(esdf); //Attach the friend to the ESD
 192
 193     Int_t ntracks = esd->GetNumberOfTracks();
 194     for (Int_t itrack=0; itrack < ntracks; itrack++) {
 195       AliESDtrack * track = esd->GetTrack(itrack);
 196       if (!track) continue;
 197
 198       if(track->GetNcls(0) < minITSpts) continue;
 199       if(cuts) {
 200         if(track->GetP()<minMom || track->GetP()>maxMom) continue;
 201         Float_t abssinphi = TMath::Abs(TMath::Sin(track->GetAlpha()+TMath::ASin(track->GetSnp())));
 202         if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
 203         Float_t sintheta = TMath::Sin(0.5*TMath::Pi()-TMath::ATan(track->GetTgl()));
 204         if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
 205       }
 206       // UInt_t status = AliESDtrack::kITSpid;
 207       // status|=AliESDtrack::kTPCpid;
 208       // status|=AliESDtrack::kTRDpid;
 209       // if ((track->GetStatus() & status) != status) continue;
 210
 211       AliTrackPoint point;
 212       array = track->GetTrackPointArray();
 213
 214       if(onlyITS) {
 215         array2 = new AliTrackPointArray(track->GetNcls(0));
 216         Int_t ipt,volId,modId,layerId;
 217         Int_t jpt=0;
 218         for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 219           array->GetPoint(point,ipt);
 220           volId = point.GetVolumeID();
 221           layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 222           if(layerId<7) {
 223             array2->AddPoint(jpt,&point);
 224             jpt++;
 225           }
 226         }
 227       } // end if(onlyITS)
 228
 229       pointsTree->Fill();
 230     }
 231   }
 232
 233   if (!pointsTree->Write()) {
 234     AliWarning("Can't write the tree with track point arrays!");
 235     return;
 236   }
 237
 238   pointsFile->Close();
 239 }
 240
 241 //______________________________________________________________________________
 242 void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
 243                                    Int_t minITSpts,Float_t maxMatchingAngle,
 244                                    Bool_t cuts,
 245                                    Float_t minMom,Float_t maxMom,
 246                                    Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
 247                                    Float_t minSinTheta,Float_t maxSinTheta)
 248 {
 249   // Analyzes and filters ESD tracks
 250   // Merges inward and outward tracks in one single track
 251   // Stores the selected track space points
 252   // into the output file
 253
 254   if (!fESDChain) return;
 255
 256   AliESDEvent *esd = new AliESDEvent();
 257   esd->ReadFromTree(fESDChain);
 258   AliESDfriend *esdf = 0;
 259   fESDChain->SetBranchStatus("ESDfriend*",1);
 260   fESDChain->SetBranchAddress("ESDfriend.",&esdf);
 261
 262   // Open the output file
 263   if (fPointsFilename.Data() == "") {
 264     AliWarning("Incorrect output filename!");
 265     return;
 266   }
 267
 268   TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
 269   if (!pointsFile || !pointsFile->IsOpen()) {
 270     AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
 271     return;
 272   }
 273
 274   TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
 275   const AliTrackPointArray *array = 0;
 276   AliTrackPointArray *array2 = 0;
 277   pointsTree->Branch("SP","AliTrackPointArray", &array2);
 278
 279   Int_t ievent = 0;
 280   while (fESDChain->GetEntry(ievent++)) {
 281     if (!esd) break;
 282
 283     esd->SetESDfriend(esdf); //Attach the friend to the ESD
 284
 285     Int_t ntracks = esd->GetNumberOfTracks();
 286     if(ntracks<2) continue;
 287     Int_t *goodtracksArray = new Int_t[ntracks];
 288     Float_t *phiArray = new Float_t[ntracks];
 289     Float_t *thetaArray = new Float_t[ntracks];
 290     Int_t ngt=0;
 291     for (Int_t itrack=0; itrack < ntracks; itrack++) {
 292       AliESDtrack * track = esd->GetTrack(itrack);
 293       if (!track) continue;
 294
 295       if(track->GetNcls(0) < minITSpts) continue;
 296       Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
 297       Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
 298       if(cuts) {
 299         if(track->GetP()<minMom || track->GetP()>maxMom) continue;
 300         Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
 301         if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
 302         Float_t sintheta = TMath::Sin(theta);
 303         if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
 304       }
 305       goodtracksArray[ngt]=itrack;
 306       phiArray[ngt]=phi;
 307       thetaArray[ngt]=theta;
 308       ngt++;
 309     }
 310
 311     if(ngt<2) {
 312       delete [] goodtracksArray; goodtracksArray=0;
 313       delete [] phiArray; phiArray=0;
 314       delete [] thetaArray; thetaArray=0;
 315       continue;
 316     }
 317
 318     // check matching of the two tracks from the muon
 319     Float_t min = 10000000.;
 320     Int_t good1 = -1, good2 = -1;
 321     for(Int_t itr1=0; itr1<ngt-1; itr1++) {
 322       for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
 323         Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
 324         if(deltatheta>maxMatchingAngle) continue;
 325         Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
 326         if(deltaphi>maxMatchingAngle) continue;
 327         printf("%f  %f     %f  %f\n",deltaphi,deltatheta,thetaArray[itr1],thetaArray[itr2]);
 328         if(deltatheta+deltaphi<min) {
 329           min=deltatheta+deltaphi;
 330           good1 = goodtracksArray[itr1];
 331           good2 = goodtracksArray[itr2];
 332         }
 333       }
 334     }
 335
 336     delete [] goodtracksArray; goodtracksArray=0;
 337     delete [] phiArray; phiArray=0;
 338     delete [] thetaArray; thetaArray=0;
 339
 340     if(good1<0) continue;
 341
 342     AliESDtrack * track1 = esd->GetTrack(good1);
 343     AliESDtrack * track2 = esd->GetTrack(good2);
 344
 345     Int_t ntotpts;
 346     if(onlyITS) {
 347       ntotpts = track1->GetNcls(0)+track2->GetNcls(0);
 348     } else {
 349       ntotpts = track1->GetTrackPointArray()->GetNPoints()+
 350                 track2->GetTrackPointArray()->GetNPoints();
 351     }
 352     array2 = new AliTrackPointArray(ntotpts);
 353     AliTrackPoint point;
 354     Int_t ipt,volId,modId,layerId;
 355     Int_t jpt=0;
 356     array = track1->GetTrackPointArray();
 357     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 358       array->GetPoint(point,ipt);
 359       if(onlyITS) {
 360         volId = point.GetVolumeID();
 361         layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 362         if(layerId>6) continue;
 363       }
 364       array2->AddPoint(jpt,&point);
 365       jpt++;
 366     }
 367     array = track2->GetTrackPointArray();
 368     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 369       array->GetPoint(point,ipt);
 370       if(onlyITS) {
 371         volId = point.GetVolumeID();
 372         layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 373         if(layerId>6) continue;
 374       }
 375       array2->AddPoint(jpt,&point);
 376       jpt++;
 377     }
 378
 379     pointsTree->Fill();
 380   }
 381
 382   if (!pointsTree->Write()) {
 383     AliWarning("Can't write the tree with track point arrays!");
 384     return;
 385   }
 386
 387   pointsFile->Close();
 388 }
 389
 390 //______________________________________________________________________________
 391 void AliAlignmentTracks::ProcessESD(TSelector *selector)
 392 {
 393   AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
 394 }
 395
 396 //______________________________________________________________________________
 397 void AliAlignmentTracks::BuildIndex()
 398 {
 399   // Build index of points tree entries
 400   // Used for access based on the volume IDs
 401   if (fIsIndexBuilt) return;
 402
 403   fIsIndexBuilt = kTRUE;
 404
 405   // Dummy object is created in order
 406   // to initialize the volume paths
 407   AliAlignObjParams alobj;
 408
 409   fPointsFile = TFile::Open(fPointsFilename);
 410   if (!fPointsFile || !fPointsFile->IsOpen()) {
 411     AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
 412     return;
 413   }
 414
 415   //  AliTrackPointArray* array = new AliTrackPointArray;
 416   AliTrackPointArray* array = 0;
 417   fPointsTree = (TTree*) fPointsFile->Get("spTree");
 418   if (!fPointsTree) {
 419     AliWarning("No pointsTree found!");
 420     return;
 421   }
 422   fPointsTree->SetBranchAddress("SP", &array);
 423
 424   Int_t nArrays = (Int_t)fPointsTree->GetEntries();
 425   for (Int_t iArray = 0; iArray < nArrays; iArray++)
 426     {
 427       fPointsTree->GetEvent(iArray);
 428       if (!array) continue;
 429       for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
 430         UShort_t volId = array->GetVolumeID()[ipoint];
 431         // check if the volId is valid
 432         if (!AliGeomManager::SymName(volId)) {
 433           AliError(Form("The volume id %d has no default volume name !",
 434                         volId));
 435           continue;
 436         }
 437         Int_t modId;
 438         Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId)
 439                       - AliGeomManager::kFirstLayer;
 440         if (!fArrayIndex[layerId][modId]) {
 441           //first entry for this volume
 442           fArrayIndex[layerId][modId] = new TArrayI(1000);
 443         }
 444         else {
 445           Int_t size = fArrayIndex[layerId][modId]->GetSize();
 446           // If needed allocate new size
 447           if (fLastIndex[layerId][modId] >= size)
 448             fArrayIndex[layerId][modId]->Set(size + 1000);
 449         }
 450
 451         // Check if the index is already filled
 452         Bool_t fillIndex = kTRUE;
 453         if (fLastIndex[layerId][modId] != 0) {
 454           if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
 455             fillIndex = kFALSE;
 456         }
 457         // Fill the index array and store last filled index
 458         if (fillIndex) {
 459           (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
 460           fLastIndex[layerId][modId]++;
 461         }
 462       }
 463     }
 464
 465 }
 466
 467 //______________________________________________________________________________
 468 void AliAlignmentTracks::InitIndex()
 469 {
 470   // Initialize the index arrays
 471   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 472   fLastIndex = new Int_t*[nLayers];
 473   fArrayIndex = new TArrayI**[nLayers];
 474   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 475     fLastIndex[iLayer] = new Int_t[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 476     fArrayIndex[iLayer] = new TArrayI*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 477     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 478       fLastIndex[iLayer][iModule] = 0;
 479       fArrayIndex[iLayer][iModule] = 0;
 480     }
 481   }
 482 }
 483
 484 //______________________________________________________________________________
 485 void AliAlignmentTracks::ResetIndex()
 486 {
 487   // Reset the value of the last filled index
 488   // Do not realocate memory
 489
 490   fIsIndexBuilt = kFALSE;
 491
 492   for (Int_t iLayer = 0; iLayer < AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer; iLayer++) {
 493     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 494       fLastIndex[iLayer][iModule] = 0;
 495     }
 496   }
 497 }
 498
 499 //______________________________________________________________________________
 500 void AliAlignmentTracks::DeleteIndex()
 501 {
 502   // Delete the index arrays
 503   // Called by the destructor
 504   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 505     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 506       if (fArrayIndex[iLayer][iModule]) {
 507         delete fArrayIndex[iLayer][iModule];
 508         fArrayIndex[iLayer][iModule] = 0;
 509       }
 510     }
 511     delete [] fLastIndex[iLayer];
 512     delete [] fArrayIndex[iLayer];
 513   }
 514   delete [] fLastIndex;
 515   delete [] fArrayIndex;
 516 }
 517
 518 //______________________________________________________________________________
 519 Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
 520 {
 521   // Read alignment object from a file
 522   // To be replaced by a call to CDB
 523   AliWarning(Form("Method not yet implemented (%s in %s) !",arrayName,alignObjFileName));
 524
 525   return kFALSE;
 526 }
 527
 528 //______________________________________________________________________________
 529 void AliAlignmentTracks::InitAlignObjs()
 530 {
 531   // Initialize the alignment objects array
 532   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 533   fAlignObjs = new AliAlignObj**[nLayers];
 534   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 535     fAlignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 536     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 537       UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
 538       fAlignObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
 539     }
 540   }
 541 }
 542
 543 //______________________________________________________________________________
 544 void AliAlignmentTracks::ResetAlignObjs()
 545 {
 546   // Reset the alignment objects array
 547   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 548     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
 549       fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
 550   }
 551 }
 552
 553 //______________________________________________________________________________
 554 void AliAlignmentTracks::DeleteAlignObjs()
 555 {
 556   // Delete the alignment objects array
 557   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 558     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
 559       if (fAlignObjs[iLayer][iModule])
 560         delete fAlignObjs[iLayer][iModule];
 561     delete [] fAlignObjs[iLayer];
 562   }
 563   delete [] fAlignObjs;
 564   fAlignObjs = 0;
 565 }
 566
 567 void AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
 568                                        AliGeomManager::ELayerID lastLayer,
 569                                        AliGeomManager::ELayerID layerRangeMin,
 570                                        AliGeomManager::ELayerID layerRangeMax,
 571                                        Int_t iterations)
 572 {
 573   // Align detector volumes within
 574   // a given layer range
 575   // (could be whole detector).
 576   // Tracks are fitted only within
 577   // the range defined by the user.
 578   Int_t nModules = 0;
 579   for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++)
 580     nModules += AliGeomManager::LayerSize(iLayer);
 581   TArrayI volIds(nModules);
 582
 583   Int_t modnum = 0;
 584   for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++) {
 585     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer); iModule++) {
 586       UShort_t volId = AliGeomManager::LayerToVolUID(iLayer,iModule);
 587       volIds.AddAt(volId,modnum);
 588       modnum++;
 589     }
 590   }
 591
 592   while (iterations > 0) {
 593     AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
 594     iterations--;
 595   }
 596 }
 597
 598 //______________________________________________________________________________
 599 void AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
 600                                     AliGeomManager::ELayerID layerRangeMin,
 601                                     AliGeomManager::ELayerID layerRangeMax,
 602                                     Int_t iterations)
 603 {
 604   // Align detector volumes within
 605   // a given layer.
 606   // Tracks are fitted only within
 607   // the range defined by the user.
 608   Int_t nModules = AliGeomManager::LayerSize(layer);
 609   TArrayI volIds(nModules);
 610   for (Int_t iModule = 0; iModule < nModules; iModule++) {
 611     UShort_t volId = AliGeomManager::LayerToVolUID(layer,iModule);
 612     volIds.AddAt(volId,iModule);
 613   }
 614
 615   while (iterations > 0) {
 616     AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
 617     iterations--;
 618   }
 619 }
 620
 621 //______________________________________________________________________________
 622 void AliAlignmentTracks::AlignVolume(UShort_t volId, UShort_t volIdFit,
 623                                      Int_t iterations)
 624 {
 625   // Align single detector volume to
 626   // another volume.
 627   // Tracks are fitted only within
 628   // the second volume.
 629   TArrayI volIds(1);
 630   volIds.AddAt(volId,0);
 631   TArrayI volIdsFit(1);
 632   volIdsFit.AddAt(volIdFit,0);
 633
 634   while (iterations > 0) {
 635     AlignVolumes(&volIds,&volIdsFit);
 636     iterations--;
 637   }
 638 }
 639
 640 //______________________________________________________________________________
 641 void AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
 642                                      AliGeomManager::ELayerID layerRangeMin,
 643                                      AliGeomManager::ELayerID layerRangeMax,
 644                                      Int_t iterations)
 645 {
 646   // Align a set of detector volumes.
 647   // Tracks are fitted only within
 648   // the range defined by the user
 649   // (by layerRangeMin and layerRangeMax)
 650   // or within the set of volidsfit
 651   // Repeat the procedure 'iterations' times
 652
 653   Int_t nVolIds = volids->GetSize();
 654   if (nVolIds == 0) {
 655     AliError("Volume IDs array is empty!");
 656     return;
 657   }
 658
 659   // Load only the tracks with at least one
 660   // space point in the set of volume (volids)
 661   BuildIndex();
 662   AliTrackPointArray **points;
 663   // Start the iterations
 664   while (iterations > 0) {
 665     Int_t nArrays = LoadPoints(volids, points);
 666     if (nArrays == 0) return;
 667
 668     AliTrackResiduals *minimizer = CreateMinimizer();
 669     minimizer->SetNTracks(nArrays);
 670     minimizer->InitAlignObj();
 671     AliTrackFitter *fitter = CreateFitter();
 672     for (Int_t iArray = 0; iArray < nArrays; iArray++) {
 673       if (!points[iArray]) continue;
 674       fitter->SetTrackPointArray(points[iArray], kFALSE);
 675       if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
 676       AliTrackPointArray *pVolId,*pTrack;
 677       fitter->GetTrackResiduals(pVolId,pTrack);
 678       minimizer->AddTrackPointArrays(pVolId,pTrack);
 679     }
 680     minimizer->Minimize();
 681
 682     // Update the alignment object(s)
 683     if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
 684       UShort_t volid = (*volids)[iVolId];
 685       Int_t iModule;
 686       AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
 687       AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
 688       *alignObj *= *minimizer->GetAlignObj();
 689       alignObj->Print("");
 690     }
 691
 692     UnloadPoints(nArrays, points);
 693
 694     iterations--;
 695   }
 696 }
 697
 698 //______________________________________________________________________________
 699 Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points)
 700 {
 701   // Load track point arrays with at least
 702   // one space point in a given set of detector
 703   // volumes (array volids).
 704   // Use the already created tree index for
 705   // fast access.
 706
 707   if (!fPointsTree) {
 708     AliError("Tree with the space point arrays not initialized!");
 709     points = 0;
 710     return 0;
 711   }
 712
 713   Int_t nVolIds = volids->GetSize();
 714   if (nVolIds == 0) {
 715     AliError("Volume IDs array is empty!");
 716     points = 0;
 717     return 0;
 718   }
 719
 720   Int_t nArrays = 0;
 721   for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
 722     UShort_t volid = (*volids)[iVolId];
 723     Int_t iModule;
 724     AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
 725
 726     // In case of empty index
 727     if (fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule] == 0) {
 728       AliWarning(Form("There are no space-points belonging to the volume which is to be aligned (Volume ID =%d)!",volid));
 729       continue;
 730     }
 731     nArrays += fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
 732   }
 733
 734   if (nArrays == 0) {
 735     AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
 736     points = 0;
 737     return 0;
 738   }
 739
 740   AliTrackPointArray* array = 0;
 741   fPointsTree->SetBranchAddress("SP", &array);
 742
 743   // Allocate the pointer to the space-point arrays
 744   points = new AliTrackPointArray*[nArrays];
 745   for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
 746
 747   // Init the array used to flag already loaded tree entries
 748   Bool_t *indexUsed = new Bool_t[(UInt_t)fPointsTree->GetEntries()];
 749   for (Int_t i = 0; i < fPointsTree->GetEntries(); i++)
 750     indexUsed[i] = kFALSE;
 751
 752   // Start the loop over the volume ids
 753   Int_t iArray = 0;
 754   for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
 755     UShort_t volid = (*volids)[iVolId];
 756     Int_t iModule;
 757     AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
 758
 759     Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
 760     TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
 761     AliTrackPoint p;
 762
 763     for (Int_t iArrayId = 0; iArrayId < nArraysId; iArrayId++) {
 764
 765       // Get tree entry
 766       Int_t entry = (*index)[iArrayId];
 767       if (indexUsed[entry] == kTRUE) continue;
 768       fPointsTree->GetEvent(entry);
 769       if (!array) {
 770         AliWarning("Wrong space point array index!");
 771         continue;
 772       }
 773       indexUsed[entry] = kTRUE;
 774
 775       // Get the space-point array
 776       Int_t nPoints = array->GetNPoints();
 777       points[iArray] = new AliTrackPointArray(nPoints);
 778       for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
 779         array->GetPoint(p,iPoint);
 780         Int_t modnum;
 781         AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
 782         // check if the layer id is valid
 783         if ((layer < AliGeomManager::kFirstLayer) ||
 784             (layer >= AliGeomManager::kLastLayer)) {
 785           AliError(Form("Layer index is invalid: %d (%d -> %d) !",
 786                         layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
 787           continue;
 788         }
 789         if ((modnum >= AliGeomManager::LayerSize(layer)) ||
 790             (modnum < 0)) {
 791           AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
 792                         layer,modnum,AliGeomManager::LayerSize(layer)));
 793           continue;
 794         }
 795
 796         // Misalignment is introduced here
 797         // Switch it off in case of real
 798         // alignment job!
 799         if (fMisalignObjs) {
 800           AliAlignObj *misalignObj = fMisalignObjs[layer-AliGeomManager::kFirstLayer][modnum];
 801           if (misalignObj)
 802             misalignObj->Transform(p);
 803         }
 804         // End of misalignment
 805
 806         AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
 807         alignObj->Transform(p);
 808         points[iArray]->AddPoint(iPoint,&p);
 809       }
 810       iArray++;
 811     }
 812   }
 813
 814
 815   delete [] indexUsed;
 816
 817   return nArrays;
 818 }
 819
 820 //______________________________________________________________________________
 821 void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
 822 {
 823   // Unload track point arrays for a given
 824   // detector volume
 825   for (Int_t iArray = 0; iArray < n; iArray++)
 826     delete points[iArray];
 827   delete [] points;
 828 }
 829
 830 //______________________________________________________________________________
 831 AliTrackFitter *AliAlignmentTracks::CreateFitter()
 832 {
 833   // Check if the user has already supplied
 834   // a track fitter object.
 835   // If not, create a default one.
 836   if (!fTrackFitter)
 837     fTrackFitter = new AliTrackFitterRieman;
 838
 839   return fTrackFitter;
 840 }
 841
 842 //______________________________________________________________________________
 843 AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
 844 {
 845   // Check if the user has already supplied
 846   // a track residuals minimizer object.
 847   // If not, create a default one.
 848   if (!fMinimizer)
 849     fMinimizer = new AliTrackResidualsChi2;
 850
 851   return fMinimizer;
 852 }
 853
 854 //______________________________________________________________________________
 855 Bool_t AliAlignmentTracks::Misalign(const char *misalignObjFileName, const char* arrayName)
 856 {
 857   // The method reads from a file a set of AliAlignObj which are
 858   // then used to apply misalignments directly on the track
 859   // space-points. The method is supposed to be used only for
 860   // fast development and debugging of the alignment algorithms.
 861   // Be careful not to use it in the case of 'real' alignment
 862   // scenario since it will bias the results.
 863
 864   // Initialize the misalignment objects array
 865   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 866   fMisalignObjs = new AliAlignObj**[nLayers];
 867   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 868     fMisalignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 869     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
 870       fMisalignObjs[iLayer][iModule] = 0x0;
 871   }
 872
 873   // Open the misliagnment file and load the array with
 874   // misalignment objects
 875   TFile* inFile = TFile::Open(misalignObjFileName,"READ");
 876   if (!inFile || !inFile->IsOpen()) {
 877     AliError(Form("Could not open misalignment file %s !",misalignObjFileName));
 878     return kFALSE;
 879   }
 880
 881   TClonesArray* array = ((TClonesArray*) inFile->Get(arrayName));
 882   if (!array) {
 883     AliError(Form("Could not find misalignment array %s in the file %s !",arrayName,misalignObjFileName));
 884     inFile->Close();
 885     return kFALSE;
 886   }
 887   inFile->Close();
 888
 889   // Store the misalignment objects for further usage
 890   Int_t nObjs = array->GetEntriesFast();
 891   AliGeomManager::ELayerID layerId; // volume layer
 892   Int_t modId; // volume ID inside the layer
 893   for(Int_t i=0; i<nObjs; i++)
 894     {
 895       AliAlignObj* alObj = (AliAlignObj*)array->UncheckedAt(i);
 896       alObj->GetVolUID(layerId,modId);
 897       fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
 898     }
 899   return kTRUE;
 900 }