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 "AliAlignObjAngles.h"
  31 #include "AliTrackFitterRieman.h"
  32 #include "AliTrackResidualsChi2.h"
  33 #include "AliESD.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   AliESD *esd = 0;
 161   fESDChain->SetBranchAddress("ESD",&esd);
 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,
 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   AliESD *esd = 0;
 257   fESDChain->SetBranchAddress("ESD",&esd);
 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     Int_t ngt=0;
 287     Int_t goodtracks[10];
 288     for (Int_t itrack=0; itrack < ntracks; itrack++) {
 289       AliESDtrack * track = esd->GetTrack(itrack);
 290       if (!track) continue;
 291
 292       if(track->GetNcls(0) < minITSpts) continue;
 293       if(cuts) {
 294         if(track->GetP()<minMom || track->GetP()>maxMom) continue;
 295         Float_t abssinphi = TMath::Abs(TMath::Sin(track->GetAlpha()+TMath::ASin(track->GetSnp())));
 296         if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
 297         Float_t sintheta = TMath::Sin(0.5*TMath::Pi()-TMath::ATan(track->GetTgl()));
 298         if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
 299       }
 300       if(ngt<10) goodtracks[ngt]=itrack;
 301       ngt++;
 302     }
 303
 304     if(ngt!=2) continue; // this can be changed to check that the two tracks match
 305
 306     AliESDtrack * track1 = esd->GetTrack(goodtracks[0]);
 307     AliESDtrack * track2 = esd->GetTrack(goodtracks[1]);
 308
 309     Int_t ntotpts;
 310     if(onlyITS) {
 311       ntotpts = track1->GetNcls(0)+track2->GetNcls(0);
 312     } else {
 313       ntotpts = track1->GetTrackPointArray()->GetNPoints()+
 314                 track2->GetTrackPointArray()->GetNPoints();
 315     }
 316     array2 = new AliTrackPointArray(ntotpts);
 317     AliTrackPoint point;
 318     Int_t ipt,volId,modId,layerId;
 319     Int_t jpt=0;
 320     array = track1->GetTrackPointArray();
 321     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 322       array->GetPoint(point,ipt);
 323       if(onlyITS) {
 324         volId = point.GetVolumeID();
 325         layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 326         if(layerId>6) continue;
 327       }
 328       array2->AddPoint(jpt,&point);
 329       jpt++;
 330     }
 331     array = track2->GetTrackPointArray();
 332     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 333       array->GetPoint(point,ipt);
 334       if(onlyITS) {
 335         volId = point.GetVolumeID();
 336         layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 337         if(layerId>6) continue;
 338       }
 339       array2->AddPoint(jpt,&point);
 340       jpt++;
 341     }
 342
 343     pointsTree->Fill();
 344   }
 345
 346   if (!pointsTree->Write()) {
 347     AliWarning("Can't write the tree with track point arrays!");
 348     return;
 349   }
 350
 351   pointsFile->Close();
 352 }
 353
 354 //______________________________________________________________________________
 355 void AliAlignmentTracks::ProcessESD(TSelector *selector)
 356 {
 357   AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
 358 }
 359
 360 //______________________________________________________________________________
 361 void AliAlignmentTracks::BuildIndex()
 362 {
 363   // Build index of points tree entries
 364   // Used for access based on the volume IDs
 365   if (fIsIndexBuilt) return;
 366
 367   fIsIndexBuilt = kTRUE;
 368
 369   // Dummy object is created in order
 370   // to initialize the volume paths
 371   AliAlignObjAngles alobj;
 372
 373   fPointsFile = TFile::Open(fPointsFilename);
 374   if (!fPointsFile || !fPointsFile->IsOpen()) {
 375     AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
 376     return;
 377   }
 378
 379   //  AliTrackPointArray* array = new AliTrackPointArray;
 380   AliTrackPointArray* array = 0;
 381   fPointsTree = (TTree*) fPointsFile->Get("spTree");
 382   if (!fPointsTree) {
 383     AliWarning("No pointsTree found!");
 384     return;
 385   }
 386   fPointsTree->SetBranchAddress("SP", &array);
 387
 388   Int_t nArrays = (Int_t)fPointsTree->GetEntries();
 389   for (Int_t iArray = 0; iArray < nArrays; iArray++)
 390     {
 391       fPointsTree->GetEvent(iArray);
 392       if (!array) continue;
 393       for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
 394         UShort_t volId = array->GetVolumeID()[ipoint];
 395         // check if the volId is valid
 396         if (!AliGeomManager::SymName(volId)) {
 397           AliError(Form("The volume id %d has no default volume name !",
 398                         volId));
 399           continue;
 400         }
 401         Int_t modId;
 402         Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId)
 403                       - AliGeomManager::kFirstLayer;
 404         if (!fArrayIndex[layerId][modId]) {
 405           //first entry for this volume
 406           fArrayIndex[layerId][modId] = new TArrayI(1000);
 407         }
 408         else {
 409           Int_t size = fArrayIndex[layerId][modId]->GetSize();
 410           // If needed allocate new size
 411           if (fLastIndex[layerId][modId] >= size)
 412             fArrayIndex[layerId][modId]->Set(size + 1000);
 413         }
 414
 415         // Check if the index is already filled
 416         Bool_t fillIndex = kTRUE;
 417         if (fLastIndex[layerId][modId] != 0) {
 418           if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
 419             fillIndex = kFALSE;
 420         }
 421         // Fill the index array and store last filled index
 422         if (fillIndex) {
 423           (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
 424           fLastIndex[layerId][modId]++;
 425         }
 426       }
 427     }
 428
 429 }
 430
 431 //______________________________________________________________________________
 432 void AliAlignmentTracks::InitIndex()
 433 {
 434   // Initialize the index arrays
 435   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 436   fLastIndex = new Int_t*[nLayers];
 437   fArrayIndex = new TArrayI**[nLayers];
 438   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 439     fLastIndex[iLayer] = new Int_t[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 440     fArrayIndex[iLayer] = new TArrayI*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 441     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 442       fLastIndex[iLayer][iModule] = 0;
 443       fArrayIndex[iLayer][iModule] = 0;
 444     }
 445   }
 446 }
 447
 448 //______________________________________________________________________________
 449 void AliAlignmentTracks::ResetIndex()
 450 {
 451   // Reset the value of the last filled index
 452   // Do not realocate memory
 453
 454   fIsIndexBuilt = kFALSE;
 455
 456   for (Int_t iLayer = 0; iLayer < AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer; iLayer++) {
 457     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 458       fLastIndex[iLayer][iModule] = 0;
 459     }
 460   }
 461 }
 462
 463 //______________________________________________________________________________
 464 void AliAlignmentTracks::DeleteIndex()
 465 {
 466   // Delete the index arrays
 467   // Called by the destructor
 468   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 469     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 470       if (fArrayIndex[iLayer][iModule]) {
 471         delete fArrayIndex[iLayer][iModule];
 472         fArrayIndex[iLayer][iModule] = 0;
 473       }
 474     }
 475     delete [] fLastIndex[iLayer];
 476     delete [] fArrayIndex[iLayer];
 477   }
 478   delete [] fLastIndex;
 479   delete [] fArrayIndex;
 480 }
 481
 482 //______________________________________________________________________________
 483 Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
 484 {
 485   // Read alignment object from a file
 486   // To be replaced by a call to CDB
 487   AliWarning(Form("Method not yet implemented (%s in %s) !",arrayName,alignObjFileName));
 488
 489   return kFALSE;
 490 }
 491
 492 //______________________________________________________________________________
 493 void AliAlignmentTracks::InitAlignObjs()
 494 {
 495   // Initialize the alignment objects array
 496   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 497   fAlignObjs = new AliAlignObj**[nLayers];
 498   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 499     fAlignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 500     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
 501       UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
 502       fAlignObjs[iLayer][iModule] = new AliAlignObjAngles(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
 503     }
 504   }
 505 }
 506
 507 //______________________________________________________________________________
 508 void AliAlignmentTracks::ResetAlignObjs()
 509 {
 510   // Reset the alignment objects array
 511   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 512     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
 513       fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
 514   }
 515 }
 516
 517 //______________________________________________________________________________
 518 void AliAlignmentTracks::DeleteAlignObjs()
 519 {
 520   // Delete the alignment objects array
 521   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 522     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
 523       if (fAlignObjs[iLayer][iModule])
 524         delete fAlignObjs[iLayer][iModule];
 525     delete [] fAlignObjs[iLayer];
 526   }
 527   delete [] fAlignObjs;
 528   fAlignObjs = 0;
 529 }
 530
 531 void AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
 532                                        AliGeomManager::ELayerID lastLayer,
 533                                        AliGeomManager::ELayerID layerRangeMin,
 534                                        AliGeomManager::ELayerID layerRangeMax,
 535                                        Int_t iterations)
 536 {
 537   // Align detector volumes within
 538   // a given layer range
 539   // (could be whole detector).
 540   // Tracks are fitted only within
 541   // the range defined by the user.
 542   Int_t nModules = 0;
 543   for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++)
 544     nModules += AliGeomManager::LayerSize(iLayer);
 545   TArrayI volIds(nModules);
 546
 547   Int_t modnum = 0;
 548   for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++) {
 549     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer); iModule++) {
 550       UShort_t volId = AliGeomManager::LayerToVolUID(iLayer,iModule);
 551       volIds.AddAt(volId,modnum);
 552       modnum++;
 553     }
 554   }
 555
 556   while (iterations > 0) {
 557     AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
 558     iterations--;
 559   }
 560 }
 561
 562 //______________________________________________________________________________
 563 void AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
 564                                     AliGeomManager::ELayerID layerRangeMin,
 565                                     AliGeomManager::ELayerID layerRangeMax,
 566                                     Int_t iterations)
 567 {
 568   // Align detector volumes within
 569   // a given layer.
 570   // Tracks are fitted only within
 571   // the range defined by the user.
 572   Int_t nModules = AliGeomManager::LayerSize(layer);
 573   TArrayI volIds(nModules);
 574   for (Int_t iModule = 0; iModule < nModules; iModule++) {
 575     UShort_t volId = AliGeomManager::LayerToVolUID(layer,iModule);
 576     volIds.AddAt(volId,iModule);
 577   }
 578
 579   while (iterations > 0) {
 580     AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
 581     iterations--;
 582   }
 583 }
 584
 585 //______________________________________________________________________________
 586 void AliAlignmentTracks::AlignVolume(UShort_t volId, UShort_t volIdFit,
 587                                      Int_t iterations)
 588 {
 589   // Align single detector volume to
 590   // another volume.
 591   // Tracks are fitted only within
 592   // the second volume.
 593   TArrayI volIds(1);
 594   volIds.AddAt(volId,0);
 595   TArrayI volIdsFit(1);
 596   volIdsFit.AddAt(volIdFit,0);
 597
 598   while (iterations > 0) {
 599     AlignVolumes(&volIds,&volIdsFit);
 600     iterations--;
 601   }
 602 }
 603
 604 //______________________________________________________________________________
 605 void AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
 606                                      AliGeomManager::ELayerID layerRangeMin,
 607                                      AliGeomManager::ELayerID layerRangeMax,
 608                                      Int_t iterations)
 609 {
 610   // Align a set of detector volumes.
 611   // Tracks are fitted only within
 612   // the range defined by the user
 613   // (by layerRangeMin and layerRangeMax)
 614   // or within the set of volidsfit
 615   // Repeat the procedure 'iterations' times
 616
 617   Int_t nVolIds = volids->GetSize();
 618   if (nVolIds == 0) {
 619     AliError("Volume IDs array is empty!");
 620     return;
 621   }
 622
 623   // Load only the tracks with at least one
 624   // space point in the set of volume (volids)
 625   BuildIndex();
 626   AliTrackPointArray **points;
 627   // Start the iterations
 628   while (iterations > 0) {
 629     Int_t nArrays = LoadPoints(volids, points);
 630     if (nArrays == 0) return;
 631
 632     AliTrackResiduals *minimizer = CreateMinimizer();
 633     minimizer->SetNTracks(nArrays);
 634     minimizer->InitAlignObj();
 635     AliTrackFitter *fitter = CreateFitter();
 636     for (Int_t iArray = 0; iArray < nArrays; iArray++) {
 637       if (!points[iArray]) continue;
 638       fitter->SetTrackPointArray(points[iArray], kFALSE);
 639       if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
 640       AliTrackPointArray *pVolId,*pTrack;
 641       fitter->GetTrackResiduals(pVolId,pTrack);
 642       minimizer->AddTrackPointArrays(pVolId,pTrack);
 643     }
 644     minimizer->Minimize();
 645
 646     // Update the alignment object(s)
 647     if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
 648       UShort_t volid = (*volids)[iVolId];
 649       Int_t iModule;
 650       AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
 651       AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
 652       *alignObj *= *minimizer->GetAlignObj();
 653       alignObj->Print("");
 654     }
 655
 656     UnloadPoints(nArrays, points);
 657
 658     iterations--;
 659   }
 660 }
 661
 662 //______________________________________________________________________________
 663 Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points)
 664 {
 665   // Load track point arrays with at least
 666   // one space point in a given set of detector
 667   // volumes (array volids).
 668   // Use the already created tree index for
 669   // fast access.
 670
 671   if (!fPointsTree) {
 672     AliError("Tree with the space point arrays not initialized!");
 673     points = 0;
 674     return 0;
 675   }
 676
 677   Int_t nVolIds = volids->GetSize();
 678   if (nVolIds == 0) {
 679     AliError("Volume IDs array is empty!");
 680     points = 0;
 681     return 0;
 682   }
 683
 684   Int_t nArrays = 0;
 685   for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
 686     UShort_t volid = (*volids)[iVolId];
 687     Int_t iModule;
 688     AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
 689
 690     // In case of empty index
 691     if (fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule] == 0) {
 692       AliWarning(Form("There are no space-points belonging to the volume which is to be aligned (Volume ID =%d)!",volid));
 693       continue;
 694     }
 695     nArrays += fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
 696   }
 697
 698   if (nArrays == 0) {
 699     AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
 700     points = 0;
 701     return 0;
 702   }
 703
 704   AliTrackPointArray* array = 0;
 705   fPointsTree->SetBranchAddress("SP", &array);
 706
 707   // Allocate the pointer to the space-point arrays
 708   points = new AliTrackPointArray*[nArrays];
 709   for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
 710
 711   // Init the array used to flag already loaded tree entries
 712   Bool_t *indexUsed = new Bool_t[(UInt_t)fPointsTree->GetEntries()];
 713   for (Int_t i = 0; i < fPointsTree->GetEntries(); i++)
 714     indexUsed[i] = kFALSE;
 715
 716   // Start the loop over the volume ids
 717   Int_t iArray = 0;
 718   for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
 719     UShort_t volid = (*volids)[iVolId];
 720     Int_t iModule;
 721     AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
 722
 723     Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
 724     TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
 725     AliTrackPoint p;
 726
 727     for (Int_t iArrayId = 0; iArrayId < nArraysId; iArrayId++) {
 728
 729       // Get tree entry
 730       Int_t entry = (*index)[iArrayId];
 731       if (indexUsed[entry] == kTRUE) continue;
 732       fPointsTree->GetEvent(entry);
 733       if (!array) {
 734         AliWarning("Wrong space point array index!");
 735         continue;
 736       }
 737       indexUsed[entry] = kTRUE;
 738
 739       // Get the space-point array
 740       Int_t nPoints = array->GetNPoints();
 741       points[iArray] = new AliTrackPointArray(nPoints);
 742       for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
 743         array->GetPoint(p,iPoint);
 744         Int_t modnum;
 745         AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
 746         // check if the layer id is valid
 747         if ((layer < AliGeomManager::kFirstLayer) ||
 748             (layer >= AliGeomManager::kLastLayer)) {
 749           AliError(Form("Layer index is invalid: %d (%d -> %d) !",
 750                         layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
 751           continue;
 752         }
 753         if ((modnum >= AliGeomManager::LayerSize(layer)) ||
 754             (modnum < 0)) {
 755           AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
 756                         layer,modnum,AliGeomManager::LayerSize(layer)));
 757           continue;
 758         }
 759
 760         // Misalignment is introduced here
 761         // Switch it off in case of real
 762         // alignment job!
 763         if (fMisalignObjs) {
 764           AliAlignObj *misalignObj = fMisalignObjs[layer-AliGeomManager::kFirstLayer][modnum];
 765           if (misalignObj)
 766             misalignObj->Transform(p);
 767         }
 768         // End of misalignment
 769
 770         AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
 771         alignObj->Transform(p);
 772         points[iArray]->AddPoint(iPoint,&p);
 773       }
 774       iArray++;
 775     }
 776   }
 777
 778
 779   delete [] indexUsed;
 780
 781   return nArrays;
 782 }
 783
 784 //______________________________________________________________________________
 785 void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
 786 {
 787   // Unload track point arrays for a given
 788   // detector volume
 789   for (Int_t iArray = 0; iArray < n; iArray++)
 790     delete points[iArray];
 791   delete [] points;
 792 }
 793
 794 //______________________________________________________________________________
 795 AliTrackFitter *AliAlignmentTracks::CreateFitter()
 796 {
 797   // Check if the user has already supplied
 798   // a track fitter object.
 799   // If not, create a default one.
 800   if (!fTrackFitter)
 801     fTrackFitter = new AliTrackFitterRieman;
 802
 803   return fTrackFitter;
 804 }
 805
 806 //______________________________________________________________________________
 807 AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
 808 {
 809   // Check if the user has already supplied
 810   // a track residuals minimizer object.
 811   // If not, create a default one.
 812   if (!fMinimizer)
 813     fMinimizer = new AliTrackResidualsChi2;
 814
 815   return fMinimizer;
 816 }
 817
 818 //______________________________________________________________________________
 819 Bool_t AliAlignmentTracks::Misalign(const char *misalignObjFileName, const char* arrayName)
 820 {
 821   // The method reads from a file a set of AliAlignObj which are
 822   // then used to apply misalignments directly on the track
 823   // space-points. The method is supposed to be used only for
 824   // fast development and debugging of the alignment algorithms.
 825   // Be careful not to use it in the case of 'real' alignment
 826   // scenario since it will bias the results.
 827
 828   // Initialize the misalignment objects array
 829   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
 830   fMisalignObjs = new AliAlignObj**[nLayers];
 831   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
 832     fMisalignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
 833     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
 834       fMisalignObjs[iLayer][iModule] = 0x0;
 835   }
 836
 837   // Open the misliagnment file and load the array with
 838   // misalignment objects
 839   TFile* inFile = TFile::Open(misalignObjFileName,"READ");
 840   if (!inFile || !inFile->IsOpen()) {
 841     AliError(Form("Could not open misalignment file %s !",misalignObjFileName));
 842     return kFALSE;
 843   }
 844
 845   TClonesArray* array = ((TClonesArray*) inFile->Get(arrayName));
 846   if (!array) {
 847     AliError(Form("Could not find misalignment array %s in the file %s !",arrayName,misalignObjFileName));
 848     inFile->Close();
 849     return kFALSE;
 850   }
 851   inFile->Close();
 852
 853   // Store the misalignment objects for further usage
 854   Int_t nObjs = array->GetEntriesFast();
 855   AliGeomManager::ELayerID layerId; // volume layer
 856   Int_t modId; // volume ID inside the layer
 857   for(Int_t i=0; i<nObjs; i++)
 858     {
 859       AliAlignObj* alObj = (AliAlignObj*)array->UncheckedAt(i);
 860       alObj->GetVolUID(layerId,modId);
 861       fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
 862     }
 863   return kTRUE;
 864 }