STEER/AliReconstruction.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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 // class for running the reconstruction                                      //
  21 //                                                                           //
  22 // Clusters and tracks are created for all detectors and all events by       //
  23 // typing:                                                                   //
  24 //                                                                           //
  25 //   AliReconstruction rec;                                                  //
  26 //   rec.Run();                                                              //
  27 //                                                                           //
  28 // The Run method returns kTRUE in case of successful execution.             //
  29 // The name of the galice file can be changed from the default               //
  30 // "galice.root" by                                                          //
  31 //                                                                           //
  32 //   rec.SetGAliceFile("...");                                               //
  33 //                                                                           //
  34 // The reconstruction can be switched on or off for individual detectors by  //
  35 //                                                                           //
  36 //   rec.SetRunReconstruction("...");                                        //
  37 //                                                                           //
  38 // The argument is a (case sensitive) string with the names of the           //
  39 // detectors separated by a space. The special string "ALL" selects all      //
  40 // available detectors. This is the default.                                 //
  41 //                                                                           //
  42 // The tracking in ITS, TPC and TRD and the creation of ESD tracks can be    //
  43 // switched off by                                                           //
  44 //                                                                           //
  45 //   rec.SetRunTracking(kFALSE);                                             //
  46 //                                                                           //
  47 // The filling of additional ESD information can be steered by               //
  48 //                                                                           //
  49 //   rec.SetFillESD("...");                                                  //
  50 //                                                                           //
  51 // Again, the string specifies the list of detectors. The default is "ALL".  //
  52 //                                                                           //
  53 // The reconstruction requires digits as input. For the creation of digits   //
  54 // have a look at the class AliSimulation.                                   //
  55 //                                                                           //
  56 ///////////////////////////////////////////////////////////////////////////////
  57
  58
  59 #include "AliReconstruction.h"
  60 #include "AliRunLoader.h"
  61 #include "AliRun.h"
  62 #include "AliModule.h"
  63 #include "AliDetector.h"
  64 #include "AliTracker.h"
  65 #include "AliESD.h"
  66 #include "AliHeader.h"
  67 #include "AliGenEventHeader.h"
  68 #include "AliESDpid.h"
  69 #include "AliMagF.h"
  70 #include <TArrayF.h>
  71
  72
  73 ClassImp(AliReconstruction)
  74
  75
  76 //_____________________________________________________________________________
  77 AliReconstruction::AliReconstruction(const char* name, const char* title) :
  78   TNamed(name, title)
  79 {
  80 // create reconstruction object with default parameters
  81
  82   Init();
  83 }
  84
  85 //_____________________________________________________________________________
  86 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
  87   TNamed(rec)
  88 {
  89 // copy constructor
  90
  91   fRunReconstruction = rec.fRunReconstruction;
  92   fRunTracking = rec.fRunTracking;
  93   fStopOnError = rec.fStopOnError;
  94
  95   fGAliceFileName = rec.fGAliceFileName;
  96
  97   fRunLoader = NULL;
  98 }
  99
 100 //_____________________________________________________________________________
 101 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
 102 {
 103 // assignment operator
 104
 105   this->~AliReconstruction();
 106   new(this) AliReconstruction(rec);
 107   return *this;
 108 }
 109
 110 //_____________________________________________________________________________
 111 AliReconstruction::~AliReconstruction()
 112 {
 113 // clean up
 114
 115 }
 116
 117 //_____________________________________________________________________________
 118 void AliReconstruction::Init()
 119 {
 120 // set default parameters
 121
 122   fRunReconstruction = "ALL";
 123   fRunTracking = kTRUE;
 124   fFillESD = "ALL";
 125   fStopOnError = kFALSE;
 126
 127   fGAliceFileName = "galice.root";
 128
 129   fRunLoader = NULL;
 130 }
 131
 132
 133 //_____________________________________________________________________________
 134 void AliReconstruction::SetGAliceFile(const char* fileName)
 135 {
 136 // set the name of the galice file
 137
 138   fGAliceFileName = fileName;
 139 }
 140
 141
 142 //_____________________________________________________________________________
 143 Bool_t AliReconstruction::Run()
 144 {
 145 // run the reconstruction
 146
 147   // open the run loader
 148   if (fRunLoader) delete fRunLoader;
 149   fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
 150   if (!fRunLoader) {
 151     Error("Run", "no run loader found in file %s",
 152           fGAliceFileName.Data());
 153     return kFALSE;
 154   }
 155   fRunLoader->LoadgAlice();
 156   gAlice = fRunLoader->GetAliRun();
 157   if (!gAlice) {
 158     Error("Run", "no gAlice object found in file %s",
 159           fGAliceFileName.Data());
 160     return kFALSE;
 161   }
 162
 163   // local reconstruction
 164   if (!fRunReconstruction.IsNull()) {
 165     if (!RunReconstruction(fRunReconstruction)) {
 166       if (fStopOnError) return kFALSE;
 167     }
 168   }
 169   if (!fRunTracking && fFillESD.IsNull()) return kTRUE;
 170
 171   // get loaders and trackers
 172   fITSLoader = fRunLoader->GetLoader("ITSLoader");
 173   if (!fITSLoader) {
 174     Error("Run", "no ITS loader found");
 175     if (fStopOnError) return kFALSE;
 176   }
 177   fITSTracker = NULL;
 178   if (gAlice->GetDetector("ITS")) {
 179     fITSTracker = gAlice->GetDetector("ITS")->CreateTracker();
 180   }
 181   if (!fITSTracker) {
 182     Error("Run", "couldn't create a tracker for ITS");
 183     if (fStopOnError) return kFALSE;
 184   }
 185
 186   fTPCLoader = fRunLoader->GetLoader("TPCLoader");
 187   if (!fTPCLoader) {
 188     Error("Run", "no TPC loader found");
 189     if (fStopOnError) return kFALSE;
 190   }
 191   fTPCTracker = NULL;
 192   if (gAlice->GetDetector("TPC")) {
 193     fTPCTracker = gAlice->GetDetector("TPC")->CreateTracker();
 194   }
 195   if (!fTPCTracker) {
 196     Error("Run", "couldn't create a tracker for TPC");
 197     if (fStopOnError) return kFALSE;
 198   }
 199
 200   fTRDLoader = fRunLoader->GetLoader("TRDLoader");
 201   if (!fTRDLoader) {
 202     Error("Run", "no TRD loader found");
 203     if (fStopOnError) return kFALSE;
 204   }
 205   fTRDTracker = NULL;
 206   if (gAlice->GetDetector("TRD")) {
 207     fTRDTracker = gAlice->GetDetector("TRD")->CreateTracker();
 208   }
 209   if (!fTRDTracker) {
 210     Error("Run", "couldn't create a tracker for TRD");
 211     if (fStopOnError) return kFALSE;
 212   }
 213
 214   fTOFLoader = fRunLoader->GetLoader("TOFLoader");
 215   if (!fTOFLoader) {
 216     Error("Run", "no TOF loader found");
 217     if (fStopOnError) return kFALSE;
 218   }
 219   fTOFTracker = NULL;
 220   if (gAlice->GetDetector("TOF")) {
 221     fTOFTracker = gAlice->GetDetector("TOF")->CreateTracker();
 222   }
 223   if (!fTOFTracker) {
 224     Error("Run", "couldn't create a tracker for TOF");
 225     if (fStopOnError) return kFALSE;
 226   }
 227
 228   // create the ESD output file
 229   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
 230   if (!file->IsOpen()) {
 231     Error("Run", "opening AliESDs.root failed");
 232     if (fStopOnError) return kFALSE;
 233   }
 234
 235   // loop over events
 236   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 237     Info("Run", "processing event %d", iEvent);
 238     AliESD* esd = new AliESD;
 239     fRunLoader->GetEvent(iEvent);
 240     esd->SetRunNumber(gAlice->GetRunNumber());
 241     esd->SetEventNumber(gAlice->GetEvNumber());
 242     esd->SetMagneticField(fRunLoader->GetAliRun()->Field()->SolenoidField());
 243
 244     // barrel tracking
 245     if (fRunTracking) {
 246       if (!RunTracking(esd)) {
 247         if (fStopOnError) return kFALSE;
 248       }
 249     }
 250
 251     // fill ESD
 252     if (!fFillESD.IsNull()) {
 253       if (!FillESD(esd, fFillESD)) {
 254         if (fStopOnError) return kFALSE;
 255       }
 256     }
 257
 258     // combined PID
 259     AliESDpid::MakePID(esd);
 260
 261     // write ESD
 262     char name[100];
 263     sprintf(name, "ESD%d", iEvent);
 264     file->cd();
 265     if (!esd->Write(name)) {
 266       Error("Run", "writing ESD failed");
 267       if (fStopOnError) return kFALSE;
 268     }
 269   }
 270
 271   file->Close();
 272
 273   return kTRUE;
 274 }
 275
 276
 277 //_____________________________________________________________________________
 278 Bool_t AliReconstruction::RunReconstruction(const TString& detectors)
 279 {
 280 // run the reconstruction
 281
 282   TStopwatch stopwatch;
 283   stopwatch.Start();
 284
 285   TString detStr = detectors;
 286   TObjArray* detArray = gAlice->Detectors();
 287   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
 288     AliModule* det = (AliModule*) detArray->At(iDet);
 289     if (!det || !det->IsActive()) continue;
 290     if (IsSelected(det->GetName(), detStr)) {
 291       Info("RunReconstruction", "running reconstruction for %s",
 292            det->GetName());
 293       TStopwatch stopwatchDet;
 294       stopwatchDet.Start();
 295       det->Reconstruct();
 296       Info("RunReconstruction", "execution time for %s:", det->GetName());
 297       stopwatchDet.Print();
 298     }
 299   }
 300
 301   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 302     Error("RunReconstruction", "the following detectors were not found: %s",
 303           detStr.Data());
 304     if (fStopOnError) return kFALSE;
 305   }
 306
 307   Info("RunReconstruction", "execution time:");
 308   stopwatch.Print();
 309
 310   return kTRUE;
 311 }
 312
 313 //_____________________________________________________________________________
 314 Bool_t AliReconstruction::RunTracking(AliESD* esd)
 315 {
 316 // run the barrel tracking
 317
 318   TStopwatch stopwatch;
 319   stopwatch.Start();
 320
 321   // get the primary vertex (from MC for the moment)
 322   TArrayF vertex(3);
 323   fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
 324   Double_t vtxPos[3] = {vertex[0], vertex[1], vertex[2]};
 325   Double_t vtxCov[6] = {
 326     0.005,
 327     0.000, 0.005,
 328     0.000, 0.000, 0.010
 329   };
 330   Double_t vtxErr[3] = {vtxCov[0], vtxCov[2], vtxCov[5]}; // diag. elements
 331   esd->SetVertex(vtxPos, vtxCov);
 332   fITSTracker->SetVertex(vtxPos, vtxErr);
 333   fTPCTracker->SetVertex(vtxPos, vtxErr);
 334   fTRDTracker->SetVertex(vtxPos, vtxErr);
 335
 336   // TPC tracking
 337   Info("RunTracking", "TPC tracking");
 338   fTPCLoader->LoadRecPoints("read");
 339   TTree* tpcTree = fTPCLoader->TreeR();
 340   if (!tpcTree) {
 341     Error("RunTracking", "Can't get the TPC cluster tree");
 342     return kFALSE;
 343   }
 344   fTPCTracker->LoadClusters(tpcTree);
 345   if (fTPCTracker->Clusters2Tracks(esd) != 0) {
 346     Error("RunTracking", "TPC Clusters2Tracks failed");
 347     return kFALSE;
 348   }
 349
 350   gAlice->GetDetector("TPC")->FillESD(esd); // preliminary PID
 351   AliESDpid::MakePID(esd);                  // for the ITS tracker
 352
 353   // ITS tracking
 354   Info("RunTracking", "ITS tracking");
 355   fITSLoader->LoadRecPoints("read");
 356   TTree* itsTree = fITSLoader->TreeR();
 357   if (!itsTree) {
 358     Error("RunTracking", "Can't get the ITS cluster tree");
 359     return kFALSE;
 360   }
 361   fITSTracker->LoadClusters(itsTree);
 362   if (fITSTracker->Clusters2Tracks(esd) != 0) {
 363     Error("RunTracking", "ITS Clusters2Tracks failed");
 364     return kFALSE;
 365   }
 366
 367   // ITS back propagation
 368   Info("RunTracking", "ITS back propagation");
 369   if (fITSTracker->PropagateBack(esd) != 0) {
 370     Error("RunTracking", "ITS backward propagation failed");
 371     return kFALSE;
 372   }
 373
 374   // TPC back propagation
 375   Info("RunTracking", "TPC back propagation");
 376   if (fTPCTracker->PropagateBack(esd) != 0) {
 377     Error("RunTracking", "TPC backward propagation failed");
 378     return kFALSE;
 379   }
 380
 381   // TRD back propagation
 382   Info("RunTracking", "TRD back propagation");
 383   fTRDLoader->LoadRecPoints("read");
 384   TTree* trdTree = fTRDLoader->TreeR();
 385   if (!trdTree) {
 386     Error("RunTracking", "Can't get the TRD cluster tree");
 387     return kFALSE;
 388   }
 389   fTRDTracker->LoadClusters(trdTree);
 390   if (fTRDTracker->PropagateBack(esd) != 0) {
 391     Error("RunTracking", "TRD backward propagation failed");
 392     return kFALSE;
 393   }
 394
 395   // TOF back propagation
 396   Info("RunTracking", "TOF back propagation");
 397   fTOFLoader->LoadDigits("read");
 398   TTree* tofTree = fTOFLoader->TreeD();
 399   if (!tofTree) {
 400     Error("RunTracking", "Can't get the TOF digits tree");
 401     return kFALSE;
 402   }
 403   fTOFTracker->LoadClusters(tofTree);
 404   if (fTOFTracker->PropagateBack(esd) != 0) {
 405     Error("RunTracking", "TOF backward propagation failed");
 406     return kFALSE;
 407   }
 408   fTOFTracker->UnloadClusters();
 409   fTOFLoader->UnloadDigits();
 410
 411   // TRD inward refit
 412   Info("RunTracking", "TRD inward refit");
 413   if (fTRDTracker->RefitInward(esd) != 0) {
 414     Error("RunTracking", "TRD inward refit failed");
 415     return kFALSE;
 416   }
 417   fTRDTracker->UnloadClusters();
 418   fTRDLoader->UnloadRecPoints();
 419
 420   // TPC inward refit
 421   Info("RunTracking", "TPC inward refit");
 422   if (fTPCTracker->RefitInward(esd) != 0) {
 423     Error("RunTracking", "TPC inward refit failed");
 424     return kFALSE;
 425   }
 426   fTPCTracker->UnloadClusters();
 427   fTPCLoader->UnloadRecPoints();
 428
 429   // ITS inward refit
 430   Info("RunTracking", "ITS inward refit");
 431   if (fITSTracker->RefitInward(esd) != 0) {
 432     Error("RunTracking", "ITS inward refit failed");
 433     return kFALSE;
 434   }
 435   fITSTracker->UnloadClusters();
 436   fITSLoader->UnloadRecPoints();
 437
 438   Info("RunTracking", "execution time:");
 439   stopwatch.Print();
 440
 441   return kTRUE;
 442 }
 443
 444 //_____________________________________________________________________________
 445 Bool_t AliReconstruction::FillESD(AliESD* esd, const TString& detectors)
 446 {
 447 // fill the event summary data
 448
 449   TStopwatch stopwatch;
 450   stopwatch.Start();
 451
 452   TString detStr = detectors;
 453   TObjArray* detArray = gAlice->Detectors();
 454   for (Int_t iDet = 0; iDet < detArray->GetEntriesFast(); iDet++) {
 455     AliModule* det = (AliModule*) detArray->At(iDet);
 456     if (!det || !det->IsActive()) continue;
 457     if (IsSelected(det->GetName(), detStr)) {
 458       Info("FillESD", "filling ESD for %s",
 459            det->GetName());
 460       det->FillESD(esd);
 461     }
 462   }
 463
 464   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 465     Error("FillESD", "the following detectors were not found: %s",
 466           detStr.Data());
 467     if (fStopOnError) return kFALSE;
 468   }
 469
 470   Info("FillESD", "execution time:");
 471   stopwatch.Print();
 472
 473   return kTRUE;
 474 }
 475
 476
 477 //_____________________________________________________________________________
 478 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
 479 {
 480 // check whether detName is contained in detectors
 481 // if yes, it is removed from detectors
 482
 483   // check if all detectors are selected
 484   if ((detectors.CompareTo("ALL") == 0) ||
 485       detectors.BeginsWith("ALL ") ||
 486       detectors.EndsWith(" ALL") ||
 487       detectors.Contains(" ALL ")) {
 488     detectors = "ALL";
 489     return kTRUE;
 490   }
 491
 492   // search for the given detector
 493   Bool_t result = kFALSE;
 494   if ((detectors.CompareTo(detName) == 0) ||
 495       detectors.BeginsWith(detName+" ") ||
 496       detectors.EndsWith(" "+detName) ||
 497       detectors.Contains(" "+detName+" ")) {
 498     detectors.ReplaceAll(detName, "");
 499     result = kTRUE;
 500   }
 501
 502   // clean up the detectors string
 503   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
 504   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
 505   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
 506
 507   return result;
 508 }