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