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       continue;
 622     }
 623     Double_t vtxPos[3];
 624     Double_t vtxErr[3]={0.005,0.005,0.010};
 625     const AliESDVertex *vertex = esd->GetVertex();
 626     vertex->GetXYZ(vtxPos);
 627     tracker->SetVertex(vtxPos,vtxErr);
 628     if(iDet != 1) {
 629       fLoader[iDet]->LoadRecPoints("read");
 630       TTree* tree = fLoader[iDet]->TreeR();
 631       if (!tree) {
 632         AliError(Form("Can't get the %s cluster tree", detName.Data()));
 633         return kFALSE;
 634       }
 635       tracker->LoadClusters(tree);
 636     }
 637     if (tracker->Clusters2Tracks(esd) != 0) {
 638       AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
 639       return kFALSE;
 640     }
 641     if(iDet != 1) {
 642       tracker->UnloadClusters();
 643     }
 644     delete tracker;
 645   }
 646
 647   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 648                stopwatch.RealTime(),stopwatch.CpuTime()));
 649
 650   return kTRUE;
 651 }
 652
 653 //_____________________________________________________________________________
 654 Bool_t AliReconstruction::RunTracking(AliESD*& esd)
 655 {
 656 // run the barrel tracking
 657
 658   TStopwatch stopwatch;
 659   stopwatch.Start();
 660
 661   AliInfo("running tracking");
 662
 663   // pass 1: TPC + ITS inwards
 664   for (Int_t iDet = 1; iDet >= 0; iDet--) {
 665     if (!fTracker[iDet]) continue;
 666     AliDebug(1, Form("%s tracking", fgkDetectorName[iDet]));
 667
 668     // load clusters
 669     fLoader[iDet]->LoadRecPoints("read");
 670     TTree* tree = fLoader[iDet]->TreeR();
 671     if (!tree) {
 672       AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
 673       return kFALSE;
 674     }
 675     fTracker[iDet]->LoadClusters(tree);
 676
 677     // run tracking
 678     if (fTracker[iDet]->Clusters2Tracks(esd) != 0) {
 679       AliError(Form("%s Clusters2Tracks failed", fgkDetectorName[iDet]));
 680       return kFALSE;
 681     }
 682     if (fCheckPointLevel > 1) {
 683       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
 684     }
 685     // preliminary PID in TPC needed by the ITS tracker
 686     if (iDet == 1) {
 687       GetReconstructor(1)->FillESD(fRunLoader, esd);
 688       GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
 689       AliESDpid::MakePID(esd);
 690     }
 691   }
 692
 693   // pass 2: ALL backwards
 694   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 695     if (!fTracker[iDet]) continue;
 696     AliDebug(1, Form("%s back propagation", fgkDetectorName[iDet]));
 697
 698     // load clusters
 699     if (iDet > 1) {     // all except ITS, TPC
 700       TTree* tree = NULL;
 701       if (iDet == 3) {   // TOF
 702         fLoader[iDet]->LoadDigits("read");
 703         tree = fLoader[iDet]->TreeD();
 704       } else {
 705         fLoader[iDet]->LoadRecPoints("read");
 706         tree = fLoader[iDet]->TreeR();
 707       }
 708       if (!tree) {
 709         AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
 710         return kFALSE;
 711       }
 712       fTracker[iDet]->LoadClusters(tree);
 713     }
 714
 715     // run tracking
 716     if (fTracker[iDet]->PropagateBack(esd) != 0) {
 717       AliError(Form("%s backward propagation failed", fgkDetectorName[iDet]));
 718       return kFALSE;
 719     }
 720     if (fCheckPointLevel > 1) {
 721       WriteESD(esd, Form("%s.back", fgkDetectorName[iDet]));
 722     }
 723
 724     // unload clusters
 725     if (iDet > 2) {     // all except ITS, TPC, TRD
 726       fTracker[iDet]->UnloadClusters();
 727       if (iDet == 3) {   // TOF
 728         fLoader[iDet]->UnloadDigits();
 729       } else {
 730         fLoader[iDet]->UnloadRecPoints();
 731       }
 732     }
 733   }
 734
 735   // pass 3: TRD + TPC + ITS refit inwards
 736   for (Int_t iDet = 2; iDet >= 0; iDet--) {
 737     if (!fTracker[iDet]) continue;
 738     AliDebug(1, Form("%s inward refit", fgkDetectorName[iDet]));
 739
 740     // run tracking
 741     if (fTracker[iDet]->RefitInward(esd) != 0) {
 742       AliError(Form("%s inward refit failed", fgkDetectorName[iDet]));
 743       return kFALSE;
 744     }
 745     if (fCheckPointLevel > 1) {
 746       WriteESD(esd, Form("%s.refit", fgkDetectorName[iDet]));
 747     }
 748
 749     // unload clusters
 750     fTracker[iDet]->UnloadClusters();
 751     fLoader[iDet]->UnloadRecPoints();
 752   }
 753
 754   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 755                stopwatch.RealTime(),stopwatch.CpuTime()));
 756
 757   return kTRUE;
 758 }
 759
 760 //_____________________________________________________________________________
 761 Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
 762 {
 763 // fill the event summary data
 764
 765   TStopwatch stopwatch;
 766   stopwatch.Start();
 767   AliInfo("filling ESD");
 768
 769   TString detStr = detectors;
 770   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 771     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
 772     AliReconstructor* reconstructor = GetReconstructor(iDet);
 773     if (!reconstructor) continue;
 774
 775     if (!ReadESD(esd, fgkDetectorName[iDet])) {
 776       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
 777       TTree* clustersTree = NULL;
 778       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
 779         fLoader[iDet]->LoadRecPoints("read");
 780         clustersTree = fLoader[iDet]->TreeR();
 781         if (!clustersTree) {
 782           AliError(Form("Can't get the %s clusters tree",
 783                         fgkDetectorName[iDet]));
 784           if (fStopOnError) return kFALSE;
 785         }
 786       }
 787       if (fRawReader && !reconstructor->HasDigitConversion()) {
 788         reconstructor->FillESD(fRawReader, clustersTree, esd);
 789       } else {
 790         TTree* digitsTree = NULL;
 791         if (fLoader[iDet]) {
 792           fLoader[iDet]->LoadDigits("read");
 793           digitsTree = fLoader[iDet]->TreeD();
 794           if (!digitsTree) {
 795             AliError(Form("Can't get the %s digits tree",
 796                           fgkDetectorName[iDet]));
 797             if (fStopOnError) return kFALSE;
 798           }
 799         }
 800         reconstructor->FillESD(digitsTree, clustersTree, esd);
 801         if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
 802       }
 803       if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
 804         fLoader[iDet]->UnloadRecPoints();
 805       }
 806
 807       if (fRawReader) {
 808         reconstructor->FillESD(fRunLoader, fRawReader, esd);
 809       } else {
 810         reconstructor->FillESD(fRunLoader, esd);
 811       }
 812       if (fCheckPointLevel > 2) WriteESD(esd, fgkDetectorName[iDet]);
 813     }
 814   }
 815
 816   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
 817     AliError(Form("the following detectors were not found: %s",
 818                   detStr.Data()));
 819     if (fStopOnError) return kFALSE;
 820   }
 821
 822   AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
 823                stopwatch.RealTime(),stopwatch.CpuTime()));
 824
 825   return kTRUE;
 826 }
 827
 828
 829 //_____________________________________________________________________________
 830 Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
 831 {
 832 // check whether detName is contained in detectors
 833 // if yes, it is removed from detectors
 834
 835   // check if all detectors are selected
 836   if ((detectors.CompareTo("ALL") == 0) ||
 837       detectors.BeginsWith("ALL ") ||
 838       detectors.EndsWith(" ALL") ||
 839       detectors.Contains(" ALL ")) {
 840     detectors = "ALL";
 841     return kTRUE;
 842   }
 843
 844   // search for the given detector
 845   Bool_t result = kFALSE;
 846   if ((detectors.CompareTo(detName) == 0) ||
 847       detectors.BeginsWith(detName+" ") ||
 848       detectors.EndsWith(" "+detName) ||
 849       detectors.Contains(" "+detName+" ")) {
 850     detectors.ReplaceAll(detName, "");
 851     result = kTRUE;
 852   }
 853
 854   // clean up the detectors string
 855   while (detectors.Contains("  ")) detectors.ReplaceAll("  ", " ");
 856   while (detectors.BeginsWith(" ")) detectors.Remove(0, 1);
 857   while (detectors.EndsWith(" ")) detectors.Remove(detectors.Length()-1, 1);
 858
 859   return result;
 860 }
 861
 862 //_____________________________________________________________________________
 863 Bool_t AliReconstruction::InitRunLoader()
 864 {
 865 // get or create the run loader
 866
 867   if (gAlice) delete gAlice;
 868   gAlice = NULL;
 869
 870   if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
 871     // load all base libraries to get the loader classes
 872     TString libs = gSystem->GetLibraries();
 873     for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
 874       TString detName = fgkDetectorName[iDet];
 875       if (detName == "HLT") continue;
 876       if (libs.Contains("lib" + detName + "base.so")) continue;
 877       gSystem->Load("lib" + detName + "base.so");
 878     }
 879     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
 880     if (!fRunLoader) {
 881       AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
 882       CleanUp();
 883       return kFALSE;
 884     }
 885     fRunLoader->CdGAFile();
 886     if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
 887       if (fRunLoader->LoadgAlice() == 0) {
 888         gAlice = fRunLoader->GetAliRun();
 889         AliTracker::SetFieldMap(gAlice->Field());
 890       }
 891     }
 892     if (!gAlice && !fRawReader) {
 893       AliError(Form("no gAlice object found in file %s",
 894                     fGAliceFileName.Data()));
 895       CleanUp();
 896       return kFALSE;
 897     }
 898
 899   } else {               // galice.root does not exist
 900     if (!fRawReader) {
 901       AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
 902       CleanUp();
 903       return kFALSE;
 904     }
 905     fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
 906                                     AliConfig::GetDefaultEventFolderName(),
 907                                     "recreate");
 908     if (!fRunLoader) {
 909       AliError(Form("could not create run loader in file %s",
 910                     fGAliceFileName.Data()));
 911       CleanUp();
 912       return kFALSE;
 913     }
 914     fRunLoader->MakeTree("E");
 915     Int_t iEvent = 0;
 916     while (fRawReader->NextEvent()) {
 917       fRunLoader->SetEventNumber(iEvent);
 918       fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(),
 919                                      iEvent, iEvent);
 920       fRunLoader->MakeTree("H");
 921       fRunLoader->TreeE()->Fill();
 922       iEvent++;
 923     }
 924     fRawReader->RewindEvents();
 925     fRunLoader->WriteHeader("OVERWRITE");
 926     fRunLoader->CdGAFile();
 927     fRunLoader->Write(0, TObject::kOverwrite);
 928 //    AliTracker::SetFieldMap(???);
 929   }
 930
 931   return kTRUE;
 932 }
 933
 934 //_____________________________________________________________________________
 935 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
 936 {
 937 // get the reconstructor object and the loader for a detector
 938
 939   if (fReconstructor[iDet]) return fReconstructor[iDet];
 940
 941   // load the reconstructor object
 942   TPluginManager* pluginManager = gROOT->GetPluginManager();
 943   TString detName = fgkDetectorName[iDet];
 944   TString recName = "Ali" + detName + "Reconstructor";
 945   if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
 946
 947   if (detName == "HLT") {
 948     if (!gROOT->GetClass("AliLevel3")) {
 949       gSystem->Load("libAliL3Src.so");
 950       gSystem->Load("libAliL3Misc.so");
 951       gSystem->Load("libAliL3Hough.so");
 952       gSystem->Load("libAliL3Comp.so");
 953     }
 954   }
 955
 956   AliReconstructor* reconstructor = NULL;
 957   // first check if a plugin is defined for the reconstructor
 958   TPluginHandler* pluginHandler =
 959     pluginManager->FindHandler("AliReconstructor", detName);
 960   // if not, add a plugin for it
 961   if (!pluginHandler) {
 962     AliDebug(1, Form("defining plugin for %s", recName.Data()));
 963     TString libs = gSystem->GetLibraries();
 964     if (libs.Contains("lib" + detName + "base.so") ||
 965         (gSystem->Load("lib" + detName + "base.so") >= 0)) {
 966       pluginManager->AddHandler("AliReconstructor", detName,
 967                                 recName, detName + "rec", recName + "()");
 968     } else {
 969       pluginManager->AddHandler("AliReconstructor", detName,
 970                                 recName, detName, recName + "()");
 971     }
 972     pluginHandler = pluginManager->FindHandler("AliReconstructor", detName);
 973   }
 974   if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
 975     reconstructor = (AliReconstructor*) pluginHandler->ExecPlugin(0);
 976   }
 977   if (reconstructor) {
 978     TObject* obj = fOptions.FindObject(detName.Data());
 979     if (obj) reconstructor->SetOption(obj->GetTitle());
 980     reconstructor->Init(fRunLoader);
 981     fReconstructor[iDet] = reconstructor;
 982   }
 983
 984   // get or create the loader
 985   if (detName != "HLT") {
 986     fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
 987     if (!fLoader[iDet]) {
 988       AliConfig::Instance()
 989         ->CreateDetectorFolders(fRunLoader->GetEventFolder(),
 990                                 detName, detName);
 991       // first check if a plugin is defined for the loader
 992       TPluginHandler* pluginHandler =
 993         pluginManager->FindHandler("AliLoader", detName);
 994       // if not, add a plugin for it
 995       if (!pluginHandler) {
 996         TString loaderName = "Ali" + detName + "Loader";
 997         AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
 998         pluginManager->AddHandler("AliLoader", detName,
 999                                   loaderName, detName + "base",
1000                                   loaderName + "(const char*, TFolder*)");
1001         pluginHandler = pluginManager->FindHandler("AliLoader", detName);
1002       }
1003       if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
1004         fLoader[iDet] =
1005           (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(),
1006                                                  fRunLoader->GetEventFolder());
1007       }
1008       if (!fLoader[iDet]) {   // use default loader
1009         fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
1010       }
1011       if (!fLoader[iDet]) {
1012         AliWarning(Form("couldn't get loader for %s", detName.Data()));
1013         if (fStopOnError) return NULL;
1014       } else {
1015         fRunLoader->AddLoader(fLoader[iDet]);
1016         fRunLoader->CdGAFile();
1017         if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
1018         fRunLoader->Write(0, TObject::kOverwrite);
1019       }
1020     }
1021   }
1022
1023   return reconstructor;
1024 }
1025
1026 //_____________________________________________________________________________
1027 Bool_t AliReconstruction::CreateVertexer()
1028 {
1029 // create the vertexer
1030
1031   fVertexer = NULL;
1032   AliReconstructor* itsReconstructor = GetReconstructor(0);
1033   if (itsReconstructor) {
1034     fVertexer = itsReconstructor->CreateVertexer(fRunLoader);
1035   }
1036   if (!fVertexer) {
1037     AliWarning("couldn't create a vertexer for ITS");
1038     if (fStopOnError) return kFALSE;
1039   }
1040
1041   return kTRUE;
1042 }
1043
1044 //_____________________________________________________________________________
1045 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
1046 {
1047 // create the trackers
1048
1049   TString detStr = detectors;
1050   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1051     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
1052     AliReconstructor* reconstructor = GetReconstructor(iDet);
1053     if (!reconstructor) continue;
1054     TString detName = fgkDetectorName[iDet];
1055     if (detName == "HLT") {
1056       fRunHLTTracking = kTRUE;
1057       continue;
1058     }
1059
1060     fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
1061     if (!fTracker[iDet] && (iDet < 7)) {
1062       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
1063       if (fStopOnError) return kFALSE;
1064     }
1065   }
1066
1067   return kTRUE;
1068 }
1069
1070 //_____________________________________________________________________________
1071 void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
1072 {
1073 // delete trackers and the run loader and close and delete the file
1074
1075   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
1076     delete fReconstructor[iDet];
1077     fReconstructor[iDet] = NULL;
1078     fLoader[iDet] = NULL;
1079     delete fTracker[iDet];
1080     fTracker[iDet] = NULL;
1081   }
1082   delete fVertexer;
1083   fVertexer = NULL;
1084
1085   delete fRunLoader;
1086   fRunLoader = NULL;
1087   delete fRawReader;
1088   fRawReader = NULL;
1089
1090   if (file) {
1091     file->Close();
1092     delete file;
1093   }
1094
1095   if (fileOld) {
1096     fileOld->Close();
1097     delete fileOld;
1098     gSystem->Unlink("AliESDs.old.root");
1099   }
1100 }
1101
1102
1103 //_____________________________________________________________________________
1104 Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
1105 {
1106 // read the ESD event from a file
1107
1108   if (!esd) return kFALSE;
1109   char fileName[256];
1110   sprintf(fileName, "ESD_%d.%d_%s.root",
1111           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1112   if (gSystem->AccessPathName(fileName)) return kFALSE;
1113
1114   AliDebug(1, Form("reading ESD from file %s", fileName));
1115   TFile* file = TFile::Open(fileName);
1116   if (!file || !file->IsOpen()) {
1117     AliError(Form("opening %s failed", fileName));
1118     delete file;
1119     return kFALSE;
1120   }
1121
1122   gROOT->cd();
1123   delete esd;
1124   esd = (AliESD*) file->Get("ESD");
1125   file->Close();
1126   delete file;
1127   return kTRUE;
1128 }
1129
1130 //_____________________________________________________________________________
1131 void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
1132 {
1133 // write the ESD event to a file
1134
1135   if (!esd) return;
1136   char fileName[256];
1137   sprintf(fileName, "ESD_%d.%d_%s.root",
1138           esd->GetRunNumber(), esd->GetEventNumber(), recStep);
1139
1140   AliDebug(1, Form("writing ESD to file %s", fileName));
1141   TFile* file = TFile::Open(fileName, "recreate");
1142   if (!file || !file->IsOpen()) {
1143     AliError(Form("opening %s failed", fileName));
1144   } else {
1145     esd->Write("ESD");
1146     file->Close();
1147   }
1148   delete file;
1149 }