]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliRun.cxx
Cleaning includes
[u/mrichter/AliRoot.git] / STEER / AliRun.cxx
index 508cc8cfc9e8e68bdd909f8716a62d5f5e49f90d..08caf547f44b22e7ccc1e62333d735f187c845f4 100644 (file)
 
 /*
 $Log$
+Revision 1.91  2002/10/29 14:26:49  hristov
+Code clean-up (F.Carminati)
+
+Revision 1.90  2002/10/22 15:02:15  alibrary
+Introducing Riostream.h
+
+Revision 1.89  2002/10/17 16:26:39  hristov
+Definition of additional particles moved to VMC (I.Hrivnacova)
+
+Revision 1.88  2002/10/14 14:57:32  hristov
+Merging the VirtualMC branch to the main development branch (HEAD)
+
+Revision 1.82.4.3  2002/10/14 09:45:57  hristov
+Updating VirtualMC to v3-09-02
+
+Revision 1.82.4.2  2002/06/10 17:54:06  hristov
+Only one SetGenEventHeader function kept
+
+Revision 1.82.4.1  2002/06/10 14:43:06  hristov
+Merged with v3-08-02
+
+Revision 1.86  2002/05/28 14:24:57  hristov
+Correct warning messages
+Revision 1.85  2002/05/24 13:29:58  hristov
+AliTrackReference added, AliDisplay modified
+
+Revision 1.84  2002/05/21 16:26:07  hristov
+Find correctly TreeK in case CONFIG_SPLIT_FILE is set (Y.Schutz)
+
+Revision 1.83  2002/04/04 13:16:17  jchudoba
+add possibility to write sdigits, digits and rec. points into separate files
+
+Revision 1.82  2002/03/12 11:06:03  morsch
+Add particle status code to argument list of SetTrack(..).
+
+Revision 1.81  2001/12/19 14:46:26  morsch
+Add possibility to disable StepManager() for each module separately.
+
+Revision 1.80  2001/10/21 18:22:55  hristov
+BranchOld replaced by Branch. It works correctly with Root 2.02.xx
+
+Revision 1.79  2001/10/09 18:00:35  hristov
+Temporary fix to provide unique event number in the simulation (J.Chudoba)
+
+Revision 1.78  2001/10/04 15:30:56  hristov
+Changes to accommodate the set of PHOS folders and tasks (Y.Schutz)
+
+Revision 1.77  2001/09/04 15:09:11  hristov
+fTreeE->Branch replaced temporary by fTreeE->BranchOld to avoid data corruption in case of many events per file
+
+Revision 1.76  2001/08/03 14:38:35  morsch
+Use original method to access TreeH.
+
+Revision 1.75  2001/07/28 10:52:27  hristov
+Event number updated correctly (M.Ivanov)
+
+Revision 1.74  2001/07/28 10:39:16  morsch
+GetEventsPerRun() method needed by afterburners added to AliRun.h
+Corresponding setters and getters have been from AliGenerator.
+
+Revision 1.73  2001/07/27 12:34:30  jchudoba
+remove the dummy argument in fStack->GetEvent call
+
+Revision 1.72  2001/07/03 08:10:57  hristov
+J.Chudoba's changes merged correctly with the HEAD
+
+Revision 1.70  2001/06/29 08:01:36  morsch
+Small correction to the previous.
+
+Revision 1.69  2001/06/28 16:27:50  morsch
+AliReco() with user control of event range.
+
+Revision 1.68  2001/06/11 13:14:40  morsch
+SetAliGenEventHeader() method added.
+
+Revision 1.67  2001/06/07 18:24:50  buncic
+Removed compilation warning in AliConfig initialisation.
+
 Revision 1.66  2001/05/22 14:32:40  hristov
 Weird inline removed
 
@@ -189,12 +268,10 @@ Introduction of the Copyright and cvs Log
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
-#include <iostream.h>
+#include <Riostream.h>
  
 #include <TFile.h>
-#include <TRandom.h>
 #include <TBRIK.h> 
-#include <TNode.h> 
 #include <TCint.h> 
 #include <TSystem.h>
 #include <TObjectTable.h>
@@ -203,7 +280,8 @@ Introduction of the Copyright and cvs Log
 #include <TROOT.h>
 #include <TBrowser.h>
 #include <TFolder.h>
-
+#include <TKey.h>
+#include <TNode.h>
 #include "TParticle.h"
 #include "AliRun.h"
 #include "AliDisplay.h"
@@ -221,53 +299,160 @@ Introduction of the Copyright and cvs Log
 #include "AliStack.h"
 #include "AliGenEventHeader.h"
 #include "AliHeader.h"
-
 #include "AliDetector.h"
+#include "AliPDG.h"
 
 AliRun *gAlice;
 
 ClassImp(AliRun)
 
-//_____________________________________________________________________________
-AliRun::AliRun()
+//_______________________________________________________________________
+AliRun::AliRun():
+  fRun(0),
+  fEvent(0),
+  fEventNrInRun(0),
+  fEventsPerRun(0),
+  fDebug(0),
+  fHeader(0),
+  fTreeD(0),
+  fTreeS(0),
+  fTreeH(0),
+  fTreeTR(0),
+  fTreeE(0),
+  fTreeR(0),
+  fModules(0),
+  fGeometry(0),
+  fDisplay(0),
+  fTimer(),
+  fField(0),
+  fMC(0),
+  fImedia(0),
+  fNdets(0),
+  fTrRmax(1.e10),
+  fTrZmax(1.e10),
+  fGenerator(0),
+  fInitDone(kFALSE),
+  fLego(0),
+  fPDGDB(0),  //Particle factory object
+  fHitLists(0),
+  fEventEnergy(0),
+  fSummEnergy(0),
+  fSum2Energy(0),
+  fConfigFunction("\0"),
+  fRandom(0),
+  fMCQA(0),
+  fTransParName("\0"),
+  fBaseFileName("\0"),
+  fStack(0),
+  fTreeDFileName(""),
+  fTreeDFile(0),
+  fTreeSFileName(""),
+  fTreeSFile(0),
+  fTreeRFileName(""),
+  fTreeRFile(0)
 {
   //
   // Default constructor for AliRun
   //
-  fHeader    = 0;
-  fRun       = 0;
-  fEvent     = 0;
-  fStack     = 0;
-  fModules   = 0;
-  fGenerator = 0;
-  fTreeD     = 0;
-  fTreeH     = 0;
-  fTreeE     = 0;
-  fTreeR     = 0;
-  fTreeS     = 0;
-  fGeometry  = 0;
-  fDisplay   = 0;
-  fField     = 0;
-  fMC        = 0;
-  fNdets     = 0;
-  fImedia    = 0;
-  fTrRmax    = 1.e10;
-  fTrZmax    = 1.e10;
-  fInitDone  = kFALSE;
-  fLego      = 0;
-  fPDGDB     = 0;        //Particle factory object!
-  fHitLists  = 0;
-  fConfigFunction    = "\0";
-  fRandom = 0;
-  fMCQA = 0;
-  fTransParName = "\0";
-  fBaseFileName = ".\0";
-  fDebug        = 0;
+}
+
+//_______________________________________________________________________
+AliRun::AliRun(const AliRun& arun):
+  TVirtualMCApplication(arun),
+  fRun(0),
+  fEvent(0),
+  fEventNrInRun(0),
+  fEventsPerRun(0),
+  fDebug(0),
+  fHeader(0),
+  fTreeD(0),
+  fTreeS(0),
+  fTreeH(0),
+  fTreeTR(0),
+  fTreeE(0),
+  fTreeR(0),
+  fModules(0),
+  fGeometry(0),
+  fDisplay(0),
+  fTimer(),
+  fField(0),
+  fMC(0),
+  fImedia(0),
+  fNdets(0),
+  fTrRmax(1.e10),
+  fTrZmax(1.e10),
+  fGenerator(0),
+  fInitDone(kFALSE),
+  fLego(0),
+  fPDGDB(0),  //Particle factory object
+  fHitLists(0),
+  fEventEnergy(0),
+  fSummEnergy(0),
+  fSum2Energy(0),
+  fConfigFunction("\0"),
+  fRandom(0),
+  fMCQA(0),
+  fTransParName("\0"),
+  fBaseFileName("\0"),
+  fStack(0),
+  fTreeDFileName(""),
+  fTreeDFile(0),
+  fTreeSFileName(""),
+  fTreeSFile(0),
+  fTreeRFileName(""),
+  fTreeRFile(0)
+{
+  //
+  // Copy constructor for AliRun
+  //
+  arun.Copy(*this);
 }
 
 //_____________________________________________________________________________
-AliRun::AliRun(const char *name, const char *title)
-  : TNamed(name,title)
+AliRun::AliRun(const char *name, const char *title):
+  TVirtualMCApplication(name,title),
+  fRun(0),
+  fEvent(0),
+  fEventNrInRun(0),
+  fEventsPerRun(0),
+  fDebug(0),
+  fHeader(new AliHeader()),
+  fTreeD(0),
+  fTreeS(0),
+  fTreeH(0),
+  fTreeTR(0),
+  fTreeE(0),
+  fTreeR(0),
+  fModules(new TObjArray(77)), // Support list for the Detectors
+  fGeometry(0),
+  fDisplay(0),
+  fTimer(),
+  fField(0),
+  fMC(gMC),
+  fImedia(new TArrayI(1000)),
+  fNdets(0),
+  fTrRmax(1.e10),
+  fTrZmax(1.e10),
+  fGenerator(0),
+  fInitDone(kFALSE),
+  fLego(0),
+  fPDGDB(TDatabasePDG::Instance()),        //Particle factory object!
+  fHitLists(new TList()),                  // Create HitLists list
+  fEventEnergy(0),
+  fSummEnergy(0),
+  fSum2Energy(0),
+  fConfigFunction("Config();"),
+  fRandom(new TRandom3()),
+  fMCQA(0),
+  fTransParName("\0"),
+  fBaseFileName("\0"),
+  fStack(new AliStack(10000)),        //Particle stack
+  fTreeDFileName(""),
+  fTreeDFile(0),
+  fTreeSFileName(""),
+  fTreeSFile(0),
+  fTreeRFileName(""),
+  fTreeRFile(0)
 {
   //
   //  Constructor for the main processor.
@@ -275,76 +460,44 @@ AliRun::AliRun(const char *name, const char *title)
   //  Creates the list of Detectors.
   //  Creates the list of particles.
   //
-  Int_t i;
-  
+
   gAlice     = this;
-  fTreeD     = 0;
-  fTreeH     = 0;
-  fTreeE     = 0;
-  fTreeR     = 0;
-  fTreeS     = 0;
-  fTrRmax    = 1.e10;
-  fTrZmax    = 1.e10;
-  fGenerator = 0;
-  fInitDone  = kFALSE;
-  fLego      = 0;
-  fField     = 0;
-  fConfigFunction    = "Config();";
 
   // Set random number generator
-  gRandom = fRandom = new TRandom3();
+  gRandom = fRandom;
 
   if (gSystem->Getenv("CONFIG_SEED")) {
-     gRandom->SetSeed((UInt_t)atoi(gSystem->Getenv("CONFIG_SEED")));
+     gRandom->SetSeed(static_cast<UInt_t>(atoi(gSystem->Getenv("CONFIG_SEED"))));
   }
-  
+
+  // Add to list of browsable  
   gROOT->GetListOfBrowsables()->Add(this,name);
-  //
-  // Particle stack
-  fStack = new AliStack(10000);
-  // create the support list for the various Detectors
-  fModules = new TObjArray(77);
-  //
+
   // Create the TNode geometry for the event display
-  
   BuildSimpleGeometry();
   
-  fHeader    = new AliHeader();
-  fRun       = 0;
-  fEvent     = 0;
-  //
-  fDisplay = 0;
-  //
   // Create default mag field
   SetField();
-  //
-  fMC      = gMC;
-  //
+
   // Prepare the tracking medium lists
-  fImedia = new TArrayI(1000);
-  for(i=0;i<1000;i++) (*fImedia)[i]=-99;
-  //
-  // Make particles
-  fPDGDB     = TDatabasePDG::Instance();        //Particle factory object!
+  for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
 
+  // Add particle list to configuration
   AliConfig::Instance()->Add(fPDGDB); 
-  //
-  // Create HitLists list
-  fHitLists  = new TList();
-  //
+
+  // Set transport parameters
   SetTransPar();
-  fBaseFileName = ".\0";
-  //
-  fDebug        = 0;
 }
 
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 AliRun::~AliRun()
 {
   //
   // Default AliRun destructor
   //
+  TFile *curfil =0;
+  if(fTreeE)curfil=fTreeE->GetCurrentFile();
   delete fImedia;
   delete fField;
   delete fMC;
@@ -354,6 +507,7 @@ AliRun::~AliRun()
   delete fLego;
   delete fTreeD;
   delete fTreeH;
+  delete fTreeTR;
   delete fTreeE;
   delete fTreeR;
   delete fTreeS;
@@ -366,29 +520,59 @@ AliRun::~AliRun()
   delete fPDGDB;
   delete fMCQA;
   delete fHeader;
+  // avoid to delete TFile objects not owned by this object
+  // avoid multiple deletions
+  if(curfil == fTreeDFile) fTreeDFile=0;
+  if(curfil == fTreeSFile) fTreeSFile=0;
+  if(curfil == fTreeRFile) fTreeRFile=0;
+  if(fTreeSFile == fTreeDFile) fTreeSFile=0;
+  if(fTreeRFile == fTreeDFile) fTreeRFile=0;
+  if(fTreeRFile == fTreeSFile) fTreeRFile=0;
+  if(fTreeDFile){
+    if(fTreeDFile->IsOpen())fTreeDFile->Close();
+    delete fTreeDFile;
+  }
+  if(fTreeSFile){
+    if(fTreeSFile->IsOpen())fTreeSFile->Close();
+    delete fTreeSFile;
+  }
+  if(fTreeRFile){
+    if(fTreeRFile->IsOpen())fTreeRFile->Close();
+    delete fTreeRFile;
+  }
+  if (gROOT->GetListOfBrowsables())
+    gROOT->GetListOfBrowsables()->Remove(this);
+
+  gAlice=0;
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
+void AliRun::Copy(AliRun &) const
+{
+  Fatal("Copy","Not implemented!\n");
+}
+
+//_______________________________________________________________________
 void AliRun::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
 {
   //
   //  Add a hit to detector id
   //
   TObjArray &dets = *fModules;
-  if(dets[id]) ((AliModule*) dets[id])->AddHit(track,vol,hits);
+  if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
 {
   //
   // Add digit to detector id
   //
   TObjArray &dets = *fModules;
-  if(dets[id]) ((AliModule*) dets[id])->AddDigit(tracks,digits);
+  if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::Browse(TBrowser *b)
 {
   //
@@ -401,6 +585,7 @@ void AliRun::Browse(TBrowser *b)
     
   if (pTreeK) b->Add(pTreeK,pTreeK->GetName());
   if (fTreeH) b->Add(fTreeH,fTreeH->GetName());
+  if (fTreeTR) b->Add(fTreeTR,fTreeH->GetName());
   if (fTreeD) b->Add(fTreeD,fTreeD->GetName());
   if (fTreeE) b->Add(fTreeE,fTreeE->GetName());
   if (fTreeR) b->Add(fTreeR,fTreeR->GetName());
@@ -408,13 +593,13 @@ void AliRun::Browse(TBrowser *b)
   
   TIter next(fModules);
   AliModule *detector;
-  while((detector = (AliModule*)next())) {
+  while((detector = dynamic_cast<AliModule*>(next()))) {
     b->Add(detector,detector->GetName());
   }
   b->Add(fMCQA,"AliMCQA");
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::Build()
 {
   //
@@ -423,7 +608,7 @@ void AliRun::Build()
   //
 }
  
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::BuildSimpleGeometry()
 {
   //
@@ -438,7 +623,7 @@ void AliRun::BuildSimpleGeometry()
   new TNode("alice","alice","S_alice");
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::CleanDetectors()
 {
   //
@@ -446,12 +631,12 @@ void AliRun::CleanDetectors()
   //
   TIter next(fModules);
   AliModule *detector;
-  while((detector = (AliModule*)next())) {
+  while((detector = dynamic_cast<AliModule*>(next()))) {
     detector->FinishEvent();
   }
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 Int_t AliRun::DistancetoPrimitive(Int_t, Int_t)
 {
   //
@@ -461,7 +646,7 @@ Int_t AliRun::DistancetoPrimitive(Int_t, Int_t)
   return 9999;
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::DumpPart (Int_t i) const
 {
   //
@@ -470,7 +655,7 @@ void AliRun::DumpPart (Int_t i) const
     fStack->DumpPart(i);
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::DumpPStack () const
 {
   //
@@ -479,7 +664,7 @@ void AliRun::DumpPStack () const
     fStack->DumpPStack();
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void  AliRun::SetField(AliMagF* magField)
 {
     // Set Magnetic Field Map
@@ -487,7 +672,7 @@ void  AliRun::SetField(AliMagF* magField)
     fField->ReadField();
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::SetField(Int_t type, Int_t version, Float_t scale,
                      Float_t maxField, char* filename)
 {
@@ -514,113 +699,7 @@ void AliRun::SetField(Int_t type, Int_t version, Float_t scale,
   }
 }
  
-//_____________________________________________________________________________
-void AliRun::PreTrack()
-{
-     TObjArray &dets = *fModules;
-     AliModule *module;
-
-     for(Int_t i=0; i<=fNdets; i++)
-       if((module = (AliModule*)dets[i]))
-        module->PreTrack();
-
-     fMCQA->PreTrack();
-}
-
-//_____________________________________________________________________________
-void AliRun::PostTrack()
-{
-     TObjArray &dets = *fModules;
-     AliModule *module;
-
-     for(Int_t i=0; i<=fNdets; i++)
-       if((module = (AliModule*)dets[i]))
-        module->PostTrack();
-}
-
-//_____________________________________________________________________________
-void AliRun::FinishPrimary()
-{
-  //
-  // Called  at the end of each primary track
-  //
-  
-  //  static Int_t count=0;
-  //  const Int_t times=10;
-  // This primary is finished, purify stack
-  fStack->PurifyKine();
-
-  TIter next(fModules);
-  AliModule *detector;
-  while((detector = (AliModule*)next())) {
-    detector->FinishPrimary();
-  }
-
-  // Write out hits if any
-  if (gAlice->TreeH()) {
-    gAlice->TreeH()->Fill();
-  }
-  
-  //
-  //  if(++count%times==1) gObjectTable->Print();
-}
-
-//_____________________________________________________________________________
-void AliRun::BeginPrimary()
-{
-  //
-  // Called  at the beginning of each primary track
-  //
-  
-  // Reset Hits info
-  gAlice->ResetHits();
-
-}
-
-//_____________________________________________________________________________
-void AliRun::FinishEvent()
-{
-  //
-  // Called at the end of the event.
-  //
-  
-  //
-  if(fLego) fLego->FinishEvent();
-
-  //Update the energy deposit tables
-  Int_t i;
-  for(i=0;i<fEventEnergy.GetSize();i++) {
-    fSummEnergy[i]+=fEventEnergy[i];
-    fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
-  }
-
-
-  
-  // Update Header information 
-
-  fHeader->SetNprimary(fStack->GetNprimary());
-  fHeader->SetNtrack(fStack->GetNtrack());  
-
-  
-  // Write out the kinematics
-  fStack->FinishEvent();
-   
-  // Write out the event Header information
-  if (fTreeE) {
-      fHeader->SetStack(fStack);
-      fTreeE->Fill();
-  }
-  
-  
-  // Write Tree headers
-  TTree* pTreeK = fStack->TreeK();
-  if (pTreeK) pTreeK->Write(0,TObject::kOverwrite);
-  if (fTreeH) fTreeH->Write(0,TObject::kOverwrite);
-  
-  ++fEvent;
-}
-
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::FinishRun()
 {
   //
@@ -633,7 +712,7 @@ void AliRun::FinishRun()
   // Clean detector information
   TIter next(fModules);
   AliModule *detector;
-  while((detector = (AliModule*)next())) {
+  while((detector = dynamic_cast<AliModule*>(next()))) {
     detector->FinishRun();
   }
   
@@ -656,24 +735,28 @@ void AliRun::FinishRun()
   if (fTreeH) {
     delete fTreeH; fTreeH = 0;
   }
+  if (fTreeTR) {
+    delete fTreeTR; fTreeTR = 0;
+  }
   if (fTreeD) {
     delete fTreeD; fTreeD = 0;
   }
   if (fTreeR) {
     delete fTreeR; fTreeR = 0;
   }
-  if (fTreeE) {
-    delete fTreeE; fTreeE = 0;
-  }
+//   if (fTreeE) {
+//     delete fTreeE; fTreeE = 0;
+//   }
   if (fTreeS) {
     delete fTreeS; fTreeS = 0;
   }
-    
+  fGenerator->FinishRun();
+  
   // Close output file
   file->Write();
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::FlagTrack(Int_t track)
 {
   // Delegate to stack
@@ -681,7 +764,7 @@ void AliRun::FlagTrack(Int_t track)
     fStack->FlagTrack(track);
 }
  
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::EnergySummary()
 {
   //
@@ -710,7 +793,7 @@ void AliRun::EnergySummary()
        } else 
          ed2=99;
        fSummEnergy[ndep]=ed;
-       fSum2Energy[ndep]=TMath::Min((Float_t) 99.,TMath::Max(ed2,kzero));
+       fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
        edtot+=ed;
        ndep++;
       }
@@ -747,25 +830,25 @@ void AliRun::EnergySummary()
   //  fSum2Energy.Set(0);
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 AliModule *AliRun::GetModule(const char *name) const
 {
   //
   // Return pointer to detector from name
   //
-  return (AliModule*)fModules->FindObject(name);
+  return dynamic_cast<AliModule*>(fModules->FindObject(name));
 }
  
-//_____________________________________________________________________________
+//_______________________________________________________________________
 AliDetector *AliRun::GetDetector(const char *name) const
 {
   //
   // Return pointer to detector from name
   //
-  return (AliDetector*)fModules->FindObject(name);
+  return dynamic_cast<AliDetector*>(fModules->FindObject(name));
 }
  
-//_____________________________________________________________________________
+//_______________________________________________________________________
 Int_t AliRun::GetModuleID(const char *name) const
 {
   //
@@ -777,7 +860,7 @@ Int_t AliRun::GetModuleID(const char *name) const
   return i;
 }
  
-//_____________________________________________________________________________
+//_______________________________________________________________________
 Int_t AliRun::GetEvent(Int_t event)
 {
   //
@@ -787,14 +870,16 @@ Int_t AliRun::GetEvent(Int_t event)
 
   // Reset existing structures
   ResetHits();
+  ResetTrackReferences();
   ResetDigits();
   ResetSDigits();
   
   // Delete Trees already connected
-  if (fTreeH) delete fTreeH;
-  if (fTreeD) delete fTreeD;
-  if (fTreeR) delete fTreeR;
-  if (fTreeS) delete fTreeS;
+  if (fTreeH) { delete fTreeH; fTreeH = 0;}
+  if (fTreeTR) { delete fTreeTR; fTreeTR = 0;}
+  if (fTreeD) { delete fTreeD; fTreeD = 0;}
+  if (fTreeR) { delete fTreeR; fTreeR = 0;}
+  if (fTreeS) { delete fTreeS; fTreeS = 0;}
 
  // Create the particle stack
   if (fHeader) delete fHeader; 
@@ -803,48 +888,111 @@ Int_t AliRun::GetEvent(Int_t event)
   // Get header from file
   if(fTreeE) {
       fTreeE->SetBranchAddress("Header", &fHeader);
-      fTreeE->GetEntry(event);
+
+      if (!fTreeE->GetEntry(event)) {
+       Error("GetEvent","Cannot find event:%d\n",event);
+       return -1;
+      }
   }  
-  else 
+  else {
       Error("GetEvent","Cannot find Header Tree (TE)\n");
+      return -1;
+  }
 
-  // Get the stack from the header
-  if (fStack) delete fStack;
-  fStack = fHeader->Stack();
-  fStack->GetEvent(event);
+  // Get the stack from the header, set fStack to 0 if it 
+  // fails to get event
   //
   TFile *file = fTreeE->GetCurrentFile();
   char treeName[20];
+  
   file->cd();
 
+  if (fStack) delete fStack;
+  fStack = fHeader->Stack();
+  if (fStack) {
+    if (!fStack->GetEvent(event)) fStack = 0;
+  }
+
   // Get Hits Tree header from file
   sprintf(treeName,"TreeH%d",event);
-  fTreeH = (TTree*)gDirectory->Get(treeName);
+  fTreeH = dynamic_cast<TTree*>(gDirectory->Get(treeName));
   if (!fTreeH) {
-    Error("GetEvent","cannot find Hits Tree for event:%d\n",event);
+      Warning("GetEvent","cannot find Hits Tree for event:%d\n",event);
   }
-  
+
+  // Get TracReferences Tree header from file
+  sprintf(treeName,"TreeTR%d",event);
+  fTreeTR = dynamic_cast<TTree*>(gDirectory->Get(treeName));
+  if (!fTreeTR) {
+    Warning("GetEvent","cannot find TrackReferences Tree for event:%d\n",event);
+  }
+
+  // get current file name and compare with names containing trees S,D,R
+  TString curfilname=static_cast<TString>(fTreeE->GetCurrentFile()->GetName());
+  if(fTreeDFileName==curfilname)fTreeDFileName="";
+  if(fTreeSFileName==curfilname)fTreeSFileName="";
+  if(fTreeRFileName==curfilname)fTreeRFileName="";
+
   // Get Digits Tree header from file
   sprintf(treeName,"TreeD%d",event);
-  fTreeD = (TTree*)gDirectory->Get(treeName);
+
+  if (!fTreeDFile && fTreeDFileName != "") {
+    InitTreeFile("D",fTreeDFileName);
+  }    
+  if (fTreeDFile) {    
+    fTreeD = dynamic_cast<TTree*>(fTreeDFile->Get(treeName));
+  } else {
+    fTreeD = dynamic_cast<TTree*>(file->Get(treeName));
+  }
   if (!fTreeD) {
     // Warning("GetEvent","cannot find Digits Tree for event:%d\n",event);
   }
+  if(fTreeDFileName != ""){
+    if(fTreeDFileName==fTreeSFileName) {
+      fTreeSFileName = "";
+      fTreeSFile = fTreeDFile;
+    }
+    if(fTreeDFileName==fTreeRFileName) {
+      fTreeRFileName = "";
+      fTreeRFile = fTreeDFile;
+    }
+  }
 
   file->cd();
 
   // Get SDigits Tree header from file
   sprintf(treeName,"TreeS%d",event);
-  fTreeS = (TTree*)gDirectory->Get(treeName);
+  if (!fTreeSFile && fTreeSFileName != "") {
+    InitTreeFile("S",fTreeSFileName);
+  } 
+  if (fTreeSFile) {
+    fTreeS = dynamic_cast<TTree*>(fTreeSFile->Get(treeName));
+  } else {
+    fTreeS = dynamic_cast<TTree*>(gDirectory->Get(treeName));
+  }
   if (!fTreeS) {
     // Warning("GetEvent","cannot find SDigits Tree for event:%d\n",event);
   }
 
+  if(fTreeSFileName != ""){
+    if(fTreeSFileName==fTreeRFileName){
+      fTreeRFileName = "";
+      fTreeRFile = fTreeSFile;
+    }
+  }
+
   file->cd();
   
   // Get Reconstruct Tree header from file
   sprintf(treeName,"TreeR%d",event);
-  fTreeR = (TTree*)gDirectory->Get(treeName);
+  if (!fTreeRFile && fTreeRFileName != "") {
+    InitTreeFile("R",fTreeRFileName);
+  } 
+  if(fTreeRFile) {
+    fTreeR = dynamic_cast<TTree*>(fTreeRFile->Get(treeName));
+  } else {
+    fTreeR = dynamic_cast<TTree*>(gDirectory->Get(treeName));
+  }
   if (!fTreeR) {
     //    printf("WARNING: cannot find Reconstructed Tree for event:%d\n",event);
   }
@@ -854,21 +1002,23 @@ Int_t AliRun::GetEvent(Int_t event)
   // Set Trees branch addresses
   TIter next(fModules);
   AliModule *detector;
-  while((detector = (AliModule*)next())) {
+  while((detector = dynamic_cast<AliModule*>(next()))) {
     detector->SetTreeAddress();
   }
 
-  return fStack->GetNtrack();
+  fEvent=event; //MI change
+
+  return fHeader->GetNtrack();
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 TGeometry *AliRun::GetGeometry()
 {
   //
   // Import Alice geometry from current file
   // Return pointer to geometry object
   //
-  if (!fGeometry) fGeometry = (TGeometry*)gDirectory->Get("AliceGeom");
+  if (!fGeometry) fGeometry = dynamic_cast<TGeometry*>(gDirectory->Get("AliceGeom"));
   //
   // Unlink and relink nodes in detectors
   // This is bad and there must be a better way...
@@ -876,12 +1026,12 @@ TGeometry *AliRun::GetGeometry()
   
   TIter next(fModules);
   AliModule *detector;
-  while((detector = (AliModule*)next())) {
+  while((detector = dynamic_cast<AliModule*>(next()))) {
     TList *dnodes=detector->Nodes();
     Int_t j;
     TNode *node, *node1;
     for ( j=0; j<dnodes->GetSize(); j++) {
-      node = (TNode*) dnodes->At(j);
+      node = dynamic_cast<TNode*>(dnodes->At(j));
       node1 = fGeometry->GetNode(node->GetName());
       dnodes->Remove(node);
       dnodes->AddAt(node1,j);
@@ -890,17 +1040,7 @@ TGeometry *AliRun::GetGeometry()
   return fGeometry;
 }
 
-//_____________________________________________________________________________
-void AliRun::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom,
-                         Float_t &e, Float_t *vpos, Float_t *polar,
-                         Float_t &tof)
-{
-  // Delegate to stack
-  //
-    fStack->GetNextTrack(mtrack, ipart, pmom, e, vpos, polar, tof);  
-}
-
-//_____________________________________________________________________________
+//_______________________________________________________________________
 Int_t AliRun::GetPrimary(Int_t track) const
 {
   //
@@ -909,90 +1049,22 @@ Int_t AliRun::GetPrimary(Int_t track) const
     return fStack->GetPrimary(track);
 }
  
-//_____________________________________________________________________________
-void AliRun::InitMC(const char *setup)
+//_______________________________________________________________________
+void AliRun::MediaTable()
 {
   //
-  // Initialize the Alice setup
+  // Built media table to get from the media number to
+  // the detector id
   //
-
-  if(fInitDone) {
-    Warning("Init","Cannot initialise AliRun twice!\n");
-    return;
-  }
-    
-  gROOT->LoadMacro(setup);
-  gInterpreter->ProcessLine(fConfigFunction.Data());
-
-
-  gMC->DefineParticles();  //Create standard MC particles
-
-  TObject *objfirst, *objlast;
-
-  fNdets = fModules->GetLast()+1;
-
-  //
-  //=================Create Materials and geometry
-  gMC->Init();
-
-  // Added also after in case of interactive initialisation of modules
-  fNdets = fModules->GetLast()+1;
-
-   TIter next(fModules);
-   AliModule *detector;
-   while((detector = (AliModule*)next())) {
-      detector->SetTreeAddress();
-      objlast = gDirectory->GetList()->Last();
-      
-      // Add Detector histograms in Detector list of histograms
-      if (objlast) objfirst = gDirectory->GetList()->After(objlast);
-      else         objfirst = gDirectory->GetList()->First();
-      while (objfirst) {
-       detector->Histograms()->Add(objfirst);
-       objfirst = gDirectory->GetList()->After(objfirst);
-      }
-   }
-   ReadTransPar(); //Read the cuts for all materials
-   
-   MediaTable(); //Build the special IMEDIA table
-   
-   //Initialise geometry deposition table
-   fEventEnergy.Set(gMC->NofVolumes()+1);
-   fSummEnergy.Set(gMC->NofVolumes()+1);
-   fSum2Energy.Set(gMC->NofVolumes()+1);
-   
-   //Compute cross-sections
-   gMC->BuildPhysics();
-   
-   //Write Geometry object to current file.
-   fGeometry->Write();
-   
-   fInitDone = kTRUE;
-
-   fMCQA = new AliMCQA(fNdets);
-
-   AliConfig::Instance();
-   //
-   // Save stuff at the beginning of the file to avoid file corruption
-   Write();
-}
-
-//_____________________________________________________________________________
-void AliRun::MediaTable()
-{
-  //
-  // Built media table to get from the media number to
-  // the detector id
-  //
-  Int_t kz, nz, idt, lz, i, k, ind;
-  //  Int_t ibeg;
-  TObjArray &dets = *gAlice->Detectors();
-  AliModule *det;
+  Int_t kz, nz, idt, lz, i, k, ind;
+  //  Int_t ibeg;
+  TObjArray &dets = *gAlice->Detectors();
+  AliModule *det;
   //
   // For all detectors
   for (kz=0;kz<fNdets;kz++) {
     // If detector is defined
-    if((det=(AliModule*) dets[kz])) {
+    if((det=dynamic_cast<AliModule*>(dets[kz]))) {
         TArrayI &idtmed = *(det->GetIdtmed()); 
         for(nz=0;nz<100;nz++) {
        // Find max and min material number
@@ -1023,7 +1095,7 @@ void AliRun::MediaTable()
   for(i=0;i<(fNdets-1)/6+1;i++) {
     for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
       ind=i*6+k;
-      det=(AliModule*)dets[ind];
+      det=dynamic_cast<AliModule*>(dets[ind]);
       if(det)
        printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
               det->HiMedium());
@@ -1034,7 +1106,7 @@ void AliRun::MediaTable()
   }
 }
 
-//____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::SetGenerator(AliGenerator *generator)
 {
   //
@@ -1043,7 +1115,7 @@ void AliRun::SetGenerator(AliGenerator *generator)
   if(!fGenerator) fGenerator = generator;
 }
 
-//____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::ResetGenerator(AliGenerator *generator)
 {
   //
@@ -1059,19 +1131,19 @@ void AliRun::ResetGenerator(AliGenerator *generator)
   fGenerator = generator;
 }
 
-//____________________________________________________________________________
-void AliRun::SetTransPar(char *filename)
+//_______________________________________________________________________
+void AliRun::SetTransPar(const char *filename)
 {
   fTransParName = filename;
 }
 
-//____________________________________________________________________________
-void AliRun::SetBaseFile(char *filename)
+//_______________________________________________________________________
+void AliRun::SetBaseFile(const char *filename)
 {
   fBaseFileName = filename;
 }
 
-//____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::ReadTransPar()
 {
   //
@@ -1182,7 +1254,7 @@ void AliRun::ReadTransPar()
 }
 
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::MakeTree(Option_t *option, const char *file)
 {
   //
@@ -1195,6 +1267,7 @@ void AliRun::MakeTree(Option_t *option, const char *file)
   // Analyse options
   const char *oK = strstr(option,"K");
   const char *oH = strstr(option,"H");
+  const char *oTR = strstr(option,"T");
   const char *oE = strstr(option,"E");
   const char *oD = strstr(option,"D");
   const char *oR = strstr(option,"R");
@@ -1209,9 +1282,10 @@ void AliRun::MakeTree(Option_t *option, const char *file)
 
   if (oE && !fTreeE) {
     fTreeE = new TTree("TE","Header");
+    //    branch = fTreeE->Branch("Header", "AliHeader", &fHeader, 4000, 0);         
     branch = fTreeE->Branch("Header", "AliHeader", &fHeader, 4000, 0);         
     branch->SetAutoDelete(kFALSE);          
-    TFolder *folder = (TFolder *)gROOT->FindObjectAny("/Folders/RunMC/Event/Header");
+    TFolder *folder = dynamic_cast<TFolder *>(gROOT->FindObjectAny("/Folders/RunMC/Event/Header"));
     if (folder) folder->Add(fHeader);
 //    branch = fTreeE->Branch("Stack","AliStack", &fStack, 4000, 0);
 //    branch->SetAutoDelete(kFALSE);        
@@ -1224,7 +1298,7 @@ void AliRun::MakeTree(Option_t *option, const char *file)
     sprintf(outFile,"%s/%s",GetBaseFile(),file);
     branch->SetFile(outFile);
     TIter next( branch->GetListOfBranches());
-    while ((branch=(TBranch*)next())) {
+    while ((branch=dynamic_cast<TBranch*>(next()))) {
        branch->SetFile(outFile);
     } 
     if (GetDebug()>1)
@@ -1239,6 +1313,14 @@ void AliRun::MakeTree(Option_t *option, const char *file)
     fTreeH->SetAutoSave(1000000000); //no autosave
     fTreeH->Write(0,TObject::kOverwrite);
   }
+
+  if (oTR && !fTreeTR) {
+    sprintf(hname,"TreeTR%d",fEvent);
+    fTreeTR = new TTree(hname,"TrackReferences");
+    fTreeTR->SetAutoSave(1000000000); //no autosave
+    fTreeTR->Write(0,TObject::kOverwrite);
+  }
+
   if (oD && !fTreeD) {
     sprintf(hname,"TreeD%d",fEvent);
     fTreeD = new TTree(hname,"Digits");
@@ -1261,77 +1343,19 @@ void AliRun::MakeTree(Option_t *option, const char *file)
   // will be in turn a subbranch of the detector master branch
   TIter next(fModules);
   AliModule *detector;
-  while((detector = (AliModule*)next())) {
+  while((detector = dynamic_cast<AliModule*>(next()))) {
      if (oH) detector->MakeBranch(option,file);
+     if (oTR) detector->MakeBranchTR(option,file);
   }
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 TParticle* AliRun::Particle(Int_t i)
 {
     return fStack->Particle(i);
 }
 
-//_____________________________________________________________________________
-void AliRun::BeginEvent()
-{
-    // Clean-up previous event
-    // Energy scores  
-    fEventEnergy.Reset();  
-    // Clean detector information
-    CleanDetectors();
-    // Reset stack info
-    fStack->Reset();
-
-  //
-  //  Reset all Detectors & kinematics & trees
-  //
-  char hname[30];
-  //
-  // Initialise event header
-  fHeader->Reset(fRun,fEvent);
-  //
-  fStack->BeginEvent(fEvent);
-
-  //
-  if(fLego) {
-    fLego->BeginEvent();
-    return;
-  }
-
-  //
-
-  ResetHits();
-  ResetDigits();
-  ResetSDigits();
-
-
-  if(fTreeH) {
-    fTreeH->Reset();
-    sprintf(hname,"TreeH%d",fEvent);
-    fTreeH->SetName(hname);
-  }
-  if(fTreeD) {
-    fTreeD->Reset();
-    sprintf(hname,"TreeD%d",fEvent);
-    fTreeD->SetName(hname);
-    fTreeD->Write(0,TObject::kOverwrite);
-  }
-  if(fTreeS) {
-    fTreeS->Reset();
-    sprintf(hname,"TreeS%d",fEvent);
-    fTreeS->SetName(hname);
-    fTreeS->Write(0,TObject::kOverwrite);
-  }
-  if(fTreeR) {
-    fTreeR->Reset();
-    sprintf(hname,"TreeR%d",fEvent);
-    fTreeR->SetName(hname);
-    fTreeR->Write(0,TObject::kOverwrite);
-  }
-}
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::ResetDigits()
 {
   //
@@ -1339,12 +1363,12 @@ void AliRun::ResetDigits()
   //
   TIter next(fModules);
   AliModule *detector;
-  while((detector = (AliModule*)next())) {
+  while((detector = dynamic_cast<AliModule*>(next()))) {
      detector->ResetDigits();
   }
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::ResetSDigits()
 {
   //
@@ -1352,12 +1376,12 @@ void AliRun::ResetSDigits()
   //
   TIter next(fModules);
   AliModule *detector;
-  while((detector = (AliModule*)next())) {
+  while((detector = dynamic_cast<AliModule*>(next()))) {
      detector->ResetSDigits();
   }
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::ResetHits()
 {
   //
@@ -1365,12 +1389,25 @@ void AliRun::ResetHits()
   //
   TIter next(fModules);
   AliModule *detector;
-  while((detector = (AliModule*)next())) {
+  while((detector = dynamic_cast<AliModule*>(next()))) {
      detector->ResetHits();
   }
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
+void AliRun::ResetTrackReferences()
+{
+  //
+  //  Reset all Detectors hits
+  //
+  TIter next(fModules);
+  AliModule *detector;
+  while((detector = dynamic_cast<AliModule*>(next()))) {
+     detector->ResetTrackReferences();
+  }
+}
+
+//_______________________________________________________________________
 void AliRun::ResetPoints()
 {
   //
@@ -1378,12 +1415,84 @@ void AliRun::ResetPoints()
   //
   TIter next(fModules);
   AliModule *detector;
-  while((detector = (AliModule*)next())) {
+  while((detector = dynamic_cast<AliModule*>(next()))) {
      detector->ResetPoints();
   }
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
+void AliRun::InitMC(const char *setup)
+{
+  //
+  // Initialize the Alice setup
+  //
+
+  if(fInitDone) {
+    Warning("Init","Cannot initialise AliRun twice!\n");
+    return;
+  }
+    
+  gROOT->LoadMacro(setup);
+  gInterpreter->ProcessLine(fConfigFunction.Data());
+
+  // Register MC in configuration 
+  AliConfig::Instance()->Add(gMC);
+  gMC->SetStack(fStack);
+
+  gMC->DefineParticles();  //Create standard MC particles
+  AliPDG::AddParticlesToPdgDataBase();  
+
+  TObject *objfirst, *objlast;
+
+  fNdets = fModules->GetLast()+1;
+
+  //
+  //=================Create Materials and geometry
+  gMC->Init();
+
+  // Added also after in case of interactive initialisation of modules
+  fNdets = fModules->GetLast()+1;
+
+   TIter next(fModules);
+   AliModule *detector;
+   while((detector = dynamic_cast<AliModule*>(next()))) {
+      detector->SetTreeAddress();
+      objlast = gDirectory->GetList()->Last();
+      
+      // Add Detector histograms in Detector list of histograms
+      if (objlast) objfirst = gDirectory->GetList()->After(objlast);
+      else         objfirst = gDirectory->GetList()->First();
+      while (objfirst) {
+        detector->Histograms()->Add(objfirst);
+        objfirst = gDirectory->GetList()->After(objfirst);
+      }
+   }
+   ReadTransPar(); //Read the cuts for all materials
+   
+   MediaTable(); //Build the special IMEDIA table
+   
+   //Initialise geometry deposition table
+   fEventEnergy.Set(gMC->NofVolumes()+1);
+   fSummEnergy.Set(gMC->NofVolumes()+1);
+   fSum2Energy.Set(gMC->NofVolumes()+1);
+   
+   //Compute cross-sections
+   gMC->BuildPhysics();
+   
+   //Write Geometry object to current file.
+   fGeometry->Write();
+   
+   fInitDone = kTRUE;
+
+   fMCQA = new AliMCQA(fNdets);
+
+   AliConfig::Instance();
+   //
+   // Save stuff at the beginning of the file to avoid file corruption
+   Write();
+}
+
+//_______________________________________________________________________
 void AliRun::RunMC(Int_t nevent, const char *setup)
 {
   //
@@ -1393,13 +1502,13 @@ void AliRun::RunMC(Int_t nevent, const char *setup)
   // a positive number of events will cause the finish routine
   // to be called
   //
-
+    fEventsPerRun = nevent;
   // check if initialisation has been done
   if (!fInitDone) InitMC(setup);
   
   // Create the Root Tree with one branch per detector
 
-   MakeTree("ESDR");
+   MakeTree("ESDRT");
 
   if (gSystem->Getenv("CONFIG_SPLIT_FILE")) {
      MakeTree("K","Kine.root");
@@ -1414,14 +1523,17 @@ void AliRun::RunMC(Int_t nevent, const char *setup)
   if(nevent>0) FinishRun();
 }
 
-//_____________________________________________________________________________
-void AliRun::RunReco(const char *selected)
+//_______________________________________________________________________
+void AliRun::RunReco(const char *selected, Int_t first, Int_t last)
 {
   //
   // Main function to be called to reconstruct Alice event
   // 
-   cout << "Found "<< gAlice->TreeE()->GetEntries() << "events" << endl; 
-   for (Int_t nevent=0; nevent<gAlice->TreeE()->GetEntries(); nevent++) {
+   cout << "Found "<< gAlice->TreeE()->GetEntries() << "events" << endl;
+   Int_t nFirst = first;
+   Int_t nLast  = (last < 0)? static_cast<Int_t>(gAlice->TreeE()->GetEntries()) : last;
+   
+   for (Int_t nevent = nFirst; nevent <= nLast; nevent++) {
      cout << "Processing event "<< nevent << endl;
      GetEvent(nevent);
      // MakeTree("R");
@@ -1429,7 +1541,7 @@ void AliRun::RunReco(const char *selected)
    }
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 
 void AliRun::Hits2Digits(const char *selected)
 {
@@ -1445,7 +1557,7 @@ void AliRun::Hits2Digits(const char *selected)
 }
 
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 
 void AliRun::Tree2Tree(Option_t *option, const char *selected)
 {
@@ -1475,7 +1587,7 @@ void AliRun::Tree2Tree(Option_t *option, const char *selected)
 
    char outFile[32];
    
-   while((detector = (AliDetector*)next())) {
+   while((detector = dynamic_cast<AliDetector*>(next()))) {
      if (selected) 
        if (strcmp(detector->GetName(),selected)) continue;
      if (detector->IsActive()){ 
@@ -1517,7 +1629,7 @@ void AliRun::Tree2Tree(Option_t *option, const char *selected)
 }
 
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min,
                     Float_t c1max,Int_t nc2,Float_t c2min,Float_t c2max,
                     Float_t rmin,Float_t rmax,Float_t zmax, AliLegoGenerator* gener)
@@ -1584,7 +1696,8 @@ void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min,
   
   //Run Lego Object
 
-  gMC->ProcessRun(nc1*nc2+1);
+  //gMC->ProcessRun(nc1*nc2+1);
+  gMC->ProcessRun(nc1*nc2);
   
   // Create only the Root event Tree
   MakeTree("E");
@@ -1597,7 +1710,7 @@ void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min,
   delete fLego; fLego=0;
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::SetConfigFunction(const char * config) 
 {
   //
@@ -1607,7 +1720,7 @@ void AliRun::SetConfigFunction(const char * config)
   fConfigFunction=config;
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::SetCurrentTrack(Int_t track)
 { 
   //
@@ -1616,32 +1729,33 @@ void AliRun::SetCurrentTrack(Int_t track)
     fStack->SetCurrentTrack(track); 
 }
  
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
                      Float_t *vpos, Float_t *polar, Float_t tof,
-                     AliMCProcess mech, Int_t &ntr, Float_t weight)
+                     AliMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
 { 
 // Delegate to stack
 //
-    fStack->SetTrack(done, parent, pdg, pmom, vpos, polar, tof,
-                    mech, ntr, weight);
+
+     fStack->SetTrack(done, parent, pdg, pmom, vpos, polar, tof,
+                    mech, ntr, weight, is);
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg,
                      Double_t px, Double_t py, Double_t pz, Double_t e,
                      Double_t vx, Double_t vy, Double_t vz, Double_t tof,
                      Double_t polx, Double_t poly, Double_t polz,
-                     AliMCProcess mech, Int_t &ntr, Float_t weight)
+                     AliMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
 { 
   // Delegate to stack
   //
     fStack->SetTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
-                  polx, poly, polz, mech, ntr, weight);
+                  polx, poly, polz, mech, ntr, weight, is);
     
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::SetHighWaterMark(const Int_t nt)
 {
     //
@@ -1649,7 +1763,7 @@ void AliRun::SetHighWaterMark(const Int_t nt)
     fStack->SetHighWaterMark(nt);
 }
 
-//_____________________________________________________________________________
+//_______________________________________________________________________
 void AliRun::KeepTrack(const Int_t track)
 { 
   //
@@ -1658,13 +1772,167 @@ void AliRun::KeepTrack(const Int_t track)
     fStack->KeepTrack(track);
 }
  
-//_____________________________________________________________________________
-void AliRun::StepManager(Int_t id) 
+// 
+// MC Application
+// 
+
+//_______________________________________________________________________
+void  AliRun::ConstructGeometry() 
+{
+  //
+  // Create modules, materials, geometry
+  //
+
+    TStopwatch stw;
+    TIter next(fModules);
+    AliModule *detector;
+    printf("Geometry creation:\n");
+    while((detector = dynamic_cast<AliModule*>(next()))) {
+      stw.Start();
+      // Initialise detector materials and geometry
+      detector->CreateMaterials();
+      detector->CreateGeometry();
+      printf("%10s R:%.2fs C:%.2fs\n",
+            detector->GetName(),stw.RealTime(),stw.CpuTime());
+    }
+}
+
+//_______________________________________________________________________
+void  AliRun::InitGeometry()
+{ 
+  //
+  // Initialize detectors and display geometry
+  //
+
+   printf("Initialisation:\n");
+    TStopwatch stw;
+    TIter next(fModules);
+    AliModule *detector;
+    while((detector = dynamic_cast<AliModule*>(next()))) {
+      stw.Start();
+      // Initialise detector and display geometry
+      detector->Init();
+      detector->BuildGeometry();
+      printf("%10s R:%.2fs C:%.2fs\n",
+            detector->GetName(),stw.RealTime(),stw.CpuTime());
+    }
+}
+
+//_______________________________________________________________________
+void  AliRun::GeneratePrimaries()
+{ 
+  //
+  // Generate primary particles and fill them in the stack.
+  //
+
+  Generator()->Generate();
+}
+
+//_______________________________________________________________________
+void AliRun::BeginEvent()
+{
+    // Clean-up previous event
+    // Energy scores  
+    fEventEnergy.Reset();  
+    // Clean detector information
+    CleanDetectors();
+    // Reset stack info
+    fStack->Reset();
+
+  //
+  //  Reset all Detectors & kinematics & trees
+  //
+  char hname[30];
+  //
+  // Initialise event header
+  fHeader->Reset(fRun,fEvent,fEventNrInRun);
+  //
+  fStack->BeginEvent(fEvent);
+
+  //
+  if(fLego) {
+    fLego->BeginEvent();
+    return;
+  }
+
+  //
+
+  ResetHits();
+  ResetTrackReferences();
+  ResetDigits();
+  ResetSDigits();
+
+
+  if(fTreeH) {
+    fTreeH->Reset();
+    sprintf(hname,"TreeH%d",fEvent);
+    fTreeH->SetName(hname);
+  }
+
+  if(fTreeTR) {
+    fTreeTR->Reset();
+    sprintf(hname,"TreeTR%d",fEvent);
+    fTreeTR->SetName(hname);
+  }
+
+  if(fTreeD) {
+    fTreeD->Reset();
+    sprintf(hname,"TreeD%d",fEvent);
+    fTreeD->SetName(hname);
+    fTreeD->Write(0,TObject::kOverwrite);
+  }
+  if(fTreeS) {
+    fTreeS->Reset();
+    sprintf(hname,"TreeS%d",fEvent);
+    fTreeS->SetName(hname);
+    fTreeS->Write(0,TObject::kOverwrite);
+  }
+  if(fTreeR) {
+    fTreeR->Reset();
+    sprintf(hname,"TreeR%d",fEvent);
+    fTreeR->SetName(hname);
+    fTreeR->Write(0,TObject::kOverwrite);
+  }
+}
+
+//_______________________________________________________________________
+void AliRun::BeginPrimary()
+{
+  //
+  // Called  at the beginning of each primary track
+  //
+  
+  // Reset Hits info
+  gAlice->ResetHits();
+  gAlice->ResetTrackReferences();
+
+}
+
+//_______________________________________________________________________
+void AliRun::PreTrack()
+{
+     TObjArray &dets = *fModules;
+     AliModule *module;
+
+     for(Int_t i=0; i<=fNdets; i++)
+       if((module = dynamic_cast<AliModule*>(dets[i])))
+        module->PreTrack();
+
+     fMCQA->PreTrack();
+}
+
+//_______________________________________________________________________
+void AliRun::Stepping() 
 {
   //
   // Called at every step during transport
   //
 
+  Int_t id = DetFromMate(gMC->GetMedium());
+  if (id < 0) return;
+
   //
   // --- If lego option, do it and leave 
   if (fLego)
@@ -1675,42 +1943,151 @@ void AliRun::StepManager(Int_t id)
     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
   
     //Call the appropriate stepping routine;
-    AliModule *det = (AliModule*)fModules->At(id);
-    if(det) {
+    AliModule *det = dynamic_cast<AliModule*>(fModules->At(id));
+    if(det && det->StepManagerIsEnabled()) {
       fMCQA->StepManager(id);
       det->StepManager();
     }
   }
 }
 
-//_____________________________________________________________________________
-void AliRun::Streamer(TBuffer &R__b)
+//_______________________________________________________________________
+void AliRun::PostTrack()
 {
-   // Stream an object of class AliRun.
+     TObjArray &dets = *fModules;
+     AliModule *module;
 
-   if (R__b.IsReading()) {
-      if (!gAlice) gAlice = this;
+     for(Int_t i=0; i<=fNdets; i++)
+       if((module = dynamic_cast<AliModule*>(dets[i])))
+        module->PostTrack();
+}
 
-      AliRun::Class()->ReadBuffer(R__b, this);
-      //
-      gROOT->GetListOfBrowsables()->Add(this,"Run");
+//_______________________________________________________________________
+void AliRun::FinishPrimary()
+{
+  //
+  // Called  at the end of each primary track
+  //
+  
+  //  static Int_t count=0;
+  //  const Int_t times=10;
+  // This primary is finished, purify stack
+  fStack->PurifyKine();
 
-      fTreeE = (TTree*)gDirectory->Get("TE");
-      if (fTreeE) {
-         fTreeE->SetBranchAddress("Header", &fHeader);
-      }
-      
-      else    Error("Streamer","cannot find Header Tree\n");
-      fTreeE->GetEntry(0);
+  TIter next(fModules);
+  AliModule *detector;
+  while((detector = dynamic_cast<AliModule*>(next()))) {
+    detector->FinishPrimary();
+  }
+
+  // Write out hits if any
+  if (gAlice->TreeH()) {
+    gAlice->TreeH()->Fill();
+  }
+
+  // Write out hits if any
+  if (gAlice->TreeTR()) {
+    gAlice->TreeTR()->Fill();
+  }
+  
+  //
+  //  if(++count%times==1) gObjectTable->Print();
+}
+
+//_______________________________________________________________________
+void AliRun::FinishEvent()
+{
+  //
+  // Called at the end of the event.
+  //
+  
+  //
+  if(fLego) fLego->FinishEvent();
+
+  //Update the energy deposit tables
+  Int_t i;
+  for(i=0;i<fEventEnergy.GetSize();i++) {
+    fSummEnergy[i]+=fEventEnergy[i];
+    fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
+  }
+
+
+  
+  // Update Header information 
+
+  fHeader->SetNprimary(fStack->GetNprimary());
+  fHeader->SetNtrack(fStack->GetNtrack());  
+
+  
+  // Write out the kinematics
+  fStack->FinishEvent();
+   
+  // Write out the event Header information
+  if (fTreeE) {
+      fHeader->SetStack(fStack);
+      fTreeE->Fill();
+  }
+  
+  
+  // Write Tree headers
+  TTree* pTreeK = fStack->TreeK();
+  if (pTreeK) pTreeK->Write(0,TObject::kOverwrite);
+  if (fTreeH) fTreeH->Write(0,TObject::kOverwrite);
+  if (fTreeTR) fTreeTR->Write(0,TObject::kOverwrite);
+  
+  ++fEvent;
+  ++fEventNrInRun;
+}
 
-      gRandom = fRandom;
-   } else {
-      AliRun::Class()->WriteBuffer(R__b, this);
+//_______________________________________________________________________
+void AliRun::Field(const Double_t* x, Double_t *b) const
+{
+   Float_t xfloat[3];
+   for (Int_t i=0; i<3; i++) xfloat[i] = x[i]; 
+
+   if (Field()) {
+         Float_t bfloat[3];
+         Field()->Field(xfloat,bfloat);
+         for (Int_t j=0; j<3; j++) b[j] = bfloat[j]; 
+   } 
+   else {
+         printf("No mag field defined!\n");
+         b[0]=b[1]=b[2]=0.;
    }
+
+}      
+
+// 
+// End of MC Application
+// 
+
+//_______________________________________________________________________
+void AliRun::Streamer(TBuffer &R__b)
+{
+  // Stream an object of class AliRun.
+
+  if (R__b.IsReading()) {
+    if (!gAlice) gAlice = this;
+
+    AliRun::Class()->ReadBuffer(R__b, this);
+    //
+    gROOT->GetListOfBrowsables()->Add(this,"Run");
+
+    fTreeE = dynamic_cast<TTree*>(gDirectory->Get("TE"));
+    if (fTreeE) {
+         fTreeE->SetBranchAddress("Header", &fHeader);
+    }      
+    else    Error("Streamer","cannot find Header Tree\n");
+    
+    fTreeE->GetEntry(0);
+    gRandom = fRandom;
+  } else {
+    AliRun::Class()->WriteBuffer(R__b, this);
+  }
 }
 
 
-//___________________________________________________________________________
+//_______________________________________________________________________
 Int_t AliRun::CurrentTrack() const {
   //
   // Returns current track
@@ -1718,7 +2095,7 @@ Int_t AliRun::CurrentTrack() const {
   return fStack->CurrentTrack();
 }
 
-//___________________________________________________________________________
+//_______________________________________________________________________
 Int_t AliRun::GetNtrack() const {
   //
   // Returns number of tracks in stack
@@ -1726,7 +2103,7 @@ Int_t AliRun::GetNtrack() const {
   return fStack->GetNtrack();
 }
 
-//___________________________________________________________________________
+//_______________________________________________________________________
 TObjArray* AliRun::Particles() {
   //
   // Returns pointer to Particles array
@@ -1734,10 +2111,272 @@ TObjArray* AliRun::Particles() {
   return fStack->Particles();
 }
 
-//___________________________________________________________________________
+//_______________________________________________________________________
 TTree* AliRun::TreeK() {
   //
   // Returns pointer to the TreeK array
   //
   return fStack->TreeK();
 }
+
+
+//_______________________________________________________________________
+void AliRun::SetGenEventHeader(AliGenEventHeader* header)
+{
+    fHeader->SetGenEventHeader(header);
+}
+
+//_______________________________________________________________________
+TFile* AliRun::InitFile(TString fileName)
+{
+// 
+// create the file where the whole tree will be saved
+//
+  TDirectory *wd = gDirectory;
+  TFile* file = TFile::Open(fileName,"update");
+  gDirectory = wd;
+  if (!file->IsOpen()) {
+    Error("Cannot open file, %s\n",fileName);
+    return 0;
+  }
+  return file;
+}
+
+//_______________________________________________________________________
+TFile* AliRun::InitTreeFile(Option_t *option, TString fileName)
+{
+  //
+  // create the file where one of the following trees will be saved
+  // trees:   S,D,R
+  // WARNING: by default these trees are saved on the file on which
+  // hits are stored. If you divert one of these trees, you cannot restore
+  // it to the original file (usually galice.root) in the same aliroot session
+  Bool_t oS = (strstr(option,"S")!=0);
+  Bool_t oR = (strstr(option,"R")!=0);
+  Bool_t oD = (strstr(option,"D")!=0);
+  Int_t choice[3]; 
+  for (Int_t i=0; i<3; i++) choice[i] = 0;
+  if(oS)choice[0] = 1; 
+  if(oD)choice[1] = 1;
+  if(oR)choice[2] = 1;
+
+  TFile *ptr=0;
+
+  if(!(oS || oR || oD))return ptr;
+
+  Int_t active[3];
+  for (Int_t i=0; i<3; i++) active[i] = 0;
+  if(fTreeSFileName != "") active[0] = 1;
+  if(fTreeDFileName != "") active[1] = 1;
+  if(fTreeDFileName != "") active[2] = 1;
+
+  Bool_t alreadyopen1 = kFALSE;
+  Bool_t alreadyopen2 = kFALSE;
+
+  if(oS){
+    // if already active and same name with non-null ptr
+    if(active[0]==1 && fileName == fTreeSFileName && fTreeSFile){
+      Warning("InitTreeFile","File %s already opened",fTreeSFileName.Data());
+      ptr = fTreeSFile;
+    }
+    else {
+      // if already active with different name with non-null ptr
+      if(active[0]==1 && fileName != fTreeSFileName && fTreeSFile){
+        // close the active files and also the other possible files in option
+        CloseTreeFile(option);
+      }
+      fTreeSFileName = fileName;
+      alreadyopen1 = 
+        (active[1] == 1 && fTreeDFileName == fTreeSFileName && fTreeDFile);
+      alreadyopen2 =
+        (active[2] == 1 && fTreeRFileName == fTreeSFileName && fTreeRFile);
+      if(!(alreadyopen1 || alreadyopen2)){
+        ptr = InitFile(fileName);
+        fTreeSFile = ptr;
+      }
+      else {
+        if(alreadyopen1){fTreeSFile = fTreeDFile; ptr = fTreeSFile;}
+        if(alreadyopen2){fTreeSFile = fTreeRFile; ptr = fTreeSFile;}
+      }
+      if(choice[1] == 1) { fTreeDFileName = fileName; fTreeDFile = ptr;}
+      if(choice[2] == 1) { fTreeRFileName = fileName; fTreeRFile = ptr;}
+    }
+    return ptr;
+  }
+
+  if(oD){
+    // if already active and same name with non-null ptr
+    if(active[1]==1 && fileName == fTreeDFileName && fTreeDFile){
+      Warning("InitTreeFile","File %s already opened",fTreeDFileName.Data());
+      ptr = fTreeDFile;
+    }
+    else {
+      // if already active with different name with non-null ptr
+      if(active[1]==1 && fileName != fTreeDFileName && fTreeDFile){
+        // close the active files and also the other possible files in option
+        CloseTreeFile(option);
+      }
+      fTreeDFileName = fileName;
+      alreadyopen1 = 
+        (active[0] == 1 && fTreeSFileName == fTreeDFileName && fTreeSFile);
+      alreadyopen2 =
+        (active[2] == 1 && fTreeRFileName == fTreeDFileName && fTreeRFile);
+      if(!(alreadyopen1 || alreadyopen2)){
+        ptr = InitFile(fileName);
+        fTreeDFile = ptr;
+      }
+      else {
+        if(alreadyopen1){fTreeDFile = fTreeSFile; ptr = fTreeDFile;}
+        if(alreadyopen2){fTreeDFile = fTreeRFile; ptr = fTreeDFile;}
+      }
+      if(choice[2] == 1) { fTreeRFileName = fileName; fTreeRFile = ptr;}
+    }
+    return ptr;
+  }
+
+  if(oR){
+    // if already active and same name with non-null ptr
+    if(active[2]==1 && fileName == fTreeRFileName && fTreeRFile){
+      Warning("InitTreeFile","File %s already opened",fTreeRFileName.Data());
+      ptr = fTreeRFile;
+    }
+    else {
+      // if already active with different name with non-null ptr
+      if(active[2]==1 && fileName != fTreeRFileName && fTreeRFile){
+        // close the active files and also the other possible files in option
+        CloseTreeFile(option);
+      }
+      fTreeRFileName = fileName;
+      alreadyopen1 = 
+        (active[1] == 1 && fTreeDFileName == fTreeRFileName && fTreeDFile);
+      alreadyopen2 =
+        (active[0]== 1 && fTreeSFileName == fTreeRFileName && fTreeSFile);
+      if(!(alreadyopen1 || alreadyopen2)){
+        ptr = InitFile(fileName);
+        fTreeRFile = ptr;
+      }
+      else {
+        if(alreadyopen1){fTreeRFile = fTreeDFile; ptr = fTreeRFile;}
+        if(alreadyopen2){fTreeRFile = fTreeSFile; ptr = fTreeRFile;}
+      }
+    }
+    return ptr;
+  }
+  return 0;
+}
+
+//_______________________________________________________________________
+void AliRun::PrintTreeFile()
+{
+  //
+  // prints the file names and pointer associated to S,D,R trees
+  //
+  cout<<"===================================================\n";
+  TFile *file = fTreeE->GetCurrentFile();
+  TString curfilname="";
+  if(file)curfilname=static_cast<TString>(file->GetName());
+  cout<<" Current tree file name: "<<curfilname<<endl;
+  cout<<"Pointer: "<<file<<endl;
+  cout<<" Tree S File name: "<<fTreeSFileName<<endl;
+  cout<<"Pointer: "<<fTreeSFile<<endl<<endl;
+  cout<<" Tree D File name: "<<fTreeDFileName<<endl;
+  cout<<"Pointer: "<<fTreeDFile<<endl<<endl;
+  cout<<" Tree R File name: "<<fTreeRFileName<<endl;
+  cout<<"Pointer: "<<fTreeRFile<<endl<<endl;
+  cout<<"===================================================\n";
+}
+//_______________________________________________________________________
+void AliRun::CloseTreeFile(Option_t *option)
+{
+  // 
+  // closes the file containing the tree specified in option
+  // (S,D,R)
+  //
+  Bool_t oS = (strstr(option,"S")!=0);
+  Bool_t oR = (strstr(option,"R")!=0);
+  Bool_t oD = (strstr(option,"D")!=0);
+  Bool_t none = !(oS || oR || oD);
+  if(none)return;
+  if(oS){
+    fTreeSFileName = "";
+    if(fTreeSFile){
+      if(!((fTreeSFile == fTreeDFile) || (fTreeSFile == fTreeRFile)) &&
+           fTreeSFile->IsOpen()){
+        fTreeSFile->Close();
+        delete fTreeSFile;
+      }
+    }
+    fTreeSFile = 0;
+  }
+  if(oD){
+    fTreeDFileName = "";
+    if(fTreeDFile){
+      if(!((fTreeDFile == fTreeRFile) || (fTreeDFile == fTreeSFile)) &&
+         fTreeDFile->IsOpen()){
+        fTreeDFile->Close();
+        delete fTreeDFile;
+      }
+    }
+    fTreeDFile = 0;
+  }
+  if(oR){
+    fTreeRFileName = "";
+    if(fTreeRFile){
+      if(!((fTreeRFile == fTreeSFile) || (fTreeRFile == fTreeDFile)) &&
+         fTreeRFile->IsOpen()){
+        fTreeRFile->Close();
+        delete fTreeRFile;
+      }
+    }
+    fTreeRFile = 0;
+  }
+}
+
+//_______________________________________________________________________
+void AliRun::MakeTree(Option_t *option, TFile *file)
+{
+  //
+  //  Create some trees in the separate file
+  //
+  const char *oD = strstr(option,"D");
+  const char *oR = strstr(option,"R");
+  const char *oS = strstr(option,"S");
+
+  TDirectory *cwd = gDirectory;
+  char hname[30];
+
+  if (oD) {
+    delete fTreeD;
+    sprintf(hname,"TreeD%d",fEvent);
+    file->cd();
+    fTreeD = static_cast<TTree*>(file->Get("hname"));
+    if (!fTreeD) {
+      fTreeD = new TTree(hname,"Digits");
+      fTreeD->Write(0,TObject::kOverwrite);
+    }
+    cwd->cd();
+  }
+  if (oS) {
+    delete fTreeS;
+    sprintf(hname,"TreeS%d",fEvent);
+    file->cd();
+    fTreeS = static_cast<TTree*>(file->Get("hname"));
+    if (!fTreeS) {
+      fTreeS = new TTree(hname,"SDigits");
+      fTreeS->Write(0,TObject::kOverwrite);
+    }
+    cwd->cd();
+  }
+
+  if (oR) {
+    delete fTreeR;
+    sprintf(hname,"TreeR%d",fEvent);
+    file->cd();
+    fTreeR = static_cast<TTree*>(file->Get("hname"));
+    if (!fTreeR) {
+      fTreeR = new TTree(hname,"RecPoint");
+      fTreeR->Write(0,TObject::kOverwrite);
+    }
+    cwd->cd();
+  }
+}