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