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