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