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   fTrackFitter(0),
  48   fMinimizer(0)
  49 {
  50   // Default constructor
  51   InitIndex();
  52   InitAlignObjs();
  53 }
  54
  55 //______________________________________________________________________________
  56 AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
  57   fESDChain(esdchain),
  58   fPointsFilename("AliTrackPoints.root"),
  59   fPointsFile(0),
  60   fPointsTree(0),
  61   fLastIndex(0),
  62   fArrayIndex(0),
  63   fIsIndexBuilt(kFALSE),
  64   fTrackFitter(0),
  65   fMinimizer(0)
  66 {
  67   // Constructor in the case
  68   // the user provides an already
  69   // built TChain with ESD trees
  70   InitIndex();
  71   InitAlignObjs();
  72 }
  73
  74
  75 //______________________________________________________________________________
  76 AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
  77   fPointsFilename("AliTrackPoints.root"),
  78   fPointsFile(0),
  79   fPointsTree(0),
  80   fLastIndex(0),
  81   fArrayIndex(0),
  82   fIsIndexBuilt(kFALSE),
  83   fTrackFitter(0),
  84   fMinimizer(0)
  85 {
  86   // Constructor in the case
  87   // the user provides a single ESD file
  88   // or a directory containing ESD files
  89   fESDChain = new TChain(esdtreename);
  90   fESDChain->Add(esdfilename);
  91
  92   InitIndex();
  93   InitAlignObjs();
  94 }
  95
  96 //______________________________________________________________________________
  97 AliAlignmentTracks::AliAlignmentTracks(const AliAlignmentTracks &alignment):
  98   TObject(alignment)
  99 {
 100   // Copy constructor
 101   // not implemented
 102   AliWarning("Copy constructor not implemented!");
 103 }
 104
 105 //______________________________________________________________________________
 106 AliAlignmentTracks& AliAlignmentTracks::operator= (const AliAlignmentTracks& alignment)
 107 {
 108   // Asignment operator
 109   // not implemented
 110   if(this==&alignment) return *this;
 111
 112   AliWarning("Asignment operator not implemented!");
 113
 114   ((TObject *)this)->operator=(alignment);
 115
 116   return *this;
 117 }
 118
 119 //______________________________________________________________________________
 120 AliAlignmentTracks::~AliAlignmentTracks()
 121 {
 122   // Destructor
 123   if (fESDChain) delete fESDChain;
 124
 125   DeleteIndex();
 126   DeleteAlignObjs();
 127
 128   delete fTrackFitter;
 129   delete fMinimizer;
 130
 131   if (fPointsFile) fPointsFile->Close();
 132 }
 133
 134 //______________________________________________________________________________
 135 void AliAlignmentTracks::AddESD(TChain *esdchain)
 136 {
 137   // Add a chain with ESD files
 138   if (fESDChain)
 139     fESDChain->Add(esdchain);
 140   else
 141     fESDChain = esdchain;
 142 }
 143
 144 //______________________________________________________________________________
 145 void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
 146 {
 147   // Add a single file or
 148   // a directory to the chain
 149   // with the ESD files
 150   if (fESDChain)
 151     fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
 152   else {
 153     fESDChain = new TChain(esdtreename);
 154     fESDChain->Add(esdfilename);
 155   }
 156 }
 157
 158 //______________________________________________________________________________
 159 void AliAlignmentTracks::ProcessESD()
 160 {
 161   // Analyzes and filters ESD tracks
 162   // Stores the selected track space points
 163   // into the output file
 164
 165   if (!fESDChain) return;
 166
 167   AliESD *esd = 0;
 168   fESDChain->SetBranchAddress("ESD",&esd);
 169
 170   // Open the output file
 171   if (fPointsFilename.Data() == "") {
 172     AliWarning("Incorrect output filename!");
 173     return;
 174   }
 175
 176   TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
 177   if (!pointsFile || !pointsFile->IsOpen()) {
 178     AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
 179     return;
 180   }
 181
 182   TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
 183   AliTrackPointArray *array = 0;
 184   pointsTree->Branch("SP","AliTrackPointArray", &array);
 185
 186   Int_t ievent = 0;
 187   while (fESDChain->GetEntry(ievent++)) {
 188     if (!esd) break;
 189     Int_t ntracks = esd->GetNumberOfTracks();
 190     for (Int_t itrack=0; itrack < ntracks; itrack++) {
 191       AliESDtrack * track = esd->GetTrack(itrack);
 192       if (!track) continue;
 193
 194       UInt_t status = AliESDtrack::kITSpid;
 195       status|=AliESDtrack::kTPCpid;
 196       status|=AliESDtrack::kTRDpid;
 197       if ((track->GetStatus() & status) != status) continue;
 198
 199       if (track->GetP() < 0.5) continue;
 200
 201       array = track->GetTrackPointArray();
 202       pointsTree->Fill();
 203     }
 204   }
 205
 206   if (!pointsTree->Write()) {
 207     AliWarning("Can't write the tree with track point arrays!");
 208     return;
 209   }
 210
 211   pointsFile->Close();
 212 }
 213
 214 //______________________________________________________________________________
 215 void AliAlignmentTracks::ProcessESD(TSelector *selector)
 216 {
 217   AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
 218 }
 219
 220 //______________________________________________________________________________
 221 void AliAlignmentTracks::BuildIndex()
 222 {
 223   // Build index of points tree entries
 224   // Used for access based on the volume IDs
 225   if (fIsIndexBuilt)
 226     ResetIndex();
 227   else
 228     fIsIndexBuilt = kTRUE;
 229
 230   TFile *fPointsFile = TFile::Open(fPointsFilename);
 231   if (!fPointsFile || !fPointsFile->IsOpen()) {
 232     AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
 233     return;
 234   }
 235
 236   //  AliTrackPointArray* array = new AliTrackPointArray;
 237   AliTrackPointArray* array = 0;
 238   TTree* pointsTree = (TTree*) fPointsFile->Get("spTree");
 239   if (!pointsTree) {
 240     AliWarning("No pointsTree found!");
 241     return;
 242   }
 243   pointsTree->SetBranchAddress("SP", &array);
 244
 245   Int_t nArrays = pointsTree->GetEntries();
 246   for (Int_t iArray = 0; iArray < nArrays; iArray++)
 247     {
 248       pointsTree->GetEvent(iArray);
 249       if (!array) continue;
 250       for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
 251         UShort_t volId = array->GetVolumeID()[ipoint];
 252         Int_t modId;
 253         Int_t layerId = AliAlignObj::VolUIDToLayer(volId,modId)
 254                       - AliAlignObj::kFirstLayer;
 255         if (!fArrayIndex[layerId][modId]) {
 256           //first entry for this volume
 257           fArrayIndex[layerId][modId] = new TArrayI(1000);
 258         }
 259         else {
 260           Int_t size = fArrayIndex[layerId][modId]->GetSize();
 261           // If needed allocate new size
 262           if (fLastIndex[layerId][modId] >= size)
 263             fArrayIndex[layerId][modId]->Set(size + 1000);
 264         }
 265
 266         // Check if the index is already filled
 267         Bool_t fillIndex = kTRUE;
 268         if (fLastIndex[layerId][modId] != 0) {
 269           if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
 270             fillIndex = kFALSE;
 271         }
 272         // Fill the index array and store last filled index
 273         if (fillIndex) {
 274           (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
 275           fLastIndex[layerId][modId]++;
 276         }
 277       }
 278     }
 279
 280 }
 281
 282 // //______________________________________________________________________________
 283 // void AliAlignmentTracks::BuildIndexLayer(AliAlignObj::ELayerID layer)
 284 // {
 285 // }
 286
 287 // //______________________________________________________________________________
 288 // void AliAlignmentTracks::BuildIndexVolume(UShort_t volid)
 289 // {
 290 // }
 291
 292 //______________________________________________________________________________
 293 void AliAlignmentTracks::InitIndex()
 294 {
 295   // Initialize the index arrays
 296   Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
 297   fLastIndex = new Int_t*[nLayers];
 298   fArrayIndex = new TArrayI**[nLayers];
 299   for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
 300     fLastIndex[iLayer] = new Int_t[AliAlignObj::LayerSize(iLayer)];
 301     fArrayIndex[iLayer] = new TArrayI*[AliAlignObj::LayerSize(iLayer)];
 302     for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
 303       fLastIndex[iLayer][iModule] = 0;
 304       fArrayIndex[iLayer][iModule] = 0;
 305     }
 306   }
 307 }
 308
 309 //______________________________________________________________________________
 310 void AliAlignmentTracks::ResetIndex()
 311 {
 312   // Reset the value of the last filled index
 313   // Do not realocate memory
 314   for (Int_t iLayer = AliAlignObj::kFirstLayer; iLayer < AliAlignObj::kLastLayer; iLayer++) {
 315     for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
 316       fLastIndex[iLayer][iModule] = 0;
 317     }
 318   }
 319 }
 320
 321 //______________________________________________________________________________
 322 void AliAlignmentTracks::DeleteIndex()
 323 {
 324   // Delete the index arrays
 325   // Called by the destructor
 326   for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
 327     for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
 328       if (fArrayIndex[iLayer][iModule]) {
 329         delete fArrayIndex[iLayer][iModule];
 330         fArrayIndex[iLayer][iModule] = 0;
 331       }
 332     }
 333     delete [] fLastIndex[iLayer];
 334     delete [] fArrayIndex[iLayer];
 335   }
 336   delete [] fLastIndex;
 337   delete [] fArrayIndex;
 338 }
 339
 340 //______________________________________________________________________________
 341 Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignobjfilename)
 342 {
 343   // Read alignment object from a file
 344   // To be replaced by a call to CDB
 345   AliWarning(Form("Method not yet implemented (%s) !",alignobjfilename));
 346
 347   return kFALSE;
 348 }
 349
 350 //______________________________________________________________________________
 351 void AliAlignmentTracks::InitAlignObjs()
 352 {
 353   // Initialize the alignment objects array
 354   Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
 355   fAlignObjs = new AliAlignObj**[nLayers];
 356   for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
 357     fAlignObjs[iLayer] = new AliAlignObj*[AliAlignObj::LayerSize(iLayer)];
 358     for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
 359       fAlignObjs[iLayer][iModule] = new AliAlignObjAngles;
 360   }
 361 }
 362
 363 //______________________________________________________________________________
 364 void AliAlignmentTracks::ResetAlignObjs()
 365 {
 366   // Reset the alignment objects array
 367   for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
 368     for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
 369       fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
 370   }
 371 }
 372
 373 //______________________________________________________________________________
 374 void AliAlignmentTracks::DeleteAlignObjs()
 375 {
 376   // Delete the alignment objects array
 377   for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
 378     for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
 379       if (fAlignObjs[iLayer][iModule])
 380         delete fAlignObjs[iLayer][iModule];
 381     delete [] fAlignObjs[iLayer];
 382   }
 383   delete [] fAlignObjs;
 384 }
 385
 386 //______________________________________________________________________________
 387 void AliAlignmentTracks::Align(Int_t iterations)
 388 {
 389   // This method is just an example
 390   // how one can user AlignLayer and
 391   // AlignVolume methods to construct
 392   // a custom alignment procedure
 393   while (iterations > 0) {
 394     // First inward pass
 395     AlignLayer(AliAlignObj::kTPC1);
 396     AlignLayer(AliAlignObj::kSSD2);
 397     AlignLayer(AliAlignObj::kSSD1);
 398     AlignLayer(AliAlignObj::kSDD2);
 399     AlignLayer(AliAlignObj::kSDD1);
 400     AlignLayer(AliAlignObj::kSPD2);
 401     AlignLayer(AliAlignObj::kSPD1);
 402     // Outward pass
 403     AlignLayer(AliAlignObj::kSPD2);
 404     AlignLayer(AliAlignObj::kSDD1);
 405     AlignLayer(AliAlignObj::kSDD2);
 406     AlignLayer(AliAlignObj::kSSD1);
 407     AlignLayer(AliAlignObj::kSSD2);
 408     AlignLayer(AliAlignObj::kTPC1);
 409     AlignLayer(AliAlignObj::kTPC2);
 410     AlignLayer(AliAlignObj::kTRD1);
 411     AlignLayer(AliAlignObj::kTRD2);
 412     AlignLayer(AliAlignObj::kTRD3);
 413     AlignLayer(AliAlignObj::kTRD4);
 414     AlignLayer(AliAlignObj::kTRD5);
 415     AlignLayer(AliAlignObj::kTRD6);
 416     AlignLayer(AliAlignObj::kTOF);
 417     // Again inward
 418     AlignLayer(AliAlignObj::kTRD6);
 419     AlignLayer(AliAlignObj::kTRD5);
 420     AlignLayer(AliAlignObj::kTRD4);
 421     AlignLayer(AliAlignObj::kTRD3);
 422     AlignLayer(AliAlignObj::kTRD2);
 423     AlignLayer(AliAlignObj::kTRD1);
 424     AlignLayer(AliAlignObj::kTPC2);
 425   }
 426 }
 427
 428 //______________________________________________________________________________
 429 void AliAlignmentTracks::AlignLayer(AliAlignObj::ELayerID layer,
 430                                     AliAlignObj::ELayerID layerRangeMin,
 431                                     AliAlignObj::ELayerID layerRangeMax,
 432                                     Int_t iterations)
 433 {
 434   // Align detector volumes within
 435   // a given layer.
 436   // Tracks are fitted only within
 437   // the range defined by the user.
 438   while (iterations > 0) {
 439     for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(layer); iModule++) {
 440       UShort_t volId = AliAlignObj::LayerToVolUID(layer,iModule);
 441       AlignVolume(volId,layerRangeMin,layerRangeMax);
 442     }
 443     iterations--;
 444   }
 445 }
 446
 447 //______________________________________________________________________________
 448 void AliAlignmentTracks::AlignVolume(UShort_t volid,
 449                                      AliAlignObj::ELayerID layerRangeMin,
 450                                      AliAlignObj::ELayerID layerRangeMax)
 451 {
 452   // Align a single detector volume.
 453   // Tracks are fitted only within
 454   // the range defined by the user.
 455
 456   // First take the alignment object to be updated
 457   Int_t iModule;
 458   AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
 459   AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
 460
 461   // Then load only the tracks with at least one
 462   // space point in the volume (volid)
 463   BuildIndex();
 464   AliTrackPointArray **points;
 465   Int_t nArrays = LoadPoints(volid, points);
 466   if (nArrays == 0) return;
 467
 468   AliTrackResiduals *minimizer = CreateMinimizer();
 469   minimizer->SetNTracks(nArrays);
 470   minimizer->SetAlignObj(alignObj);
 471   AliTrackFitter *fitter = CreateFitter();
 472   for (Int_t iArray = 0; iArray < nArrays; iArray++) {
 473     fitter->SetTrackPointArray(points[iArray], kFALSE);
 474     AliTrackPointArray *pVolId = 0, *pTrack = 0;
 475     fitter->Fit(volid,pVolId,pTrack,layerRangeMin,layerRangeMax);
 476     minimizer->AddTrackPointArrays(pVolId,pTrack);
 477   }
 478   minimizer->Minimize();
 479
 480   UnloadPoints(nArrays, points);
 481 }
 482
 483 //______________________________________________________________________________
 484 Int_t AliAlignmentTracks::LoadPoints(UShort_t volid, AliTrackPointArray** &points)
 485 {
 486   // Load track point arrays with at least
 487   // one space point in a given detector
 488   // volume (volid).
 489   // Use the already created tree index for
 490   // fast access.
 491   Int_t iModule;
 492   AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
 493   Int_t nArrays = fLastIndex[iLayer][iModule];
 494
 495   // In case of empty index
 496   if (nArrays == 0) {
 497     points = 0;
 498     return 0;
 499   }
 500
 501   if (!fPointsTree) {
 502     AliWarning("Tree with the space point arrays not initialized!");
 503     points = 0;
 504     return 0;
 505   }
 506
 507   AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
 508   TGeoHMatrix m;
 509   alignObj->GetMatrix(m);
 510   Double_t *rot = m.GetRotationMatrix();
 511   Double_t *tr  = m.GetTranslation();
 512
 513   AliTrackPointArray* array = 0;
 514   fPointsTree->SetBranchAddress("SP", &array);
 515
 516   points = new AliTrackPointArray*[nArrays];
 517   TArrayI *index = fArrayIndex[iLayer][iModule];
 518   AliTrackPoint p;
 519   Float_t xyz[3],cov[6];
 520   Float_t newxyz[3];
 521   for (Int_t iArray = 0; iArray < nArrays; iArray++) {
 522     fPointsTree->GetEvent((*index)[iArray]);
 523     if (!array) {
 524       AliWarning("Wrong space point array index!");
 525       continue;
 526     }
 527     Int_t nPoints = array->GetNPoints();
 528     points[iArray] = new AliTrackPointArray(nPoints);
 529     for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
 530       array->GetPoint(p,iPoint);
 531       p.GetXYZ(xyz,cov);
 532       for (Int_t i = 0; i < 3; i++)
 533         newxyz[i] = tr[i]
 534                   + xyz[0]*rot[3*i]
 535                   + xyz[1]*rot[3*i+1]
 536                   + xyz[2]*rot[3*i+2];
 537       p.SetXYZ(newxyz,cov);
 538       points[iArray]->AddPoint(iPoint,&p);
 539     }
 540   }
 541
 542   return nArrays;
 543 }
 544
 545 //______________________________________________________________________________
 546 void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
 547 {
 548   // Unload track point arrays for a given
 549   // detector volume
 550   for (Int_t iArray = 0; iArray < n; iArray++)
 551     delete points[iArray];
 552   delete [] points;
 553 }
 554
 555 //______________________________________________________________________________
 556 AliTrackFitter *AliAlignmentTracks::CreateFitter()
 557 {
 558   // Check if the user has already supplied
 559   // a track fitter object.
 560   // If not, create a default one.
 561   if (!fTrackFitter)
 562     fTrackFitter = new AliTrackFitterRieman;
 563
 564   return fTrackFitter;
 565 }
 566
 567 //______________________________________________________________________________
 568 AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
 569 {
 570   // Check if the user has already supplied
 571   // a track residuals minimizer object.
 572   // If not, create a default one.
 573   if (!fMinimizer)
 574     fMinimizer = new AliTrackResidualsChi2;
 575
 576   return fMinimizer;
 577 }