]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliReconstruction.cxx
Correct overloading of virtual functions in the derived classes (icc)
[u/mrichter/AliRoot.git] / STEER / AliReconstruction.cxx
index da311d5e8039ffe230986201b1f2cbfcdc7c393a..ff8731f5f873f95529b3ae15f6f97354278be48b 100644 (file)
 // - raw data DATE file: DATE file name, any other non-empty string          //
 // - MC root files     : empty string, default                               //
 //                                                                           //
+// By default all events are reconstructed. The reconstruction can be        //
+// limited to a range of events by giving the index of the first and the     //
+// last event as an argument to the Run method or by calling                 //
+//                                                                           //
+//   rec.SetEventRange(..., ...);                                            //
+//                                                                           //
+// The index -1 (default) can be used for the last event to indicate no      //
+// upper limit of the event range.                                           //
+//                                                                           //
 // The name of the galice file can be changed from the default               //
 // "galice.root" by passing it as argument to the AliReconstruction          //
 // constructor or by                                                         //
 //                                                                           //
 //   rec.SetRunTracking("...");                                              //
 //                                                                           //
+// Uniform/nonuniform field tracking switches (default: uniform field)       //
+//                                                                           //
+//   rec.SetUniformFieldTracking();  ( rec.SetNonuniformFieldTracking(); )   //
+//                                                                           //
 // The filling of additional ESD information can be steered by               //
 //                                                                           //
 //   rec.SetFillESD("...");                                                  //
 #include "AliRawReaderFile.h"
 #include "AliRawReaderDate.h"
 #include "AliRawReaderRoot.h"
-#include "AliTracker.h"
 #include "AliESD.h"
 #include "AliESDVertex.h"
+#include "AliTracker.h"
 #include "AliVertexer.h"
 #include "AliHeader.h"
 #include "AliGenEventHeader.h"
+#include "AliPID.h"
 #include "AliESDpid.h"
-#include "AliMagF.h"
+//#include "AliMagF.h"
+
+
+
+#include "AliRunTag.h"
+//#include "AliLHCTag.h"
+#include "AliDetectorTag.h"
+#include "AliEventTag.h"
+
+
 
 ClassImp(AliReconstruction)
 
@@ -128,11 +151,15 @@ AliReconstruction::AliReconstruction(const char* gAliceFilename,
   TNamed(name, title),
 
   fRunLocalReconstruction("ALL"),
+  fUniformField(kTRUE),
   fRunVertexFinder(kTRUE),
+  fRunHLTTracking(kFALSE),
   fRunTracking("ALL"),
   fFillESD("ALL"),
   fGAliceFileName(gAliceFilename),
   fInput(""),
+  fFirstEvent(0),
+  fLastEvent(-1),
   fStopOnError(kFALSE),
   fCheckPointLevel(0),
   fOptions(),
@@ -149,6 +176,7 @@ AliReconstruction::AliReconstruction(const char* gAliceFilename,
     fLoader[iDet] = NULL;
     fTracker[iDet] = NULL;
   }
+  AliPID pid;
 }
 
 //_____________________________________________________________________________
@@ -156,11 +184,15 @@ AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
   TNamed(rec),
 
   fRunLocalReconstruction(rec.fRunLocalReconstruction),
+  fUniformField(rec.fUniformField),
   fRunVertexFinder(rec.fRunVertexFinder),
+  fRunHLTTracking(rec.fRunHLTTracking),
   fRunTracking(rec.fRunTracking),
   fFillESD(rec.fFillESD),
   fGAliceFileName(rec.fGAliceFileName),
   fInput(rec.fInput),
+  fFirstEvent(rec.fFirstEvent),
+  fLastEvent(rec.fLastEvent),
   fStopOnError(rec.fStopOnError),
   fCheckPointLevel(0),
   fOptions(),
@@ -222,7 +254,8 @@ void AliReconstruction::SetOption(const char* detector, const char* option)
 
 
 //_____________________________________________________________________________
-Bool_t AliReconstruction::Run(const char* input)
+Bool_t AliReconstruction::Run(const char* input,
+                             Int_t firstEvent, Int_t lastEvent)
 {
 // run the reconstruction
 
@@ -238,23 +271,8 @@ Bool_t AliReconstruction::Run(const char* input)
     fRawReader->SelectEvents(7);
   }
 
-  // open the run loader
-  fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
-  if (!fRunLoader) {
-    AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
-    CleanUp();
-    return kFALSE;
-  }
-  fRunLoader->LoadgAlice();
-  AliRun* aliRun = fRunLoader->GetAliRun();
-  if (!aliRun) {
-    AliError(Form("no gAlice object found in file %s",
-                  fGAliceFileName.Data()));
-    CleanUp();
-    return kFALSE;
-  }
-  gAlice = aliRun;
-  AliTracker::SetFieldMap(gAlice->Field());
+  // get the run loader
+  if (!InitRunLoader()) return kFALSE;
 
   // local reconstruction
   if (!fRunLocalReconstruction.IsNull()) {
@@ -262,8 +280,8 @@ Bool_t AliReconstruction::Run(const char* input)
       if (fStopOnError) {CleanUp(); return kFALSE;}
     }
   }
-  if (!fRunVertexFinder && fRunTracking.IsNull() && 
-      fFillESD.IsNull()) return kTRUE;
+//  if (!fRunVertexFinder && fRunTracking.IsNull() && 
+//      fFillESD.IsNull()) return kTRUE;
 
   // get vertexer
   if (fRunVertexFinder && !CreateVertexer()) {
@@ -273,7 +291,7 @@ Bool_t AliReconstruction::Run(const char* input)
     }
   }
 
-  // get loaders and trackers
+  // get trackers
   if (!fRunTracking.IsNull() && !CreateTrackers(fRunTracking)) {
     if (fStopOnError) {
       CleanUp(); 
@@ -281,50 +299,108 @@ Bool_t AliReconstruction::Run(const char* input)
     }      
   }
 
+  // get the possibly already existing ESD file and tree
+  AliESD* esd = new AliESD; AliESD* hltesd = new AliESD;
+  TFile* fileOld = NULL;
+  TTree* treeOld = NULL; TTree *hlttreeOld = NULL;
+  if (!gSystem->AccessPathName("AliESDs.root")){
+    gSystem->CopyFile("AliESDs.root", "AliESDs.old.root", kTRUE);
+    fileOld = TFile::Open("AliESDs.old.root");
+    if (fileOld && fileOld->IsOpen()) {
+      treeOld = (TTree*) fileOld->Get("esdTree");
+      if (treeOld) treeOld->SetBranchAddress("ESD", &esd);
+      hlttreeOld = (TTree*) fileOld->Get("HLTesdTree");
+      if (hlttreeOld) hlttreeOld->SetBranchAddress("ESD", &hltesd);
+    }
+  }
+
   // create the ESD output file and tree
   TFile* file = TFile::Open("AliESDs.root", "RECREATE");
   if (!file->IsOpen()) {
     AliError("opening AliESDs.root failed");
-    if (fStopOnError) {CleanUp(file); return kFALSE;}    
+    if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}    
   }
-  AliESD* esd = new AliESD;
   TTree* tree = new TTree("esdTree", "Tree with ESD objects");
   tree->Branch("ESD", "AliESD", &esd);
-  delete esd;
+  TTree* hlttree = new TTree("HLTesdTree", "Tree with HLT ESD objects");
+  hlttree->Branch("ESD", "AliESD", &hltesd);
+  delete esd; delete hltesd;
+  esd = NULL; hltesd = NULL;
   gROOT->cd();
 
   // loop over events
   if (fRawReader) fRawReader->RewindEvents();
+  
   for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
+    if (fRawReader) fRawReader->NextEvent();
+    if ((iEvent < firstEvent) || ((lastEvent >= 0) && (iEvent > lastEvent))) {
+      // copy old ESD to the new one
+      if (treeOld) {
+       treeOld->SetBranchAddress("ESD", &esd);
+       treeOld->GetEntry(iEvent);
+      }
+      tree->Fill();
+      if (hlttreeOld) {
+       hlttreeOld->SetBranchAddress("ESD", &hltesd);
+       hlttreeOld->GetEntry(iEvent);
+      }
+      hlttree->Fill();
+      continue;
+    }
+
     AliInfo(Form("processing event %d", iEvent));
     fRunLoader->GetEvent(iEvent);
-    if (fRawReader) fRawReader->NextEvent();
 
     char fileName[256];
     sprintf(fileName, "ESD_%d.%d_final.root", 
-           aliRun->GetRunNumber(), aliRun->GetEvNumber());
+           fRunLoader->GetHeader()->GetRun(), 
+           fRunLoader->GetHeader()->GetEventNrInRun());
     if (!gSystem->AccessPathName(fileName)) continue;
 
-    esd = new AliESD;
-    esd->SetRunNumber(aliRun->GetRunNumber());
-    esd->SetEventNumber(aliRun->GetEvNumber());
-    esd->SetMagneticField(aliRun->Field()->SolenoidField());
+    // local reconstruction
+    if (!fRunLocalReconstruction.IsNull()) {
+      if (!RunLocalEventReconstruction(fRunLocalReconstruction)) {
+       if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
+      }
+    }
+
+    esd = new AliESD; hltesd = new AliESD;
+    esd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
+    hltesd->SetRunNumber(fRunLoader->GetHeader()->GetRun());
+    esd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
+    hltesd->SetEventNumber(fRunLoader->GetHeader()->GetEventNrInRun());
+    if (gAlice) {
+      esd->SetMagneticField(gAlice->Field()->SolenoidField());
+      hltesd->SetMagneticField(gAlice->Field()->SolenoidField());
+    } else {
+      // ???
+    }
 
     // vertex finder
     if (fRunVertexFinder) {
       if (!ReadESD(esd, "vertex")) {
        if (!RunVertexFinder(esd)) {
-         if (fStopOnError) {CleanUp(file); return kFALSE;}
+         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
        }
        if (fCheckPointLevel > 0) WriteESD(esd, "vertex");
       }
     }
 
+    // HLT tracking
+    if (!fRunTracking.IsNull()) {
+      if (fRunHLTTracking) {
+       hltesd->SetVertex(esd->GetVertex());
+       if (!RunHLTTracking(hltesd)) {
+         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
+       }
+      }
+    }
+
     // barrel tracking
     if (!fRunTracking.IsNull()) {
       if (!ReadESD(esd, "tracking")) {
        if (!RunTracking(esd)) {
-         if (fStopOnError) {CleanUp(file); return kFALSE;}
+         if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
        }
        if (fCheckPointLevel > 0) WriteESD(esd, "tracking");
       }
@@ -333,7 +409,7 @@ Bool_t AliReconstruction::Run(const char* input)
     // fill ESD
     if (!fFillESD.IsNull()) {
       if (!FillESD(esd, fFillESD)) {
-       if (fStopOnError) {CleanUp(file); return kFALSE;}
+       if (fStopOnError) {CleanUp(file, fileOld); return kFALSE;}
       }
     }
 
@@ -343,14 +419,23 @@ Bool_t AliReconstruction::Run(const char* input)
 
     // write ESD
     tree->Fill();
+    // write HLT ESD
+    hlttree->Fill();
 
-    if (fCheckPointLevel > 0) WriteESD(esd, "final");
-    delete esd;
+    if (fCheckPointLevel > 0)  WriteESD(esd, "final"); 
+    delete esd; delete hltesd;
+    esd = NULL; hltesd = NULL;
   }
 
   file->cd();
   tree->Write();
-  CleanUp(file);
+  hlttree->Write();
+
+  // Create tags for the events in the ESD tree (the ESD tree is always present)
+  // In case of empty events the tags will contain dummy values
+  CreateTag(file);
+  CleanUp(file, fileOld);
 
   return kTRUE;
 }
@@ -369,6 +454,7 @@ Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
     if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
     AliReconstructor* reconstructor = GetReconstructor(iDet);
     if (!reconstructor) continue;
+    if (reconstructor->HasLocalReconstruction()) continue;
 
     AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
     TStopwatch stopwatchDet;
@@ -379,8 +465,9 @@ Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
     } else {
       reconstructor->Reconstruct(fRunLoader);
     }
-    AliInfo(Form("execution time for %s:", fgkDetectorName[iDet]));
-    ToAliInfo(stopwatchDet.Print());
+    AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
+                fgkDetectorName[iDet],
+                stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
   }
 
   if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
@@ -389,8 +476,82 @@ Bool_t AliReconstruction::RunLocalReconstruction(const TString& detectors)
     if (fStopOnError) return kFALSE;
   }
 
-  AliInfo("execution time:");
-  ToAliInfo(stopwatch.Print());
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
+
+  return kTRUE;
+}
+
+//_____________________________________________________________________________
+Bool_t AliReconstruction::RunLocalEventReconstruction(const TString& detectors)
+{
+// run the local reconstruction
+
+  TStopwatch stopwatch;
+  stopwatch.Start();
+
+  TString detStr = detectors;
+  for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
+    if (!IsSelected(fgkDetectorName[iDet], detStr)) continue;
+    AliReconstructor* reconstructor = GetReconstructor(iDet);
+    if (!reconstructor) continue;
+    AliLoader* loader = fLoader[iDet];
+
+    // conversion of digits
+    if (fRawReader && reconstructor->HasDigitConversion()) {
+      AliInfo(Form("converting raw data digits into root objects for %s", 
+                  fgkDetectorName[iDet]));
+      TStopwatch stopwatchDet;
+      stopwatchDet.Start();
+      loader->LoadDigits("update");
+      loader->CleanDigits();
+      loader->MakeDigitsContainer();
+      TTree* digitsTree = loader->TreeD();
+      reconstructor->ConvertDigits(fRawReader, digitsTree);
+      loader->WriteDigits("OVERWRITE");
+      loader->UnloadDigits();
+      AliInfo(Form("Execution time for %s: R:%.2fs C:%.2fs",
+                  fgkDetectorName[iDet],
+                  stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
+    }
+
+    // local reconstruction
+    if (!reconstructor->HasLocalReconstruction()) continue;
+    AliInfo(Form("running reconstruction for %s", fgkDetectorName[iDet]));
+    TStopwatch stopwatchDet;
+    stopwatchDet.Start();
+    loader->LoadRecPoints("update");
+    loader->CleanRecPoints();
+    loader->MakeRecPointsContainer();
+    TTree* clustersTree = loader->TreeR();
+    if (fRawReader && !reconstructor->HasDigitConversion()) {
+      reconstructor->Reconstruct(fRawReader, clustersTree);
+    } else {
+      loader->LoadDigits("read");
+      TTree* digitsTree = loader->TreeD();
+      if (!digitsTree) {
+       AliError(Form("Can't get the %s digits tree", fgkDetectorName[iDet]));
+       if (fStopOnError) return kFALSE;
+      } else {
+       reconstructor->Reconstruct(digitsTree, clustersTree);
+      }
+      loader->UnloadDigits();
+    }
+    loader->WriteRecPoints("OVERWRITE");
+    loader->UnloadRecPoints();
+    AliDebug(1,Form("Execution time for %s: R:%.2fs C:%.2fs",
+                   fgkDetectorName[iDet],
+                   stopwatchDet.RealTime(),stopwatchDet.CpuTime()));
+  }
+
+  if ((detStr.CompareTo("ALL") != 0) && !detStr.IsNull()) {
+    AliError(Form("the following detectors were not found: %s",
+                  detStr.Data()));
+    if (fStopOnError) return kFALSE;
+  }
+  
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
 
   return kTRUE;
 }
@@ -414,7 +575,9 @@ Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
 
   if (fVertexer) {
     AliInfo("running the ITS vertex finder");
+    if (fLoader[0]) fLoader[0]->LoadRecPoints();
     vertex = fVertexer->FindVertexForCurrentEvent(fRunLoader->GetEventNumber());
+    if (fLoader[0]) fLoader[0]->UnloadRecPoints();
     if(!vertex){
       AliWarning("Vertex not found");
       vertex = new AliESDVertex();
@@ -441,8 +604,68 @@ Bool_t AliReconstruction::RunVertexFinder(AliESD*& esd)
   }  
   delete vertex;
 
-  AliInfo("execution time:");
-  ToAliInfo(stopwatch.Print());
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
+
+  return kTRUE;
+}
+
+//_____________________________________________________________________________
+Bool_t AliReconstruction::RunHLTTracking(AliESD*& esd)
+{
+// run the HLT barrel tracking
+
+  TStopwatch stopwatch;
+  stopwatch.Start();
+
+  if (!fRunLoader) {
+    AliError("Missing runLoader!");
+    return kFALSE;
+  }
+
+  AliInfo("running HLT tracking");
+
+  // Get a pointer to the HLT reconstructor
+  AliReconstructor *reconstructor = GetReconstructor(fgkNDetectors-1);
+  if (!reconstructor) return kFALSE;
+
+  // TPC + ITS
+  for (Int_t iDet = 1; iDet >= 0; iDet--) {
+    TString detName = fgkDetectorName[iDet];
+    AliDebug(1, Form("%s HLT tracking", detName.Data()));
+    reconstructor->SetOption(detName.Data());
+    AliTracker *tracker = reconstructor->CreateTracker(fRunLoader);
+    if (!tracker) {
+      AliWarning(Form("couldn't create a HLT tracker for %s", detName.Data()));
+      if (fStopOnError) return kFALSE;
+      continue;
+    }
+    Double_t vtxPos[3];
+    Double_t vtxErr[3]={0.005,0.005,0.010};
+    const AliESDVertex *vertex = esd->GetVertex();
+    vertex->GetXYZ(vtxPos);
+    tracker->SetVertex(vtxPos,vtxErr);
+    if(iDet != 1) {
+      fLoader[iDet]->LoadRecPoints("read");
+      TTree* tree = fLoader[iDet]->TreeR();
+      if (!tree) {
+       AliError(Form("Can't get the %s cluster tree", detName.Data()));
+       return kFALSE;
+      }
+      tracker->LoadClusters(tree);
+    }
+    if (tracker->Clusters2Tracks(esd) != 0) {
+      AliError(Form("HLT %s Clusters2Tracks failed", fgkDetectorName[iDet]));
+      return kFALSE;
+    }
+    if(iDet != 1) {
+      tracker->UnloadClusters();
+    }
+    delete tracker;
+  }
+
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
 
   return kTRUE;
 }
@@ -479,6 +702,12 @@ Bool_t AliReconstruction::RunTracking(AliESD*& esd)
     if (fCheckPointLevel > 1) {
       WriteESD(esd, Form("%s.tracking", fgkDetectorName[iDet]));
     }
+    // preliminary PID in TPC needed by the ITS tracker
+    if (iDet == 1) {
+      GetReconstructor(1)->FillESD(fRunLoader, esd);
+      GetReconstructor(1)->FillESD((TTree*)NULL, (TTree*)NULL, esd);
+      AliESDpid::MakePID(esd);
+    }
   }
 
   // pass 2: ALL backwards
@@ -489,13 +718,8 @@ Bool_t AliReconstruction::RunTracking(AliESD*& esd)
     // load clusters
     if (iDet > 1) {     // all except ITS, TPC
       TTree* tree = NULL;
-      if (iDet == 3) {   // TOF
-       fLoader[iDet]->LoadDigits("read");
-       tree = fLoader[iDet]->TreeD();
-      } else {
-       fLoader[iDet]->LoadRecPoints("read");
-       tree = fLoader[iDet]->TreeR();
-      }
+      fLoader[iDet]->LoadRecPoints("read");
+      tree = fLoader[iDet]->TreeR();
       if (!tree) {
        AliError(Form("Can't get the %s cluster tree", fgkDetectorName[iDet]));
        return kFALSE;
@@ -515,11 +739,7 @@ Bool_t AliReconstruction::RunTracking(AliESD*& esd)
     // unload clusters
     if (iDet > 2) {     // all except ITS, TPC, TRD
       fTracker[iDet]->UnloadClusters();
-      if (iDet == 3) {   // TOF
-       fLoader[iDet]->UnloadDigits();
-      } else {
-       fLoader[iDet]->UnloadRecPoints();
-      }
+      fLoader[iDet]->UnloadRecPoints();
     }
   }
 
@@ -542,8 +762,8 @@ Bool_t AliReconstruction::RunTracking(AliESD*& esd)
     fLoader[iDet]->UnloadRecPoints();
   }
 
-  AliInfo("execution time:");
-  ToAliInfo(stopwatch.Print());
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
 
   return kTRUE;
 }
@@ -565,6 +785,36 @@ Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
 
     if (!ReadESD(esd, fgkDetectorName[iDet])) {
       AliDebug(1, Form("filling ESD for %s", fgkDetectorName[iDet]));
+      TTree* clustersTree = NULL;
+      if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
+       fLoader[iDet]->LoadRecPoints("read");
+       clustersTree = fLoader[iDet]->TreeR();
+       if (!clustersTree) {
+         AliError(Form("Can't get the %s clusters tree", 
+                       fgkDetectorName[iDet]));
+         if (fStopOnError) return kFALSE;
+       }
+      }
+      if (fRawReader && !reconstructor->HasDigitConversion()) {
+        reconstructor->FillESD(fRawReader, clustersTree, esd);
+      } else {
+       TTree* digitsTree = NULL;
+       if (fLoader[iDet]) {
+         fLoader[iDet]->LoadDigits("read");
+         digitsTree = fLoader[iDet]->TreeD();
+         if (!digitsTree) {
+           AliError(Form("Can't get the %s digits tree", 
+                         fgkDetectorName[iDet]));
+           if (fStopOnError) return kFALSE;
+         }
+       }
+       reconstructor->FillESD(digitsTree, clustersTree, esd);
+       if (fLoader[iDet]) fLoader[iDet]->UnloadDigits();
+      }
+      if (reconstructor->HasLocalReconstruction() && fLoader[iDet]) {
+       fLoader[iDet]->UnloadRecPoints();
+      }
+
       if (fRawReader) {
         reconstructor->FillESD(fRunLoader, fRawReader, esd);
       } else {
@@ -580,8 +830,8 @@ Bool_t AliReconstruction::FillESD(AliESD*& esd, const TString& detectors)
     if (fStopOnError) return kFALSE;
   }
 
-  AliInfo("execution time:");
-  ToAliInfo(stopwatch.Print());
+  AliInfo(Form("Execution time: R:%.2fs C:%.2fs",
+              stopwatch.RealTime(),stopwatch.CpuTime()));
 
   return kTRUE;
 }
@@ -620,10 +870,87 @@ Bool_t AliReconstruction::IsSelected(TString detName, TString& detectors) const
   return result;
 }
 
+//_____________________________________________________________________________
+Bool_t AliReconstruction::InitRunLoader()
+{
+// get or create the run loader
+
+  if (gAlice) delete gAlice;
+  gAlice = NULL;
+
+  if (!gSystem->AccessPathName(fGAliceFileName.Data())) { // galice.root exists
+    // load all base libraries to get the loader classes
+    TString libs = gSystem->GetLibraries();
+    for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
+      TString detName = fgkDetectorName[iDet];
+      if (detName == "HLT") continue;
+      if (libs.Contains("lib" + detName + "base.so")) continue;
+      gSystem->Load("lib" + detName + "base.so");
+    }
+    fRunLoader = AliRunLoader::Open(fGAliceFileName.Data());
+    if (!fRunLoader) {
+      AliError(Form("no run loader found in file %s", fGAliceFileName.Data()));
+      CleanUp();
+      return kFALSE;
+    }
+    fRunLoader->CdGAFile();
+    if (gFile->GetKey(AliRunLoader::GetGAliceName())) {
+      if (fRunLoader->LoadgAlice() == 0) {
+       gAlice = fRunLoader->GetAliRun();
+       AliTracker::SetFieldMap(gAlice->Field(),fUniformField);
+       AliExternalTrackParam::SetFieldMap(gAlice->Field());
+       if(fUniformField)
+         AliExternalTrackParam::SetUniformFieldTracking();
+       else
+         AliExternalTrackParam::SetNonuniformFieldTracking();
+      }
+    }
+    if (!gAlice && !fRawReader) {
+      AliError(Form("no gAlice object found in file %s",
+                   fGAliceFileName.Data()));
+      CleanUp();
+      return kFALSE;
+    }
+
+  } else {               // galice.root does not exist
+    if (!fRawReader) {
+      AliError(Form("the file %s does not exist", fGAliceFileName.Data()));
+      CleanUp();
+      return kFALSE;
+    }
+    fRunLoader = AliRunLoader::Open(fGAliceFileName.Data(),
+                                   AliConfig::GetDefaultEventFolderName(),
+                                   "recreate");
+    if (!fRunLoader) {
+      AliError(Form("could not create run loader in file %s", 
+                   fGAliceFileName.Data()));
+      CleanUp();
+      return kFALSE;
+    }
+    fRunLoader->MakeTree("E");
+    Int_t iEvent = 0;
+    while (fRawReader->NextEvent()) {
+      fRunLoader->SetEventNumber(iEvent);
+      fRunLoader->GetHeader()->Reset(fRawReader->GetRunNumber(), 
+                                    iEvent, iEvent);
+      fRunLoader->MakeTree("H");
+      fRunLoader->TreeE()->Fill();
+      iEvent++;
+    }
+    fRawReader->RewindEvents();
+    fRunLoader->WriteHeader("OVERWRITE");
+    fRunLoader->CdGAFile();
+    fRunLoader->Write(0, TObject::kOverwrite);
+//    AliTracker::SetFieldMap(???);
+  }
+
+  return kTRUE;
+}
+
 //_____________________________________________________________________________
 AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
 {
-// get the reconstructor object for a detector
+// get the reconstructor object and the loader for a detector
 
   if (fReconstructor[iDet]) return fReconstructor[iDet];
 
@@ -631,7 +958,7 @@ AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
   TPluginManager* pluginManager = gROOT->GetPluginManager();
   TString detName = fgkDetectorName[iDet];
   TString recName = "Ali" + detName + "Reconstructor";
-  if (!gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
+  if (gAlice && !gAlice->GetDetector(detName) && (detName != "HLT")) return NULL;
 
   if (detName == "HLT") {
     if (!gROOT->GetClass("AliLevel3")) {
@@ -646,10 +973,12 @@ AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
   // first check if a plugin is defined for the reconstructor
   TPluginHandler* pluginHandler = 
     pluginManager->FindHandler("AliReconstructor", detName);
-  // if not, but the reconstructor class is implemented, add a plugin for it
-  if (!pluginHandler && gROOT->GetClass(recName.Data())) {
+  // if not, add a plugin for it
+  if (!pluginHandler) {
     AliDebug(1, Form("defining plugin for %s", recName.Data()));
-    if (gSystem->Load("lib" + detName + "base.so") == 0) {
+    TString libs = gSystem->GetLibraries();
+    if (libs.Contains("lib" + detName + "base.so") ||
+       (gSystem->Load("lib" + detName + "base.so") >= 0)) {
       pluginManager->AddHandler("AliReconstructor", detName, 
                                recName, detName + "rec", recName + "()");
     } else {
@@ -664,9 +993,49 @@ AliReconstructor* AliReconstruction::GetReconstructor(Int_t iDet)
   if (reconstructor) {
     TObject* obj = fOptions.FindObject(detName.Data());
     if (obj) reconstructor->SetOption(obj->GetTitle());
+    reconstructor->Init(fRunLoader);
     fReconstructor[iDet] = reconstructor;
   }
 
+  // get or create the loader
+  if (detName != "HLT") {
+    fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
+    if (!fLoader[iDet]) {
+      AliConfig::Instance()
+       ->CreateDetectorFolders(fRunLoader->GetEventFolder(), 
+                               detName, detName);
+      // first check if a plugin is defined for the loader
+      TPluginHandler* pluginHandler = 
+       pluginManager->FindHandler("AliLoader", detName);
+      // if not, add a plugin for it
+      if (!pluginHandler) {
+       TString loaderName = "Ali" + detName + "Loader";
+       AliDebug(1, Form("defining plugin for %s", loaderName.Data()));
+       pluginManager->AddHandler("AliLoader", detName, 
+                                 loaderName, detName + "base", 
+                                 loaderName + "(const char*, TFolder*)");
+       pluginHandler = pluginManager->FindHandler("AliLoader", detName);
+      }
+      if (pluginHandler && (pluginHandler->LoadPlugin() == 0)) {
+       fLoader[iDet] = 
+         (AliLoader*) pluginHandler->ExecPlugin(2, detName.Data(), 
+                                                fRunLoader->GetEventFolder());
+      }
+      if (!fLoader[iDet]) {   // use default loader
+       fLoader[iDet] = new AliLoader(detName, fRunLoader->GetEventFolder());
+      }
+      if (!fLoader[iDet]) {
+       AliWarning(Form("couldn't get loader for %s", detName.Data()));
+       if (fStopOnError) return NULL;
+      } else {
+       fRunLoader->AddLoader(fLoader[iDet]);
+       fRunLoader->CdGAFile();
+       if (gFile && !gFile->IsWritable()) gFile->ReOpen("UPDATE");
+       fRunLoader->Write(0, TObject::kOverwrite);
+      }
+    }
+  }
+      
   return reconstructor;
 }
 
@@ -691,7 +1060,7 @@ Bool_t AliReconstruction::CreateVertexer()
 //_____________________________________________________________________________
 Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
 {
-// get the loaders and create the trackers
+// create the trackers
 
   TString detStr = detectors;
   for (Int_t iDet = 0; iDet < fgkNDetectors; iDet++) {
@@ -699,17 +1068,15 @@ Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
     AliReconstructor* reconstructor = GetReconstructor(iDet);
     if (!reconstructor) continue;
     TString detName = fgkDetectorName[iDet];
-    if (detName == "HLT") continue;
-    fLoader[iDet] = fRunLoader->GetLoader(detName + "Loader");
-    if (!fLoader[iDet]) {
-      AliWarning(Form("no %s loader found", detName.Data()));
+    if (detName == "HLT") {
+      fRunHLTTracking = kTRUE;
+      continue;
+    }
+
+    fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
+    if (!fTracker[iDet] && (iDet < 7)) {
+      AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
       if (fStopOnError) return kFALSE;
-    } else {
-      fTracker[iDet] = reconstructor->CreateTracker(fRunLoader);
-      if (!fTracker[iDet] && (iDet < 7)) {
-       AliWarning(Form("couldn't create a tracker for %s", detName.Data()));
-       if (fStopOnError) return kFALSE;
-      }
     }
   }
 
@@ -717,7 +1084,7 @@ Bool_t AliReconstruction::CreateTrackers(const TString& detectors)
 }
 
 //_____________________________________________________________________________
-void AliReconstruction::CleanUp(TFile* file)
+void AliReconstruction::CleanUp(TFile* file, TFile* fileOld)
 {
 // delete trackers and the run loader and close and delete the file
 
@@ -740,6 +1107,12 @@ void AliReconstruction::CleanUp(TFile* file)
     file->Close();
     delete file;
   }
+
+  if (fileOld) {
+    fileOld->Close();
+    delete fileOld;
+    gSystem->Unlink("AliESDs.old.root");
+  }
 }
 
 
@@ -754,6 +1127,7 @@ Bool_t AliReconstruction::ReadESD(AliESD*& esd, const char* recStep) const
          esd->GetRunNumber(), esd->GetEventNumber(), recStep);
   if (gSystem->AccessPathName(fileName)) return kFALSE;
 
+  AliInfo(Form("reading ESD from file %s", fileName));
   AliDebug(1, Form("reading ESD from file %s", fileName));
   TFile* file = TFile::Open(fileName);
   if (!file || !file->IsOpen()) {
@@ -790,3 +1164,251 @@ void AliReconstruction::WriteESD(AliESD* esd, const char* recStep) const
   }
   delete file;
 }
+
+
+
+
+//_____________________________________________________________________________
+void AliReconstruction::CreateTag(TFile* file)
+{
+  // Creates the tags for all the events in a given ESD file
+  Int_t ntrack;
+  Int_t nProtons, nKaons, nPions, nMuons, nElectrons;
+  Int_t nPos, nNeg, nNeutr;
+  Int_t nK0s, nNeutrons, nPi0s, nGamas;
+  Int_t nCh1GeV, nCh3GeV, nCh10GeV;
+  Int_t nMu1GeV, nMu3GeV, nMu10GeV;
+  Int_t nEl1GeV, nEl3GeV, nEl10GeV;
+  Float_t maxPt = .0, meanPt = .0, totalP = .0;
+
+  AliRunTag *tag = new AliRunTag();
+  AliDetectorTag *detTag = new AliDetectorTag();
+  AliEventTag *evTag = new AliEventTag();
+  TTree ttag("T","A Tree with event tags");
+  TBranch * btag = ttag.Branch("AliTAG", "AliRunTag", &tag);
+  btag->SetCompressionLevel(9);
+
+  AliInfo(Form("Creating the tags......."));   
+  
+  if (!file || !file->IsOpen()) {
+    AliError(Form("opening failed"));
+    delete file;
+    return ;
+  }
+
+  TTree *t = (TTree*) file->Get("esdTree");
+  TBranch * b = t->GetBranch("ESD");
+  AliESD *esd = 0;
+  b->SetAddress(&esd);
+
+  tag->SetRunId(esd->GetRunNumber());
+
+  Int_t firstEvent = 0,lastEvent = 0;
+  Int_t iNumberOfEvents = b->GetEntries();
+  for (Int_t iEventNumber = 0; iEventNumber < iNumberOfEvents; iEventNumber++)
+    {
+      ntrack = 0;
+      nPos = 0;
+      nNeg = 0;
+      nNeutr =0;
+      nK0s = 0;
+      nNeutrons = 0;
+      nPi0s = 0;
+      nGamas = 0;
+      nProtons = 0;
+      nKaons = 0;
+      nPions = 0;
+      nMuons = 0;
+      nElectrons = 0;    
+      nCh1GeV = 0;
+      nCh3GeV = 0;
+      nCh10GeV = 0;
+      nMu1GeV = 0;
+      nMu3GeV = 0;
+      nMu10GeV = 0;
+      nEl1GeV = 0;
+      nEl3GeV = 0;
+      nEl10GeV = 0;
+      maxPt = .0;
+      meanPt = .0;
+      totalP = .0;
+
+      b->GetEntry(iEventNumber);
+      const AliESDVertex * vertexIn = esd->GetVertex();
+
+      for (Int_t iTrackNumber = 0; iTrackNumber < esd->GetNumberOfTracks(); iTrackNumber++)
+       {
+         AliESDtrack * esdTrack = esd->GetTrack(iTrackNumber);
+         UInt_t status = esdTrack->GetStatus();
+         
+         //select only tracks with ITS refit
+         if ((status&AliESDtrack::kITSrefit)==0) continue;
+         
+         //select only tracks with TPC refit-->remove extremely high Pt tracks
+         if ((status&AliESDtrack::kTPCrefit)==0) continue;
+         
+         //select only tracks with the "combined PID"
+         if ((status&AliESDtrack::kESDpid)==0) continue;
+                 Double_t p[3];
+         esdTrack->GetPxPyPz(p);
+         Double_t momentum = sqrt(pow(p[0],2) + pow(p[1],2) + pow(p[2],2));
+         Double_t fPt = sqrt(pow(p[0],2) + pow(p[1],2));
+         totalP += momentum;
+         meanPt += fPt;
+         if(fPt > maxPt)
+           maxPt = fPt;
+         
+         if(esdTrack->GetSign() > 0)
+           {
+             nPos++;
+             if(fPt > 1.0)
+               nCh1GeV++;
+             if(fPt > 3.0)
+               nCh3GeV++;
+             if(fPt > 10.0)
+               nCh10GeV++;
+           }
+         if(esdTrack->GetSign() < 0)
+           {
+             nNeg++;
+             if(fPt > 1.0)
+               nCh1GeV++;
+             if(fPt > 3.0)
+               nCh3GeV++;
+             if(fPt > 10.0)
+               nCh10GeV++;
+           }
+         if(esdTrack->GetSign() == 0)
+           nNeutr++;
+         
+         //PID
+         Double_t prob[10];
+         esdTrack->GetESDpid(prob);
+                   
+         //K0s
+         if ((prob[8]>prob[7])&&(prob[8]>prob[6])&&(prob[8]>prob[5])&&(prob[8]>prob[4])&&(prob[8]>prob[3])&&(prob[8]>prob[2])&&(prob[8]>prob[1])&&(prob[8]>prob[0]))
+           nK0s++;
+         //neutrons
+         if ((prob[7]>prob[8])&&(prob[7]>prob[6])&&(prob[7]>prob[5])&&(prob[7]>prob[4])&&(prob[7]>prob[3])&&(prob[7]>prob[2])&&(prob[7]>prob[1])&&(prob[7]>prob[0]))
+           nNeutrons++; 
+         //pi0s
+         if ((prob[6]>prob[8])&&(prob[6]>prob[7])&&(prob[6]>prob[5])&&(prob[6]>prob[4])&&(prob[6]>prob[3])&&(prob[6]>prob[2])&&(prob[6]>prob[1])&&(prob[6]>prob[0]))
+           nPi0s++;
+         //gamas
+         if ((prob[5]>prob[8])&&(prob[5]>prob[7])&&(prob[5]>prob[6])&&(prob[5]>prob[4])&&(prob[5]>prob[3])&&(prob[5]>prob[2])&&(prob[5]>prob[1])&&(prob[5]>prob[0]))
+           nGamas++;
+         //protons
+         if ((prob[4]>prob[8])&&(prob[4]>prob[7])&&(prob[4]>prob[6])&&(prob[4]>prob[5])&&(prob[4]>prob[3])&&(prob[4]>prob[2])&&(prob[4]>prob[1])&&(prob[4]>prob[0]))
+           nProtons++;
+         //kaons
+         if ((prob[3]>prob[8])&&(prob[3]>prob[7])&&(prob[3]>prob[6])&&(prob[3]>prob[5])&&(prob[3]>prob[4])&&(prob[3]>prob[2])&&(prob[3]>prob[1])&&(prob[3]>prob[0]))
+           nKaons++;
+         //kaons
+         if ((prob[2]>prob[8])&&(prob[2]>prob[7])&&(prob[2]>prob[6])&&(prob[2]>prob[5])&&(prob[2]>prob[4])&&(prob[2]>prob[3])&&(prob[2]>prob[1])&&(prob[2]>prob[0]))
+           nPions++; 
+         //muons
+         if ((prob[1]>prob[8])&&(prob[1]>prob[7])&&(prob[1]>prob[6])&&(prob[1]>prob[5])&&(prob[1]>prob[4])&&(prob[1]>prob[3])&&(prob[1]>prob[2])&&(prob[1]>prob[0]))
+           {
+             nMuons++;
+             if(fPt > 1.0)
+               nMu1GeV++;
+             if(fPt > 3.0)
+               nMu3GeV++;
+             if(fPt > 10.0)
+               nMu10GeV++;
+           }
+         //electrons
+         if ((prob[0]>prob[8])&&(prob[0]>prob[7])&&(prob[0]>prob[6])&&(prob[0]>prob[5])&&(prob[0]>prob[4])&&(prob[0]>prob[3])&&(prob[0]>prob[2])&&(prob[0]>prob[1]))
+           {
+             nElectrons++;
+             if(fPt > 1.0)
+               nEl1GeV++;
+             if(fPt > 3.0)
+               nEl3GeV++;
+             if(fPt > 10.0)
+               nEl10GeV++;
+           }
+         
+         
+         
+         ntrack++;
+       }//track loop
+      // Fill the event tags 
+      if(ntrack != 0)
+       meanPt = meanPt/ntrack;
+      
+      evTag->SetEventId(iEventNumber+1);
+      evTag->SetVertexX(vertexIn->GetXv());
+      evTag->SetVertexY(vertexIn->GetYv());
+      evTag->SetVertexZ(vertexIn->GetZv());
+      
+      evTag->SetT0VertexZ(esd->GetT0zVertex());
+      
+      evTag->SetTrigger(esd->GetTrigger());
+      
+      evTag->SetZDCNeutronEnergy(esd->GetZDCNEnergy());
+      evTag->SetZDCProtonEnergy(esd->GetZDCPEnergy());
+      evTag->SetZDCEMEnergy(esd->GetZDCEMEnergy());
+      evTag->SetNumOfParticipants(esd->GetZDCParticipants());
+      
+      
+      evTag->SetNumOfTracks(esd->GetNumberOfTracks());
+      evTag->SetNumOfPosTracks(nPos);
+      evTag->SetNumOfNegTracks(nNeg);
+      evTag->SetNumOfNeutrTracks(nNeutr);
+      
+      evTag->SetNumOfV0s(esd->GetNumberOfV0s());
+      evTag->SetNumOfCascades(esd->GetNumberOfCascades());
+      evTag->SetNumOfKinks(esd->GetNumberOfKinks());
+      evTag->SetNumOfPMDTracks(esd->GetNumberOfPmdTracks());
+      
+      evTag->SetNumOfProtons(nProtons);
+      evTag->SetNumOfKaons(nKaons);
+      evTag->SetNumOfPions(nPions);
+      evTag->SetNumOfMuons(nMuons);
+      evTag->SetNumOfElectrons(nElectrons);
+      evTag->SetNumOfPhotons(nGamas);
+      evTag->SetNumOfPi0s(nPi0s);
+      evTag->SetNumOfNeutrons(nNeutrons);
+      evTag->SetNumOfKaon0s(nK0s);
+      
+      evTag->SetNumOfChargedAbove1GeV(nCh1GeV);
+      evTag->SetNumOfChargedAbove3GeV(nCh3GeV);
+      evTag->SetNumOfChargedAbove10GeV(nCh10GeV);
+      evTag->SetNumOfMuonsAbove1GeV(nMu1GeV);
+      evTag->SetNumOfMuonsAbove3GeV(nMu3GeV);
+      evTag->SetNumOfMuonsAbove10GeV(nMu10GeV);
+      evTag->SetNumOfElectronsAbove1GeV(nEl1GeV);
+      evTag->SetNumOfElectronsAbove3GeV(nEl3GeV);
+      evTag->SetNumOfElectronsAbove10GeV(nEl10GeV);
+      
+      evTag->SetNumOfPHOSTracks(esd->GetNumberOfPHOSParticles());
+      evTag->SetNumOfEMCALTracks(esd->GetNumberOfEMCALParticles());
+      
+      evTag->SetTotalMomentum(totalP);
+      evTag->SetMeanPt(meanPt);
+      evTag->SetMaxPt(maxPt);
+  
+      tag->AddEventTag(*evTag);
+    }
+  lastEvent = iNumberOfEvents;
+       
+  ttag.Fill();
+  tag->Clear();
+
+  char fileName[256];
+  sprintf(fileName, "Run%d.Event%d_%d.ESD.tag.root", 
+         tag->GetRunId(),firstEvent,lastEvent );
+  AliInfo(Form("writing tags to file %s", fileName));
+  AliDebug(1, Form("writing tags to file %s", fileName));
+  TFile* ftag = TFile::Open(fileName, "recreate");
+  ftag->cd();
+  ttag.Write();
+  ftag->Close();
+  file->cd();
+  delete tag;
+  delete detTag;
+  delete evTag;
+}
+