]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliMC.cxx
Increasing the class version number.
[u/mrichter/AliRoot.git] / STEER / AliMC.cxx
index 9651d78d2e849a9865abff2a0789855c0e3ca2bb..f18473ffadeefdc00988a945c8575dbc9472cc77 100644 (file)
 // Author: F.Carminati
 //         Federico.Carminati@cern.ch
 
+#include <string.h>
+
 #include <RVersion.h>
-#include <TBrowser.h>
+#include <TArrayI.h>
 #include <TClonesArray.h>
+#include <TFile.h>
+#include <TGeoGlobalMagField.h>
 #include <TGeoManager.h>
+#include <TParticle.h>
+#include <TROOT.h>
 #include <TStopwatch.h>
 #include <TSystem.h>
 #include <TVirtualMC.h>
-#include <TParticle.h>
-#include <TROOT.h>
-#include <TFile.h>
  
-#include "AliLog.h"
+#include "AliCDBEntry.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
 #include "AliDetector.h"
 #include "AliGenerator.h"
+#include "AliGeomManager.h"
 #include "AliHeader.h"
+#include "AliHit.h"
 #include "AliLego.h"
+#include "AliLog.h"
 #include "AliMC.h"
-#include "AliMCQA.h"
+#include "AliMagF.h"
 #include "AliRun.h"
-#include "AliHit.h"
+#include "AliSimulation.h"
 #include "AliStack.h"
-#include "AliMagF.h"
 #include "AliTrackReference.h"
-#include "AliSimulation.h"
-#include "AliGeomManager.h"
-#include "AliCDBManager.h"
-#include "AliCDBStorage.h"
-#include "AliCDBEntry.h"
-
 
 ClassImp(AliMC)
 
@@ -66,12 +67,11 @@ AliMC::AliMC() :
   fDecayPdg(0),
   fImedia(0),
   fTransParName("\0"),
-  fMCQA(0),
   fHitLists(0),
   fTmpTreeTR(0),
   fTmpFileTR(0),
-  fTrackReferences(0),
-  fTmpTrackReferences(0)
+  fTrackReferences(),
+  fTmpTrackReferences()
 
 {
   //default constructor
@@ -92,12 +92,11 @@ AliMC::AliMC(const char *name, const char *title) :
   fDecayPdg(0),
   fImedia(new TArrayI(1000)),
   fTransParName("\0"),
-  fMCQA(0),
   fHitLists(new TList()),
   fTmpTreeTR(0),
   fTmpFileTR(0),
-  fTrackReferences(new TClonesArray("AliTrackReference", 100)),
-  fTmpTrackReferences(new TClonesArray("AliTrackReference", 100))
+  fTrackReferences("AliTrackReference", 100),
+  fTmpTrackReferences("AliTrackReference", 100)
 {
   //constructor
   // Set transport parameters
@@ -107,61 +106,14 @@ AliMC::AliMC(const char *name, const char *title) :
   for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
 }
 
-//_______________________________________________________________________
-AliMC::AliMC(const AliMC &mc) :
-  TVirtualMCApplication(mc),
-  fGenerator(0),
-  fEventEnergy(0),
-  fSummEnergy(0),
-  fSum2Energy(0),
-  fTrRmax(1.e10),
-  fTrZmax(1.e10),
-  fRDecayMax(1.e10),
-  fRDecayMin(-1.),
-  fDecayPdg(0),
-  fImedia(0),
-  fTransParName("\0"),
-  fMCQA(0),
-  fHitLists(0),
-  fTmpTreeTR(0),
-  fTmpFileTR(0),
-  fTrackReferences(0),
-  fTmpTrackReferences(0)
-{
-  //
-  // Copy constructor for AliMC
-  //
-  mc.Copy(*this);
-}
-
 //_______________________________________________________________________
 AliMC::~AliMC()
 {
   //destructor
   delete fGenerator;
   delete fImedia;
-  delete fMCQA;
   delete fHitLists;
   // Delete track references
-  if (fTrackReferences) {
-    fTrackReferences->Delete();
-    delete fTrackReferences;
-    fTrackReferences     = 0;
-  }
-
-if (fTmpTrackReferences) {
-    fTmpTrackReferences->Delete();
-    delete fTmpTrackReferences;
-    fTmpTrackReferences     = 0;
-  }
-
-}
-
-//_______________________________________________________________________
-void AliMC::Copy(TObject &) const
-{
-  //dummy Copy function
-  AliFatal("Not implemented!");
 }
 
 //_______________________________________________________________________
@@ -174,8 +126,8 @@ void  AliMC::ConstructGeometry()
   // at InitGeometry().
   //
 
-  if(gAlice->IsRootGeometry()){ //load geometry either from CDB or from file
-    if(gAlice->IsGeomFromCDB()){
+  if(AliSimulation::Instance()->IsGeometryFromFile()){ //load geometry either from CDB or from file
+    if(IsGeometryFromCDB()){
       AliInfo("Loading geometry from CDB default storage");
       AliCDBPath path("GRP","Geometry","Data");
       AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
@@ -185,7 +137,7 @@ void  AliMC::ConstructGeometry()
       if (!gGeoManager) AliFatal("TGeoManager object not found in the specified CDB entry!");
     }else{
       // Load geometry
-      const char *geomfilename = gAlice->GetGeometryFileName();
+      const char *geomfilename = AliSimulation::Instance()->GetGeometryFile();
       if(gSystem->ExpandPathName(geomfilename)){
        AliInfo(Form("Loading geometry from file:\n %40s",geomfilename));
        TGeoManager::Import(geomfilename);
@@ -194,8 +146,10 @@ void  AliMC::ConstructGeometry()
        return;
       }
     }
+    gMC->SetRootGeometry();
   }else{
     // Create modules, materials, geometry
+    if (!gGeoManager) new TGeoManager("ALICE", "ALICE geometry");
     TStopwatch stw;
     TIter next(gAlice->Modules());
     AliModule *detector;
@@ -215,16 +169,19 @@ void  AliMC::ConstructGeometry()
 //_______________________________________________________________________
 Bool_t  AliMC::MisalignGeometry() 
 {
-// Call misalignment code if AliSimulation object was defined.
-
-   if(!gAlice->IsRootGeometry()){
-     //Set alignable volumes for the whole geometry
-     SetAllAlignableVolumes();
-   }
-   // Misalign geometry via AliSimulation instance
-   if (!AliSimulation::GetInstance()) return kFALSE;
-   AliGeomManager::SetGeometry(gGeoManager);
-   return AliSimulation::GetInstance()->MisalignGeometry(gAlice->GetRunLoader());
+  // Call misalignment code if AliSimulation object was defined.
+  
+  if(!AliSimulation::Instance()->IsGeometryFromFile()){
+    //Set alignable volumes for the whole geometry
+    SetAllAlignableVolumes();
+  }
+  // Misalign geometry via AliSimulation instance
+  if (!AliSimulation::Instance()) return kFALSE;
+  AliGeomManager::SetGeometry(gGeoManager);
+  if(!AliGeomManager::CheckSymNamesLUT("ALL"))
+    AliFatal("Current loaded geometry differs in the definition of symbolic names!");
+  
+  return AliSimulation::Instance()->MisalignGeometry(AliRunLoader::Instance());
 }   
 
 //_______________________________________________________________________
@@ -238,6 +195,8 @@ void  AliMC::ConstructOpGeometry()
   AliModule *detector;
   AliInfo("Optical properties definition");
   while((detector = dynamic_cast<AliModule*>(next()))) {
+    // Initialise detector geometry
+    if(AliSimulation::Instance()->IsGeometryFromFile()) detector->CreateMaterials();
     // Initialise detector optical properties
     detector->DefineOpticalProperties();
   }  
@@ -256,14 +215,30 @@ void  AliMC::InitGeometry()
   AliModule *detector;
   while((detector = dynamic_cast<AliModule*>(next()))) {
     stw.Start();
-    // Initialise detector geometry
-    if(gAlice->IsRootGeometry()) detector->CreateMaterials();
     detector->Init();
     AliInfo(Form("%10s R:%.2fs C:%.2fs",
                 detector->GetName(),stw.RealTime(),stw.CpuTime()));
   }
 }
 
+//_______________________________________________________________________
+void AliMC::SetGeometryFromCDB()
+{
+  // Set the loading of geometry from cdb instead of creating it
+  // A default CDB storage needs to be set before this method is called
+  if(AliCDBManager::Instance()->IsDefaultStorageSet() &&
+     AliCDBManager::Instance()->GetRun() >= 0)
+    AliSimulation::Instance()->SetGeometryFile("*OCDB*");
+  else
+    AliError("Loading of geometry from CDB ignored. First set a default CDB storage!");
+}
+
+//_______________________________________________________________________
+Bool_t AliMC::IsGeometryFromCDB() const
+{
+  return (strcmp(AliSimulation::Instance()->GetGeometryFile(),"*OCDB*")==0);
+}
+
 //_______________________________________________________________________
 void  AliMC::SetAllAlignableVolumes()
 { 
@@ -346,14 +321,13 @@ void AliMC::BeginPrimary()
 void AliMC::PreTrack()
 {
   // Actions before the track's transport
+
      TObjArray &dets = *gAlice->Modules();
      AliModule *module;
 
      for(Int_t i=0; i<=gAlice->GetNdets(); i++)
        if((module = dynamic_cast<AliModule*>(dets[i])))
         module->PreTrack();
-
-     fMCQA->PreTrack();
 }
 
 //_______________________________________________________________________
@@ -362,7 +336,6 @@ void AliMC::Stepping()
   //
   // Called at every step during transport
   //
-    
   Int_t id = DetFromMate(gMC->CurrentMedium());
   if (id < 0) return;
 
@@ -380,23 +353,25 @@ void AliMC::Stepping()
   
   //
   // --- If lego option, do it and leave 
-  if (gAlice->Lego())
-    gAlice->Lego()->StepManager();
+  if (AliSimulation::Instance()->Lego())
+    AliSimulation::Instance()->Lego()->StepManager();
   else {
     Int_t copy;
     //Update energy deposition tables
     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
     //
     // write tracke reference for track which is dissapearing - MI
-    if (gMC->IsTrackDisappeared()) {      
-       if (gMC->Etot()>0.05) AddTrackReference(GetCurrentTrackNumber(), 
+
+    if (gMC->IsTrackDisappeared() && !(gMC->IsTrackAlive())) {      
+       if (gMC->Etot() > 0.05) AddTrackReference(GetCurrentTrackNumber(), 
                                                AliTrackReference::kDisappeared);
+       
+
     }
-  
+
     //Call the appropriate stepping routine;
     AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
     if(det && det->StepManagerIsEnabled()) {
-      if(AliLog::GetGlobalDebugLevel()>0) fMCQA->StepManager(id);
       det->StepManager();
     }
   }
@@ -414,7 +389,7 @@ void AliMC::EnergySummary()
   Float_t ed, ed2;
   Int_t kn, i, left, j, id;
   const Float_t kzero=0;
-  Int_t ievent=gAlice->GetRunLoader()->GetHeader()->GetEvent()+1;
+  Int_t ievent=AliRunLoader::Instance()->GetHeader()->GetEvent()+1;
   //
   // Energy loss information
   if(ievent) {
@@ -482,7 +457,7 @@ void AliMC::BeginEvent()
   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
     
-  AliRunLoader *runloader=gAlice->GetRunLoader();
+  AliRunLoader *runloader=AliRunLoader::Instance();
 
   /*******************************/    
   /*   Clean after eventual      */
@@ -504,14 +479,14 @@ void AliMC::BeginEvent()
       runloader->MakeStack();//or make a new one
   
   
-  if(gAlice->Lego() == 0x0)
+  if(AliSimulation::Instance()->Lego() == 0x0)
   { 
       AliDebug(1, "fRunLoader->MakeTree(K)");
       runloader->MakeTree("K");
   }
   
   AliDebug(1, "gMC->SetStack(fRunLoader->Stack())");
-  gMC->SetStack(gAlice->GetRunLoader()->Stack());//Was in InitMC - but was moved here 
+  gMC->SetStack(runloader->Stack());//Was in InitMC - but was moved here 
                                      //because we don't have guarantee that 
                                      //stack pointer is not going to change from event to event
                         //since it bellobgs to header and is obtained via RunLoader
@@ -519,13 +494,13 @@ void AliMC::BeginEvent()
   //  Reset all Detectors & kinematics & make/reset trees
   //
     
-  runloader->GetHeader()->Reset(gAlice->GetRunNumber(),gAlice->GetEvNumber(),
+  runloader->GetHeader()->Reset(AliCDBManager::Instance()->GetRun(),gAlice->GetEvNumber(),
                                gAlice->GetEventNrInRun());
 //  fRunLoader->WriteKinematics("OVERWRITE");  is there any reason to rewrite here since MakeTree does so
 
-  if(gAlice->Lego()) 
+  if(AliSimulation::Instance()->Lego()) 
   {
-      gAlice->Lego()->BeginEvent();
+      AliSimulation::Instance()->Lego()->BeginEvent();
       return;
   }
   
@@ -562,10 +537,37 @@ void AliMC::ResetHits()
   }
 }
 
+//_______________________________________________________________________
+void AliMC::ResetDigits()
+{
+  //
+  //  Reset all Detectors digits
+  //
+  TIter next(gAlice->Modules());
+  AliModule *detector;
+  while((detector = dynamic_cast<AliModule*>(next()))) {
+     detector->ResetDigits();
+  }
+}
+
+//_______________________________________________________________________
+void AliMC::ResetSDigits()
+{
+  //
+  //  Reset all Detectors digits
+  //
+  TIter next(gAlice->Modules());
+  AliModule *detector;
+  while((detector = dynamic_cast<AliModule*>(next()))) {
+     detector->ResetSDigits();
+  }
+}
+
 //_______________________________________________________________________
 void AliMC::PostTrack()
 {
   // Posts tracks for each module
+
   TObjArray &dets = *gAlice->Modules();
   AliModule *module;
   
@@ -580,18 +582,17 @@ void AliMC::FinishPrimary()
   //
   // Called  at the end of each primary track
   //
-  AliRunLoader *runloader=gAlice->GetRunLoader();
+
+  AliRunLoader *runloader=AliRunLoader::Instance();
   //  static Int_t count=0;
   //  const Int_t times=10;
   // This primary is finished, purify stack
 #if ROOT_VERSION_CODE > 262152
   if (!(gMC->SecondariesAreOrdered())) {
-      runloader->Stack()->ReorderKine();
-      RemapHits();
+      if (runloader->Stack()->ReorderKine()) RemapHits();
   }
 #endif
-  runloader->Stack()->PurifyKine();
-  RemapHits();
+  if (runloader->Stack()->PurifyKine()) RemapHits();
   
   TIter next(gAlice->Modules());
   AliModule *detector;
@@ -613,7 +614,7 @@ void AliMC::RemapHits()
 {
 //    
 // Remaps the track labels of the hits
-    AliRunLoader *runloader=gAlice->GetRunLoader();
+    AliRunLoader *runloader=AliRunLoader::Instance();
     AliStack* stack = runloader->Stack();
     TList* hitLists = GetHitLists();
     TIter next(hitLists);
@@ -631,12 +632,14 @@ void AliMC::RemapHits()
     // This for detectors which have a special mapping mechanism
     // for hits, such as TPC and TRD
     //
+
     
     TObjArray* modules = gAlice->Modules();
     TIter nextmod(modules);
-    AliDetector *detector;
-    while((detector = dynamic_cast<AliDetector*>(nextmod()))) {
-       detector->RemapTrackHitIDs(stack->TrackLabelMap());
+    AliModule *module;
+    while((module = (AliModule*) nextmod())) {
+       AliDetector* det = dynamic_cast<AliDetector*> (module);
+       if (det) det->RemapTrackHitIDs(stack->TrackLabelMap());
     }
     //
     RemapTrackReferencesIDs(stack->TrackLabelMap());
@@ -648,10 +651,8 @@ void AliMC::FinishEvent()
   //
   // Called at the end of the event.
   //
-  
-  //
     
-  if(gAlice->Lego()) gAlice->Lego()->FinishEvent();
+  if(AliSimulation::Instance()->Lego()) AliSimulation::Instance()->Lego()->FinishEvent();
 
   TIter next(gAlice->Modules());
   AliModule *detector;
@@ -667,7 +668,7 @@ void AliMC::FinishEvent()
     fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
    }
 
-  AliRunLoader *runloader=gAlice->GetRunLoader();
+  AliRunLoader *runloader=AliRunLoader::Instance();
 
   AliHeader* header = runloader->GetHeader();
   AliStack* stack = runloader->Stack();
@@ -681,10 +682,10 @@ void AliMC::FinishEvent()
   header->SetNtrack(stack->GetNtrack());  
 
   // Write out the kinematics
-  if (!gAlice->Lego()) stack->FinishEvent();
+  if (!AliSimulation::Instance()->Lego()) stack->FinishEvent();
 
   // Synchronize the TreeTR with TreeK
-  ReorderAndExpandTreeTR();
+  if (fTmpTreeTR) ReorderAndExpandTreeTR();
    
   // Write out the event Header information
   TTree* treeE = runloader->TreeE();
@@ -698,7 +699,7 @@ void AliMC::FinishEvent()
     AliError("Can not get TreeE from RL");
    }
   
-  if(gAlice->Lego() == 0x0)
+  if(AliSimulation::Instance()->Lego() == 0x0)
    {
      runloader->WriteKinematics("OVERWRITE");
      runloader->WriteTrackRefs("OVERWRITE");
@@ -714,13 +715,6 @@ void AliMC::FinishEvent()
   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 }
 
-//_______________________________________________________________________
-void AliMC::Field(const Double_t* x, Double_t* b) const
-{
-  // Calculates field "b" at point "x"
-    gAlice->Field(x,b);
-}
-    
 //_______________________________________________________________________
 void AliMC::Init()
 {
@@ -745,9 +739,6 @@ void AliMC::Init()
    fSummEnergy.Set(gMC->NofVolumes()+1);
    fSum2Energy.Set(gMC->NofVolumes()+1);
 
-   //
-   fMCQA = new AliMCQA(gAlice->GetNdets());
-
    // Register MC in configuration 
    AliConfig::Instance()->Add(gMC);
 
@@ -824,12 +815,12 @@ void AliMC::ReadTransPar()
 
 
   const Int_t kncuts=10;
-  const Int_t knflags=11;
+  const Int_t knflags=12;
   const Int_t knpars=kncuts+knflags;
   const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
                               "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
                               "BREM","COMP","DCAY","DRAY","HADR","LOSS",
-                              "MULS","PAIR","PHOT","RAYL"};
+                              "MULS","PAIR","PHOT","RAYL","STRA"};
   char line[256];
   char detName[7];
   char* filtmp;
@@ -852,9 +843,9 @@ void AliMC::ReadTransPar()
     for(i=0;i<kncuts;i++) cut[i]=-99;
     for(i=0;i<knflags;i++) flag[i]=-99;
     itmed=0;
-    for(i=0;i<256;i++) line[i]='\0';
+    memset(line,0,256);
     // Read up to the end of line excluded
-    iret=fscanf(lun,"%[^\n]",line);
+    iret=fscanf(lun,"%255[^\n]",line);
     if(iret<0) {
       //End of file
       fclose(lun);
@@ -865,10 +856,10 @@ void AliMC::ReadTransPar()
     if(!iret) continue;
     if(line[0]=='*') continue;
     // Read the numbers
-    iret=sscanf(line,"%s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d",
+    iret=sscanf(line,"%s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d %d",
                detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
                &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
-               &flag[8],&flag[9],&flag[10]);
+               &flag[8],&flag[9],&flag[10],&flag[11]);
     if(!iret) continue;
     if(iret<0) {
       //reading error
@@ -923,18 +914,6 @@ void AliMC::SetTransPar(const char *filename)
   fTransParName = filename;
 }
 
-//_______________________________________________________________________
-void AliMC::Browse(TBrowser *b)
-{
-  //
-  // Called when the item "Run" is clicked on the left pane
-  // of the Root browser.
-  // It displays the Root Trees and all detectors.
-  //
-  //detectors are in folders anyway
-  b->Add(fMCQA,"AliMCQA");
-}
-
 //_______________________________________________________________________
 void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
 {
@@ -942,7 +921,7 @@ void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
   //  Add a hit to detector id
   //
   TObjArray &dets = *gAlice->Modules();
-  if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
+  if(dets[id]) static_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
 }
 
 //_______________________________________________________________________
@@ -952,7 +931,7 @@ void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
   // Add digit to detector id
   //
   TObjArray &dets = *gAlice->Modules();
-  if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
+  if(dets[id]) static_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
 }
 
 //_______________________________________________________________________
@@ -960,7 +939,7 @@ Int_t AliMC::GetCurrentTrackNumber() const {
   //
   // Returns current track
   //
-  return gAlice->GetRunLoader()->Stack()->GetCurrentTrackNumber();
+  return AliRunLoader::Instance()->Stack()->GetCurrentTrackNumber();
 }
 
 //_______________________________________________________________________
@@ -969,7 +948,7 @@ void AliMC::DumpPart (Int_t i) const
   //
   // Dumps particle i in the stack
   //
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
    if (runloader->Stack())
     runloader->Stack()->DumpPart(i);
 }
@@ -980,7 +959,7 @@ void AliMC::DumpPStack () const
   //
   // Dumps the particle stack
   //
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
    if (runloader->Stack())
     runloader->Stack()->DumpPStack();
 }
@@ -991,7 +970,7 @@ Int_t AliMC::GetNtrack() const {
   // Returns number of tracks in stack
   //
   Int_t ntracks = -1;
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
    if (runloader->Stack())
      ntracks = runloader->Stack()->GetNtrack();
    return ntracks;
@@ -1004,7 +983,7 @@ Int_t AliMC::GetPrimary(Int_t track) const
   // return number of primary that has generated track
   //
   Int_t nprimary = -999;
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
   if (runloader->Stack())
     nprimary = runloader->Stack()->GetPrimary(track);
   return nprimary;
@@ -1015,7 +994,7 @@ TParticle* AliMC::Particle(Int_t i) const
 {
   // Returns the i-th particle from the stack taking into account
   // the remaping done by PurifyKine
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
   if (runloader)
    if (runloader->Stack())
     return runloader->Stack()->Particle(i);
@@ -1023,11 +1002,11 @@ TParticle* AliMC::Particle(Int_t i) const
 }
 
 //_______________________________________________________________________
-TObjArray* AliMC::Particles() const {
+const TObjArray* AliMC::Particles() const {
   //
   // Returns pointer to Particles array
   //
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
   if (runloader)
    if (runloader->Stack())
     return runloader->Stack()->Particles();
@@ -1035,13 +1014,13 @@ TObjArray* AliMC::Particles() const {
 }
 
 //_______________________________________________________________________
-void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
-                      Float_t *vpos, Float_t *polar, Float_t tof,
+void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom,
+                      const Float_t *vpos, const Float_t *polar, Float_t tof,
                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
 { 
 // Delegate to stack
 //
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
   if (runloader)
     if (runloader->Stack())
       runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
@@ -1057,7 +1036,7 @@ void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg,
 { 
   // Delegate to stack
   //
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
   if (runloader)
     if (runloader->Stack())
       runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
@@ -1069,7 +1048,7 @@ void AliMC::SetHighWaterMark(Int_t nt) const
 {
     //
     // Set high water mark for last track in event
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
   if (runloader)
     if (runloader->Stack())
       runloader->Stack()->SetHighWaterMark(nt);
@@ -1081,7 +1060,7 @@ void AliMC::KeepTrack(Int_t track) const
   //
   // Delegate to stack
   //
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
   if (runloader)
     if (runloader->Stack())
       runloader->Stack()->KeepTrack(track);
@@ -1092,7 +1071,7 @@ void AliMC::FlagTrack(Int_t track) const
 {
   // Delegate to stack
   //
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
   if (runloader)
     if (runloader->Stack())
       runloader->Stack()->FlagTrack(track);
@@ -1104,7 +1083,7 @@ void AliMC::SetCurrentTrack(Int_t track) const
   //
   // Set current track number
   //
-  AliRunLoader * runloader = gAlice->GetRunLoader();
+  AliRunLoader * runloader = AliRunLoader::Instance();
   if (runloader)
     if (runloader->Stack())
       runloader->Stack()->SetCurrentTrack(track); 
@@ -1115,16 +1094,11 @@ AliTrackReference*  AliMC::AddTrackReference(Int_t label, Int_t id)
 {
   //
   // add a trackrefernce to the list
-  if (!fTrackReferences) {
-    AliError("Container trackrefernce not active");
-    return NULL;
-  }
   Int_t primary = GetPrimary(label);
   Particle(primary)->SetBit(kKeepBit);
 
-  Int_t nref = fTmpTrackReferences->GetEntriesFast();
-  TClonesArray &lref = *fTmpTrackReferences;
-  return new(lref[nref]) AliTrackReference(label, id);
+  Int_t nref = fTmpTrackReferences.GetEntriesFast();
+  return new(fTmpTrackReferences[nref]) AliTrackReference(label, id);
 }
 
 
@@ -1135,31 +1109,33 @@ void AliMC::ResetTrackReferences()
   //
   //  Reset all  references
   //
-  if (fTmpTrackReferences)   fTmpTrackReferences->Clear();
+    fTmpTrackReferences.Clear();
 }
 
-void AliMC::RemapTrackReferencesIDs(Int_t *map)
+//_______________________________________________________________________
+void AliMC::RemapTrackReferencesIDs(const Int_t *map)
 {
   // 
   // Remapping track reference
   // Called at finish primary
   //
-  if (!fTmpTrackReferences) return;
-  Int_t nEntries = fTmpTrackReferences->GetEntries();
+    
+  Int_t nEntries = fTmpTrackReferences.GetEntries();
   for (Int_t i=0; i < nEntries; i++){
-      AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTmpTrackReferences->UncheckedAt(i));
+      AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTmpTrackReferences.UncheckedAt(i));
       if (ref) {
          Int_t newID = map[ref->GetTrack()];
          if (newID>=0) ref->SetTrack(newID);
          else {
              ref->SetBit(kNotDeleted,kFALSE);
-             fTmpTrackReferences->RemoveAt(i);  
+             fTmpTrackReferences.RemoveAt(i);  
          }      
       } // if ref
   }
-  fTmpTrackReferences->Compress();
+  fTmpTrackReferences.Compress();
 }
 
+//_______________________________________________________________________
 void AliMC::FixParticleDecaytime()
 {
     //
@@ -1174,7 +1150,7 @@ void AliMC::FixParticleDecaytime()
     // Transverse velocity 
     Double_t vt    = p.Pt() / p.E();
     
-    if ((b = gAlice->Field()->SolenoidField()) > 0.) {     // [kG]
+    if ((b = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField()) > 0.) {     // [kG]
 
        // Radius of helix
        
@@ -1187,12 +1163,13 @@ void AliMC::FixParticleDecaytime()
        // Maximum and minimum decay time
        //
        // Check for curlers first
-       if (fRDecayMax * fRDecayMax / rho / rho / 2. > 1.) return;
+       const Double_t kOvRhoSqr2 = 1./(rho*TMath::Sqrt(2.));
+       if (fRDecayMax * kOvRhoSqr2 > 1.) return;
        
        //
  
-       tmax  = TMath::ACos(1. - fRDecayMax * fRDecayMax / rho / rho / 2.) / omega;   // [ct]
-       tmin  = TMath::ACos(1. - fRDecayMin * fRDecayMin / rho / rho / 2.) / omega;   // [ct]
+       tmax  = TMath::ACos((1.-fRDecayMax*kOvRhoSqr2)*(1.+fRDecayMax*kOvRhoSqr2)) / omega;   // [ct]
+       tmin  = TMath::ACos((1.-fRDecayMin*kOvRhoSqr2)*(1.+fRDecayMin*kOvRhoSqr2)) / omega;   // [ct]
     } else {
        tmax =  fRDecayMax / vt;                                                      // [ct] 
        tmin =  fRDecayMin / vt;                                                      // [ct]
@@ -1213,29 +1190,26 @@ void AliMC::MakeTmpTrackRefsTree()
     // Make the temporary track reference tree
     fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
     fTmpTreeTR = new TTree("TreeTR", "Track References");
-    if (!fTmpTrackReferences)  fTmpTrackReferences = new TClonesArray("AliTrackReference", 100);
-    fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTmpTrackReferences, 4000);
+    TClonesArray* pRef = &fTmpTrackReferences;
+    fTmpTreeTR->Branch("TrackReferences", &pRef, 4000);
 }
 
+//_______________________________________________________________________
 void AliMC::ReorderAndExpandTreeTR()
 {
 //
 //  Reorder and expand the temporary track reference tree in order to match the kinematics tree
 //
 
-    AliRunLoader *rl = gAlice->GetRunLoader();
+    AliRunLoader *rl = AliRunLoader::Instance();
 //
 //  TreeTR
     AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
     rl->MakeTrackRefsContainer(); 
     TTree * treeTR = rl->TreeTR();
-    if (treeTR){
        // make branch for central track references
-       if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference",0);
-       TBranch *branch;
-       branch = treeTR->Branch("TrackReferences",&fTrackReferences);
-       branch->SetAddress(&fTrackReferences);
-    }
+       TClonesArray* pRef = &fTrackReferences;
+       treeTR->Branch("TrackReferences", &pRef);
 
     AliStack* stack  = rl->Stack();
     Int_t np = stack->GetNprimary();
@@ -1246,16 +1220,15 @@ void AliMC::ReorderAndExpandTreeTR()
     Int_t it = 0;
     for (Int_t ip = np - 1; ip > -1; ip--) {
        TParticle *part = stack->Particle(ip);
-//     printf("Particle %5d %5d %5d %5d %5d \n", ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), part->GetLastDaughter());
+       //printf("Particle %5d %5d %5d %5d %5d \n", ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), part->GetLastDaughter());
        
        // Skip primaries that have not been transported
        Int_t dau1  = part->GetFirstDaughter();
        Int_t dau2  = -1;
-       // if ((dau1 > -1 && dau1 < np) || part->GetStatusCode() > 1) continue;
        if (!part->TestBit(kTransportBit)) continue;
        //
        fTmpTreeTR->GetEntry(it++);
-       Int_t nh = fTmpTrackReferences->GetEntries();
+       Int_t nh = fTmpTrackReferences.GetEntries();
        // Determine range of secondaries produced by this primary
        if (dau1 > -1) {
            Int_t inext = ip - 1;
@@ -1263,7 +1236,8 @@ void AliMC::ReorderAndExpandTreeTR()
                if (inext >= 0) {
                    part = stack->Particle(inext);
                    dau2 =  part->GetFirstDaughter();
-                   if (dau2 == -1 || dau2 < np) {
+                   if (!(part->TestBit(kTransportBit)) || dau2 == -1 || dau2 < np) {
+//                 if (dau2 == -1 || dau2 < np) {
                        dau2 = -1;
                    } else {
                        dau2--;
@@ -1275,12 +1249,11 @@ void AliMC::ReorderAndExpandTreeTR()
            } // find upper bound
        }  // dau2 < 0
 //     printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
-       
        // 
        // Loop over reference hits and find secondary label
        for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
            for (Int_t ih = 0; ih < nh; ih++) {
-               AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences->At(ih);
+               AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences.At(ih);
                Int_t label = tr->Label();
                // Skip primaries
                if (label == ip) continue;
@@ -1288,13 +1261,12 @@ void AliMC::ReorderAndExpandTreeTR()
                    AliWarning(Form("Track Reference Label out of range !: %5d %5d %5d \n", label, dau1, dau2));
                if (label == id) {
                    // secondary found
-                   Int_t nref =  fTrackReferences->GetEntriesFast();
-                   TClonesArray &lref = *fTrackReferences;
-                   new(lref[nref]) AliTrackReference(*tr);
+                   Int_t nref =  fTrackReferences.GetEntriesFast();
+                   new(fTrackReferences[nref]) AliTrackReference(*tr);
                }
            } // hits
            treeTR->Fill();
-           fTrackReferences->Clear();
+           fTrackReferences.Clear();
            ifills++;
        } // daughters
     } // tracks
@@ -1308,21 +1280,20 @@ void AliMC::ReorderAndExpandTreeTR()
        {
            // Skip particles that have not been transported
            fTmpTreeTR->GetEntry(it--);
-           Int_t nh = fTmpTrackReferences->GetEntries();
+           Int_t nh = fTmpTrackReferences.GetEntries();
            // 
            // Loop over reference hits and find primary labels
            for (Int_t ih = 0; ih < nh; ih++) {
-               AliTrackReference* tr = (AliTrackReference*)  fTmpTrackReferences->At(ih);
+               AliTrackReference* tr = (AliTrackReference*)  fTmpTrackReferences.At(ih);
                Int_t label = tr->Label();
                if (label == ip) {
-                   Int_t nref = fTrackReferences->GetEntriesFast();
-                   TClonesArray &lref = *fTrackReferences;
-                   new(lref[nref]) AliTrackReference(*tr);
+                   Int_t nref = fTrackReferences.GetEntriesFast();
+                   new(fTrackReferences[nref]) AliTrackReference(*tr);
                }
            } 
        }
        treeTR->Fill();
-       fTrackReferences->Clear();
+       fTrackReferences.Clear();
        ifills++;
     } // tracks
     // Check
@@ -1333,8 +1304,7 @@ void AliMC::ReorderAndExpandTreeTR()
     delete fTmpTreeTR;
     fTmpFileTR->Close();
     delete fTmpFileTR;
-    delete fTmpTrackReferences;
-    fTmpTrackReferences = 0;
+    fTmpTrackReferences.Clear();
     gSystem->Exec("rm -rf TrackRefsTmp.root");
 }