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 //                                                                           //
  30 // If the input to the reconstruction are not simulated digits but raw data, //
  31 // this can be specified by an argument of the Run method or by the method   //
  32 //                                                                           //
  33 //   rec.SetInput("...");                                                    //
  34 //                                                                           //
  35 // The input formats and the corresponding argument are:                     //
  36 // - DDL raw data files: directory name, ends with "/"                       //
  37 // - raw data root file: root file name, extension ".root"                   //
  38 // - raw data DATE file: DATE file name, any other non-empty string          //
  39 // - MC root files     : empty string, default                               //
  40 //                                                                           //
  41 // By default all events are reconstructed. The reconstruction can be        //
  42 // limited to a range of events by giving the index of the first and the     //
  43 // last event as an argument to the Run method or by calling                 //
  44 //                                                                           //
  45 //   rec.SetEventRange(..., ...);                                            //
  46 //                                                                           //
  47 // The index -1 (default) can be used for the last event to indicate no      //
  48 // upper limit of the event range.                                           //
  49 //                                                                           //
  50 // The name of the galice file can be changed from the default               //
  51 // "galice.root" by passing it as argument to the AliReconstruction          //
  52 // constructor or by                                                         //
  53 //                                                                           //
  54 //   rec.SetGAliceFile("...");                                               //
  55 //                                                                           //
  56 // The local reconstruction can be switched on or off for individual         //
  57 // detectors by                                                              //
  58 //                                                                           //
  59 //   rec.SetRunLocalReconstruction("...");                                   //
  60 //                                                                           //
  61 // The argument is a (case sensitive) string with the names of the           //
  62 // detectors separated by a space. The special string "ALL" selects all      //
  63 // available detectors. This is the default.                                 //
  64 //                                                                           //
  65 // The reconstruction of the primary vertex position can be switched off by  //
  66 //                                                                           //
  67 //   rec.SetRunVertexFinder(kFALSE);                                         //
  68 //                                                                           //
  69 // The tracking and the creation of ESD tracks can be switched on for        //
  70 // selected detectors by                                                     //
  71 //                                                                           //
  72 //   rec.SetRunTracking("...");                                              //
  73 //                                                                           //
  74 // Uniform/nonuniform field tracking switches (default: uniform field)       //
  75 //                                                                           //
  76 //   rec.SetUniformFieldTracking();  ( rec.SetNonuniformFieldTracking(); )   //
  77 //                                                                           //
  78 // The filling of additional ESD information can be steered by               //
  79 //                                                                           //
  80 //   rec.SetFillESD("...");                                                  //
  81 //                                                                           //
  82 // Again, for both methods the string specifies the list of detectors.       //
  83 // The default is "ALL".                                                     //
  84 //                                                                           //
  85 // The call of the shortcut method                                           //
  86 //                                                                           //
  87 //   rec.SetRunReconstruction("...");                                        //
  88 //                                                                           //
  89 // is equivalent to calling SetRunLocalReconstruction, SetRunTracking and    //
  90 // SetFillESD with the same detector selecting string as argument.           //
  91 //                                                                           //
  92 // The reconstruction requires digits or raw data as input. For the creation //
  93 // of digits and raw data have a look at the class AliSimulation.            //
  94 //                                                                           //
  95 // For debug purposes the method SetCheckPointLevel can be used. If the      //
  96 // argument is greater than 0, files with ESD events will be written after   //
  97 // selected steps of the reconstruction for each event:                      //
  98 //   level 1: after tracking and after filling of ESD (final)                //
  99 //   level 2: in addition after each tracking step                           //
 100 //   level 3: in addition after the filling of ESD for each detector         //
 101 // If a final check point file exists for an event, this event will be       //
 102 // skipped in the reconstruction. The tracking and the filling of ESD for    //
 103 // a detector will be skipped as well, if the corresponding check point      //
 104 // file exists. The ESD event will then be loaded from the file instead.     //
 105 //                                                                           //
 106 ///////////////////////////////////////////////////////////////////////////////
 107
 108 #include <TArrayF.h>
 109 #include <TFile.h>
 110 #include <TSystem.h>
 111 #include <TROOT.h>
 112 #include <TPluginManager.h>
 113 #include <TStopwatch.h>
 114
 115 #include "AliReconstruction.h"
 116 #include "AliReconstructor.h"
 117 #include "AliLog.h"
 118 #include "AliRunLoader.h"
 119 #include "AliRun.h"
 120 #include "AliRawReaderFile.h"
 121 #include "AliRawReaderDate.h"
 122 #include "AliRawReaderRoot.h"
 123 #include "AliESD.h"
 124 #include "AliESDVertex.h"
 125 #include "AliTracker.h"
 126 #include "AliVertexer.h"
 127 #include "AliHeader.h"
 128 #include "AliGenEventHeader.h"
 129 #include "AliPID.h"
 130 #include "AliESDpid.h"
 131 #include "AliMagF.h"
 132
 133 ClassImp(AliReconstruction)
 134
 135
 136 //_____________________________________________________________________________
 137 const char* AliReconstruction::fgkDetectorName[AliReconstruction::fgkNDetectors] = {"ITS", "TPC", "TRD", "TOF", "PHOS", "RICH", "EMCAL", "MUON", "FMD", "ZDC", "PMD", "START", "VZERO", "CRT", "HLT"};
 138
 139 //_____________________________________________________________________________
 140 AliReconstruction::AliReconstruction(const char* gAliceFilename,
 141                                      const char* name, const char* title) :
 142   TNamed(name, title),
 143
 144   fRunLocalReconstruction("ALL"),
 145   fUniformField(kTRUE),
 146   fRunVertexFinder(kTRUE),
 147   fRunHLTTracking(kFALSE),
 148   fRunTracking("ALL"),
 149   fFillESD("ALL"),
 150   fGAliceFileName(gAliceFilename),
 151   fInput(""),
 152   fFirstEvent(0),
 153   fLastEvent(-1),
 154   fStopOnError(kFALSE),
 155   fCheckPointLevel(0),
 156   fOptions(),
 157
 158   fRunLoader(NULL),
 159   fRawReader(NULL),
 160
 161   fVertexer(NULL)
 162 {
 163 // create reconstruction object with default parameters
 164
 165   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 166     fReconstructor[iDet] = NULL;
 167     fLoader[iDet] = NULL;
 168     fTracker[iDet] = NULL;
 169   }
 170   AliPID pid;
 171 }
 172
 173 //_____________________________________________________________________________
 174 AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
 175   TNamed(rec),
 176
 177   fRunLocalReconstruction(rec.fRunLocalReconstruction),
 178   fUniformField(rec.fUniformField),
 179   fRunVertexFinder(rec.fRunVertexFinder),
 180   fRunHLTTracking(rec.fRunHLTTracking),
 181   fRunTracking(rec.fRunTracking),
 182   fFillESD(rec.fFillESD),
 183   fGAliceFileName(rec.fGAliceFileName),
 184   fInput(rec.fInput),
 185   fFirstEvent(rec.fFirstEvent),
 186   fLastEvent(rec.fLastEvent),
 187   fStopOnError(rec.fStopOnError),
 188   fCheckPointLevel(0),
 189   fOptions(),
 190
 191   fRunLoader(NULL),
 192   fRawReader(NULL),
 193
 194   fVertexer(NULL)
 195 {
 196 // copy constructor
 197
 198   for (Int_t i = 0; i < fOptions.GetEntriesFast(); i++) {
 199     if (rec.fOptions[i]) fOptions.Add(rec.fOptions[i]->Clone());
 200   }
 201   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 202     fReconstructor[iDet] = NULL;
 203     fLoader[iDet] = NULL;
 204     fTracker[iDet] = NULL;
 205   }
 206 }
 207
 208 //_____________________________________________________________________________
 209 AliReconstruction& AliReconstruction::operator = (const AliReconstruction& rec)
 210 {
 211 // assignment operator
 212
 213   this->~AliReconstruction();
 214   new(this) AliReconstruction(rec);
 215   return *this;
 216 }
 217
 218 //_____________________________________________________________________________
 219 AliReconstruction::~AliReconstruction()
 220 {
 221 // clean up
 222
 223   CleanUp();
 224   fOptions.Delete();
 225 }
 226
 227
 228 //_____________________________________________________________________________
 229 void AliReconstruction::SetGAliceFile(const char* fileName)
 230 {
 231 // set the name of the galice file
 232
 233   fGAliceFileName = fileName;
 234 }
 235
 236 //_____________________________________________________________________________
 237 void AliReconstruction::SetOption(const char* detector, const char* option)
 238 {
 239 // set options for the reconstruction of a detector
 240
 241   TObject* obj = fOptions.FindObject(detector);
 242   if (obj) fOptions.Remove(obj);
 243   fOptions.Add(new TNamed(detector, option));
 244 }
 245
 246
 247 //_____________________________________________________________________________
 248 Bool_t AliReconstruction::Run(const char* input,
 249                               Int_t firstEvent, Int_t lastEvent)
 250 {
 251 // run the reconstruction
 252
 253   // set the input
 254   if (!input) input = fInput.Data();
 255   TString fileName(input);
 256   if (fileName.EndsWith("/")) {
 257     fRawReader = new AliRawReaderFile(fileName);
 258   } else if (fileName.EndsWith(".root")) {
 259     fRawReader = new AliRawReaderRoot(fileName);
 260   } else if (!fileName.IsNull()) {
 261     fRawReader = new AliRawReaderDate(fileName);
 262     fRawReader->SelectEvents(7);
 263   }
 264
 265   // get the run loader
 266   if (!InitRunLoader()) return kFALSE;
 267
 268   // local reconstruction
 269   if (!fRunLocalReconstruction.IsNull()) {
 270     if (!RunLocalReconstruction(fRunLocalReconstruction)) {
 271       if (fStopOnError) {CleanUp(); return kFALSE;}
 272     }
 273   }
 274 //  if (!fRunVertexFinder && fRunTracking.IsNull() &&
 275 //      fFillESD.IsNull()) return kTRUE;
 276
 277   // get vertexer
 278   if (fRunVertexFinder && !CreateVertexer()) {
 279     if (fStopOnError) {
 280       CleanUp();
 281       return kFALSE;
 282     }
 283   }
 284
 285   // get trackers
 286   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
 287     if (fStopOnError) {
 288       CleanUp();
 289       return kFALSE;
 290     }
 291   }
 292
 293   // get the possibly already existing ESD file and tree
 294   AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
 295   TFile* fileOld = NULL;
 296   TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
 297   if (!gSystem->AccessPathName("AliESDs.root")){
 298     gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
 299     fileOld = TFile::Open("AliESDs.old.root");
 300     if (fileOld && fileOld->IsOpen()) {
 301       treeOld = (TTree*) fileOld->Get("esdTree");
 302       if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
 303       hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
 304       if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
 305     }
 306   }
 307
 308   // create the ESD output file and tree
 309   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
 310   if (!file->IsOpen()) {
 311     AliError("opening AliESDs.root failed");
 312     if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 313   }
 314   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
 315   tree->Branch("ESD", "AliESD", &esd);
 316   TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
 317   hlttree->Branch("ESD", "AliESD", &hltesd);
 318   delete esd; delete hltesd;
 319   esd = NULL; hltesd = NULL;
 320   gROOT->cd();
 321
 322   // loop over events
 323   if (fRawReader) fRawReader->RewindEvents();
 324
 325   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 326     if (fRawReader) fRawReader->NextEvent();
 327     if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
 328       // copy old ESD to the new one
 329       if (treeOld) {
 330         treeOld->SetBranchAddress("ESD", &esd);
 331         treeOld->GetEntry(iEvent);
 332       }
 333       tree->Fill();
 334       if (hlttreeOld) {
 335         hlttreeOld->SetBranchAddress("ESD", &hltesd);
 336         hlttreeOld->GetEntry(iEvent);
 337       }
 338       hlttree->Fill();
 339       continue;
 340     }
 341
 342     AliInfo(Form("processing event %d", iEvent));
 343     fRunLoader->GetEvent(iEvent);
 344
 345     char fileName[256];
 346     sprintf(fileName, "ESD_%d.%d_final.root",
 347             fRunLoader->GetHeader()->GetRun(),
 348             fRunLoader->GetHeader()->GetEventNrInRun());
 349     if (!gSystem->AccessPathName(fileName)) continue;
 350
 351     // local reconstruction
 352     if (!fRunLocalReconstruction.IsNull()) {
 353       if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
 354         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 355       }
 356     }
 357
 358     esd = new AliESD; hltesd = new AliESD;
 359     esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 360     hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
 361     esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
 362     hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
 363     if (gAlice) {
 364       esd->SetMagneticField(gAlice->Field()->SolenoidField());
 365       hltesd->SetMagneticField(gAlice->Field()->SolenoidField());
 366     } else {
 367       // ???
 368     }
 369
 370     // vertex finder
 371     if (fRunVertexFinder) {
 372       if (!ReadESD(esd, "vertex")) {
 373         if (!RunVertexFinder(esd)) {
 374           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 375         }
 376         if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
 377       }
 378     }
 379
 380     // HLT tracking
 381     if (!fRunTracking.IsNull()) {
 382       if (fRunHLTTracking) {
 383         hltesd->SetVertex(esd->GetVertex());
 384         if (!RunHLTTracking(hltesd)) {
 385           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 386         }
 387       }
 388     }
 389
 390     // barrel tracking
 391     if (!fRunTracking.IsNull()) {
 392       if (!ReadESD(esd, "tracking")) {
 393         if (!RunTracking(esd)) {
 394           if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 395         }
 396         if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
 397       }
 398     }
 399
 400     // fill ESD
 401     if (!fFillESD.IsNull()) {
 402       if (!FillESD(esd, fFillESD)) {
 403         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
 404       }
 405     }
 406
 407     // combined PID
 408     AliESDpid::MakePID(esd);
 409     if (fCheckPointLevel > 1) WriteESD(esd, "PID");
 410
 411     // write ESD
 412     tree->Fill();
 413     // write HLT ESD
 414     hlttree->Fill();
 415
 416     if (fCheckPointLevel > 0) WriteESD(esd, "final");
 417     delete esd; delete hltesd;
 418     esd = NULL; hltesd = NULL;
 419   }
 420
 421   file->cd();
 422   tree->Write();
 423   hlttree->Write();
 424   CleanUp(file, fileOld);
 425
 426   return kTRUE;
 427 }
 428
 429
 430 //_____________________________________________________________________________
 431 Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
 432 {
 433 // run the local reconstruction
 434
 435   TStopwatch stopwatch;
 436   stopwatch.Start();
 437
 438   TString detStr = detectors;
 439   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 440     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 441     AliReconstructor* reconstructor = GetReconstructor(iDet);
 442     if (!reconstructor) continue;
 443     if (reconstructor->HasLocalReconstruction()) continue;
 444
 445     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 446     TStopwatch stopwatchDet;
 447     stopwatchDet.Start();
 448     if (fRawReader) {
 449       fRawReader->RewindEvents();
 450       reconstructor->Reconstruct(fRunLoader, fRawReader);
 451     } else {
 452       reconstructor->Reconstruct(fRunLoader);
 453     }
 454     AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 455                  fgkDetectorName[iDet],
 456                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 457   }
 458
 459   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 460     AliError(Form("the following detectors were not found: %s",
 461                   detStr.Data()));
 462     if (fStopOnError) return kFALSE;
 463   }
 464
 465   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 466                stopwatch.RealTime(),stopwatch.CpuTime()));
 467
 468   return kTRUE;
 469 }
 470
 471 //_____________________________________________________________________________
 472 Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
 473 {
 474 // run the local reconstruction
 475
 476   TStopwatch stopwatch;
 477   stopwatch.Start();
 478
 479   TString detStr = detectors;
 480   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 481     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 482     AliReconstructor* reconstructor = GetReconstructor(iDet);
 483     if (!reconstructor) continue;
 484     AliLoader* loader = fLoader[iDet];
 485
 486     // conversion of digits
 487     if (fRawReader && reconstructor->HasDigitConversion()) {
 488       AliInfo(Form("converting raw data digits into root objects for %s",
 489                    fgkDetectorName[iDet]));
 490       TStopwatch stopwatchDet;
 491       stopwatchDet.Start();
 492       loader->LoadDigits("update");
 493       loader->CleanDigits();
 494       loader->MakeDigitsContainer();
 495       TTree* digitsTree = loader->TreeD();
 496       reconstructor->ConvertDigits(fRawReader, digitsTree);
 497       loader->WriteDigits("OVERWRITE");
 498       loader->UnloadDigits();
 499       AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
 500                    fgkDetectorName[iDet],
 501                    stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 502     }
 503
 504     // local reconstruction
 505     if (!reconstructor->HasLocalReconstruction()) continue;
 506     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
 507     TStopwatch stopwatchDet;
 508     stopwatchDet.Start();
 509     loader->LoadRecPoints("update");
 510     loader->CleanRecPoints();
 511     loader->MakeRecPointsContainer();
 512     TTree* clustersTree = loader->TreeR();
 513     if (fRawReader && !reconstructor->HasDigitConversion()) {
 514       reconstructor->Reconstruct(fRawReader, clustersTree);
 515     } else {
 516       loader->LoadDigits("read");
 517       TTree* digitsTree = loader->TreeD();
 518       if (!digitsTree) {
 519         AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
 520         if (fStopOnError) return kFALSE;
 521       } else {
 522         reconstructor->Reconstruct(digitsTree, clustersTree);
 523       }
 524       loader->UnloadDigits();
 525     }
 526     loader->WriteRecPoints("OVERWRITE");
 527     loader->UnloadRecPoints();
 528     AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
 529                     fgkDetectorName[iDet],
 530                     stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
 531   }
 532
 533   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 534     AliError(Form("the following detectors were not found: %s",
 535                   detStr.Data()));
 536     if (fStopOnError) return kFALSE;
 537   }
 538
 539   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 540                stopwatch.RealTime(),stopwatch.CpuTime()));
 541
 542   return kTRUE;
 543 }
 544
 545 //_____________________________________________________________________________
 546 Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
 547 {
 548 // run the barrel tracking
 549
 550   TStopwatch stopwatch;
 551   stopwatch.Start();
 552
 553   AliESDVertex* vertex = NULL;
 554   Double_t vtxPos[3] = {0, 0, 0};
 555   Double_t vtxErr[3] = {0.07, 0.07, 0.1};
 556   TArrayF mcVertex(3);
 557   if (fRunLoader->GetHeader() && fRunLoader->GetHeader()->GenEventHeader()) {
 558     fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 559     for (Int_t i = 0; i < 3; i++) vtxPos[i] = mcVertex[i];
 560   }
 561
 562   if (fVertexer) {
 563     AliInfo("running the ITS vertex finder");
 564     if (fLoader[0]) fLoader[0]->LoadRecPoints();
 565     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
 566     if (fLoader[0]) fLoader[0]->UnloadRecPoints();
 567     if(!vertex){
 568       AliWarning("Vertex not found");
 569       vertex = new AliESDVertex();
 570     }
 571     else {
 572       vertex->SetTruePos(vtxPos);  // store also the vertex from MC
 573     }
 574
 575   } else {
 576     AliInfo("getting the primary vertex from MC");
 577     vertex = new AliESDVertex(vtxPos, vtxErr);
 578   }
 579
 580   if (vertex) {
 581     vertex->GetXYZ(vtxPos);
 582     vertex->GetSigmaXYZ(vtxErr);
 583   } else {
 584     AliWarning("no vertex reconstructed");
 585     vertex = new AliESDVertex(vtxPos, vtxErr);
 586   }
 587   esd->SetVertex(vertex);
 588   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 589     if (fTracker[iDet]) fTracker[iDet]->SetVertex(vtxPos, vtxErr);
 590   }
 591   delete vertex;
 592
 593   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 594                stopwatch.RealTime(),stopwatch.CpuTime()));
 595
 596   return kTRUE;
 597 }
 598
 599 //_____________________________________________________________________________
 600 Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
 601 {
 602 // run the HLT barrel tracking
 603
 604   TStopwatch stopwatch;
 605   stopwatch.Start();
 606
 607   if (!fRunLoader) {
 608     AliError("Missing runLoader!");
 609     return kFALSE;
 610   }
 611
 612   AliInfo("running HLT tracking");
 613
 614   // Get a pointer to the HLT reconstructor
 615   AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
 616   if (!reconstructor) return kFALSE;
 617
 618   // TPC + ITS
 619   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 620     TString detName = fgkDetectorName[iDet];
 621     AliDebug(1, Form("%s HLT tracking", detName.Data()));
 622     reconstructor->SetOption(detName.Data());
 623     AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
 624     if (!tracker) {
 625       AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
 626       if (fStopOnError) return kFALSE;
 627       continue;
 628     }
 629     Double_t vtxPos[3];
 630     Double_t vtxErr[3]={0.005,0.005,0.010};
 631     const AliESDVertex *vertex = esd->GetVertex();
 632     vertex->GetXYZ(vtxPos);
 633     tracker->SetVertex(vtxPos,vtxErr);
 634     if(iDet != 1) {
 635       fLoader[iDet]->LoadRecPoints("read");
 636       TTree* tree = fLoader[iDet]->TreeR();
 637       if (!tree) {
 638         AliError(Form("Can't get the %s cluster tree", detName.Data()));
 639         return kFALSE;
 640       }
 641       tracker->LoadClusters(tree);
 642     }
 643     if (tracker->Clusters2Tracks(esd) != 0) {
 644       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
 645       return kFALSE;
 646     }
 647     if(iDet != 1) {
 648       tracker->UnloadClusters();
 649     }
 650     delete tracker;
 651   }
 652
 653   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 654                stopwatch.RealTime(),stopwatch.CpuTime()));
 655
 656   return kTRUE;
 657 }
 658
 659 //_____________________________________________________________________________
 660 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
 661 {
 662 // run the barrel tracking
 663
 664   TStopwatch stopwatch;
 665   stopwatch.Start();
 666
 667   AliInfo("running tracking");
 668
 669   // pass 1: TPC + ITS inwards
 670   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 671     if (!fTracker[iDet]) continue;
 672     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
 673
 674     // load clusters
 675     fLoader[iDet]->LoadRecPoints("read");
 676     TTree* tree = fLoader[iDet]->TreeR();
 677     if (!tree) {
 678       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
 679       return kFALSE;
 680     }
 681     fTracker[iDet]->LoadClusters(tree);
 682
 683     // run tracking
 684     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
 685       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
 686       return kFALSE;
 687     }
 688     if (fCheckPointLevel > 1) {
 689       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
 690     }
 691     // preliminary PID in TPC needed by the ITS tracker
 692     if (iDet == 1) {
 693       GetReconstructor(1)->FillESD(fRunLoader, esd);
 694       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
 695       AliESDpid::MakePID(esd);
 696     }
 697   }
 698
 699   // pass 2: ALL backwards
 700   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 701     if (!fTracker[iDet]) continue;
 702     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
 703
 704     // load clusters
 705     if (iDet > 1) {     // all except ITS, TPC
 706       TTree* tree = NULL;
 707       if (iDet == 3) {   // TOF
 708         fLoader[iDet]->LoadDigits("read");
 709         tree = fLoader[iDet]->TreeD();
 710       } else {
 711         fLoader[iDet]->LoadRecPoints("read");
 712         tree = fLoader[iDet]->TreeR();
 713       }
 714       if (!tree) {
 715         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
 716         return kFALSE;
 717       }
 718       fTracker[iDet]->LoadClusters(tree);
 719     }
 720
 721     // run tracking
 722     if (fTracker[iDet]->PropagateBack(esd) != 0) {
 723       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
 724       return kFALSE;
 725     }
 726     if (fCheckPointLevel > 1) {
 727       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
 728     }
 729
 730     // unload clusters
 731     if (iDet > 2) {     // all except ITS, TPC, TRD
 732       fTracker[iDet]->UnloadClusters();
 733       if (iDet == 3) {   // TOF
 734         fLoader[iDet]->UnloadDigits();
 735       } else {
 736         fLoader[iDet]->UnloadRecPoints();
 737       }
 738     }
 739   }
 740
 741   // pass 3: TRD + TPC + ITS refit inwards
 742   for (Int_t iDet = 2; iDet >= 0; iDet--) {
 743     if (!fTracker[iDet]) continue;
 744     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
 745
 746     // run tracking
 747     if (fTracker[iDet]->RefitInward(esd) != 0) {
 748       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
 749       return kFALSE;
 750     }
 751     if (fCheckPointLevel > 1) {
 752       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
 753     }
 754
 755     // unload clusters
 756     fTracker[iDet]->UnloadClusters();
 757     fLoader[iDet]->UnloadRecPoints();
 758   }
 759
 760   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 761                stopwatch.RealTime(),stopwatch.CpuTime()));
 762
 763   return kTRUE;
 764 }
 765
 766 //_____________________________________________________________________________
 767 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
 768 {
 769 // fill the event summary data
 770
 771   TStopwatch stopwatch;
 772   stopwatch.Start();
 773   AliInfo("filling ESD");
 774
 775   TString detStr = detectors;
 776   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 777     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 778     AliReconstructor* reconstructor = GetReconstructor(iDet);
 779     if (!reconstructor) continue;
 780
 781     if (!ReadESD(esd, fgkDetectorName[iDet])) {
 782       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
 783       TTree* clustersTree = NULL;
 784       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
 785         fLoader[iDet]->LoadRecPoints("read");
 786         clustersTree = fLoader[iDet]->TreeR();
 787         if (!clustersTree) {
 788           AliError(Form("Can't get the %s clusters tree",
 789                         fgkDetectorName[iDet]));
 790           if (fStopOnError) return kFALSE;
 791         }
 792       }
 793       if (fRawReader && !reconstructor->HasDigitConversion()) {
 794         reconstructor->FillESD(fRawReader, clustersTree, esd);
 795       } else {
 796         TTree* digitsTree = NULL;
 797         if (fLoader[iDet]) {
 798           fLoader[iDet]->LoadDigits("read");
 799           digitsTree = fLoader[iDet]->TreeD();
 800           if (!digitsTree) {
 801             AliError(Form("Can't get the %s digits tree",
 802                           fgkDetectorName[iDet]));
 803             if (fStopOnError) return kFALSE;
 804           }
 805         }
 806         reconstructor->FillESD(digitsTree, clustersTree, esd);
 807         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
 808       }
 809       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
 810         fLoader[iDet]->UnloadRecPoints();
 811       }
 812
 813       if (fRawReader) {
 814         reconstructor->FillESD(fRunLoader, fRawReader, esd);
 815       } else {
 816         reconstructor->FillESD(fRunLoader, esd);
 817       }
 818       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
 819     }
 820   }
 821
 822   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 823     AliError(Form("the following detectors were not found: %s",
 824                   detStr.Data()));
 825     if (fStopOnError) return kFALSE;
 826   }
 827
 828   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 829                stopwatch.RealTime(),stopwatch.CpuTime()));
 830
 831   return kTRUE;
 832 }
 833
 834
 835 //_____________________________________________________________________________
 836 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
 837 {
 838 // check whether detName is contained in detectors
 839 // if yes, it is removed from detectors
 840
 841   // check if all detectors are selected
 842   if ((detectors.CompareTo("ALL") == 0) ||
 843       detectors.BeginsWith("ALL ") ||
 844       detectors.EndsWith(" ALL") ||
 845       detectors.Contains(" ALL ")) {
 846     detectors = "ALL";
 847     return kTRUE;
 848   }
 849
 850   // search for the given detector
 851   Bool_t result = kFALSE;
 852   if ((detectors.CompareTo(detName) == 0) ||
 853       detectors.BeginsWith(detName+" ") ||
 854       detectors.EndsWith(" "+detName) ||
 855       detectors.Contains(" "+detName+" ")) {
 856     detectors.ReplaceAll(detName, "");
 857     result = kTRUE;
 858   }
 859
 860   // clean up the detectors string
 861   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
 862   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
 863   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
 864
 865   return result;
 866 }
 867
 868 //_____________________________________________________________________________
 869 Bool_t AliReconstruction::InitRunLoader()
 870 {
 871 // get or create the run loader
 872
 873   if (gAlice) delete gAlice;
 874   gAlice = NULL;
 875
 876   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
 877     // load all base libraries to get the loader classes
 878     TString libs = gSystem->GetLibraries();
 879     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 880       TString detName = fgkDetectorName[iDet];
 881       if (detName == "HLT") continue;
 882       if (libs.Contains("lib" + detName + "base.so")) continue;
 883       gSystem->Load("lib" + detName + "base.so");
 884     }
 885     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
 886     if (!fRunLoader) {
 887       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
 888       CleanUp();
 889       return kFALSE;
 890     }
 891     fRunLoader->CdGAFile();
 892     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
 893       if (fRunLoader->LoadgAlice() == 0) {
 894         gAlice = fRunLoader->GetAliRun();
 895         AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
 896       }
 897     }
 898     if (!gAlice && !fRawReader) {
 899       AliError(Form("no gAlice object found in file %s",
 900                     fGAliceFileName.Data()));
 901       CleanUp();
 902       return kFALSE;
 903     }
 904
 905   } else {               // galice.root does not exist
 906     if (!fRawReader) {
 907       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
 908       CleanUp();
 909       return kFALSE;
 910     }
 911     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
 912                                     AliConfig::GetDefaultEventFolderName(),
 913                                     "recreate");
 914     if (!fRunLoader) {
 915       AliError(Form("could not create run loader in file %s",
 916                     fGAliceFileName.Data()));
 917       CleanUp();
 918       return kFALSE;
 919     }
 920     fRunLoader->MakeTree("E");
 921     Int_t iEvent = 0;
 922     while (fRawReader->NextEvent()) {
 923       fRunLoader->SetEventNumber(iEvent);
 924       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
 925                                      iEvent, iEvent);
 926       fRunLoader->MakeTree("H");
 927       fRunLoader->TreeE()->Fill();
 928       iEvent++;
 929     }
 930     fRawReader->RewindEvents();
 931     fRunLoader->WriteHeader("OVERWRITE");
 932     fRunLoader->CdGAFile();
 933     fRunLoader->Write(0, TObject::kOverwrite);
 934 //    AliTracker::SetFieldMap(???);
 935   }
 936
 937   return kTRUE;
 938 }
 939
 940 //_____________________________________________________________________________
 941 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
 942 {
 943 // get the reconstructor object and the loader for a detector
 944
 945   if (fReconstructor[iDet]) return fReconstructor[iDet];
 946
 947   // load the reconstructor object
 948   TPluginManager* pluginManager = gROOT->GetPluginManager();
 949   TString detName = fgkDetectorName[iDet];
 950   TString recName = "Ali" + detName + "Reconstructor";
 951   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
 952
 953   if (detName == "HLT") {
 954     if (!gROOT->GetClass("AliLevel3")) {
 955       gSystem->Load("libAliL3Src.so");
 956       gSystem->Load("libAliL3Misc.so");
 957       gSystem->Load("libAliL3Hough.so");
 958       gSystem->Load("libAliL3Comp.so");
 959     }
 960   }
 961
 962   AliReconstructor* reconstructor = NULL;
 963   // first check if a plugin is defined for the reconstructor
 964   TPluginHandler* pluginHandler =
 965     pluginManager->FindHandler("AliReconstructor", detName);
 966   // if not, add a plugin for it
 967   if (!pluginHandler) {
 968     AliDebug(1, Form("defining plugin for %s", recName.Data()));
 969     TString libs = gSystem->GetLibraries();
 970     if (libs.Contains("lib" + detName + "base.so") ||
 971         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
 972       pluginManager->AddHandler("AliReconstructor", detName,
 973                                 recName, detName + "rec", recName + "()");
 974     } else {
 975       pluginManager->AddHandler("AliReconstructor", detName,
 976                                 recName, detName, recName + "()");
 977     }
 978     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
 979   }
 980   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
 981     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
 982   }
 983   if (reconstructor) {
 984     TObject* obj = fOptions.FindObject(detName.Data());
 985     if (obj) reconstructor->SetOption(obj->GetTitle());
 986     reconstructor->Init(fRunLoader);
 987     fReconstructor[iDet] = reconstructor;
 988   }
 989
 990   // get or create the loader
 991   if (detName != "HLT") {
 992     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
 993     if (!fLoader[iDet]) {
 994       AliConfig::Instance()
 995         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
 996                                 detName, detName);
 997       // first check if a plugin is defined for the loader
 998       TPluginHandler* pluginHandler =
 999         pluginManager->FindHandler("AliLoader", detName);
1000       // if not, add a plugin for it
1001       if (!pluginHandler) {
1002         TString loaderName = "Ali" + detName + "Loader";
1003         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
1004         pluginManager->AddHandler("AliLoader", detName,
1005                                   loaderName, detName + "base",
1006                                   loaderName + "(const char*, TFolder*)");
1007         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1008       }
1009       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1010         fLoader[iDet] =
1011           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1012                                                  fRunLoader->GetEventFolder());
1013       }
1014       if (!fLoader[iDet]) {   // use default loader
1015         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1016       }
1017       if (!fLoader[iDet]) {
1018         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1019         if (fStopOnError) return NULL;
1020       } else {
1021         fRunLoader->AddLoader(fLoader[iDet]);
1022         fRunLoader->CdGAFile();
1023         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1024         fRunLoader->Write(0, TObject::kOverwrite);
1025       }
1026     }
1027   }
1028
1029   return reconstructor;
1030 }
1031
1032 //_____________________________________________________________________________
1033 Bool_t AliReconstruction::CreateVertexer()
1034 {
1035 // create the vertexer
1036
1037   fVertexer = NULL;
1038   AliReconstructor* itsReconstructor = GetReconstructor(0);
1039   if (itsReconstructor) {
1040     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1041   }
1042   if (!fVertexer) {
1043     AliWarning("couldn't create a vertexer for ITS");
1044     if (fStopOnError) return kFALSE;
1045   }
1046
1047   return kTRUE;
1048 }
1049
1050 //_____________________________________________________________________________
1051 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1052 {
1053 // create the trackers
1054
1055   TString detStr = detectors;
1056   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1057     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1058     AliReconstructor* reconstructor = GetReconstructor(iDet);
1059     if (!reconstructor) continue;
1060     TString detName = fgkDetectorName[iDet];
1061     if (detName == "HLT") {
1062       fRunHLTTracking = kTRUE;
1063       continue;
1064     }
1065
1066     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1067     if (!fTracker[iDet] && (iDet < 7)) {
1068       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1069       if (fStopOnError) return kFALSE;
1070     }
1071   }
1072
1073   return kTRUE;
1074 }
1075
1076 //_____________________________________________________________________________
1077 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1078 {
1079 // delete trackers and the run loader and close and delete the file
1080
1081   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1082     delete fReconstructor[iDet];
1083     fReconstructor[iDet] = NULL;
1084     fLoader[iDet] = NULL;
1085     delete fTracker[iDet];
1086     fTracker[iDet] = NULL;
1087   }
1088   delete fVertexer;
1089   fVertexer = NULL;
1090
1091   delete fRunLoader;
1092   fRunLoader = NULL;
1093   delete fRawReader;
1094   fRawReader = NULL;
1095
1096   if (file) {
1097     file->Close();
1098     delete file;
1099   }
1100
1101   if (fileOld) {
1102     fileOld->Close();
1103     delete fileOld;
1104     gSystem->Unlink("AliESDs.old.root");
1105   }
1106 }
1107
1108
1109 //_____________________________________________________________________________
1110 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1111 {
1112 // read the ESD event from a file
1113
1114   if (!esd) return kFALSE;
1115   char fileName[256];
1116   sprintf(fileName, "ESD_%d.%d_%s.root",
1117           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1118   if (gSystem->AccessPathName(fileName)) return kFALSE;
1119
1120   AliDebug(1, Form("reading ESD from file %s", fileName));
1121   TFile* file = TFile::Open(fileName);
1122   if (!file || !file->IsOpen()) {
1123     AliError(Form("opening %s failed", fileName));
1124     delete file;
1125     return kFALSE;
1126   }
1127
1128   gROOT->cd();
1129   delete esd;
1130   esd = (AliESD*) file->Get("ESD");
1131   file->Close();
1132   delete file;
1133   return kTRUE;
1134 }
1135
1136 //_____________________________________________________________________________
1137 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1138 {
1139 // write the ESD event to a file
1140
1141   if (!esd) return;
1142   char fileName[256];
1143   sprintf(fileName, "ESD_%d.%d_%s.root",
1144           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1145
1146   AliDebug(1, Form("writing ESD to file %s", fileName));
1147   TFile* file = TFile::Open(fileName, "recreate");
1148   if (!file || !file->IsOpen()) {
1149     AliError(Form("opening %s failed", fileName));
1150   } else {
1151     esd->Write("ESD");
1152     file->Close();
1153   }
1154   delete file;
1155 }