Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / STEER / STEER / AliMC.cxx
index 83dae2f..335cf13 100644 (file)
@@ -34,6 +34,7 @@
 #include <TStopwatch.h>
 #include <TSystem.h>
 #include <TVirtualMC.h>
+#include <TMCVerbose.h>
 #include <TTree.h>
  
 #include "AliCDBEntry.h"
 #include "AliSimulation.h"
 #include "AliStack.h"
 #include "AliTrackReference.h"
+#include "AliTransportMonitor.h"
 
+using std::endl;
+using std::cout;
 ClassImp(AliMC)
 
 //_______________________________________________________________________
 AliMC::AliMC() :
   fGenerator(0),
+  fSaveRndmStatus(kFALSE),
+  fSaveRndmEventStatus(kFALSE),
+  fReadRndmStatus(kFALSE),
+  fUseMonitoring(kFALSE),
+  fRndmFileName("random.root"),
   fEventEnergy(0),
   fSummEnergy(0),
   fSum2Energy(0),
@@ -68,6 +77,7 @@ AliMC::AliMC() :
   fDecayPdg(0),
   fImedia(0),
   fTransParName("\0"),
+  fMonitor(0),
   fHitLists(0),
   fTmpTreeTR(0),
   fTmpFileTR(0),
@@ -83,6 +93,11 @@ AliMC::AliMC() :
 AliMC::AliMC(const char *name, const char *title) :
   TVirtualMCApplication(name, title),
   fGenerator(0),
+  fSaveRndmStatus(kFALSE),
+  fSaveRndmEventStatus(kFALSE),
+  fReadRndmStatus(kFALSE),
+  fUseMonitoring(kFALSE),
+  fRndmFileName("random.root"),
   fEventEnergy(0),
   fSummEnergy(0),
   fSum2Energy(0),
@@ -93,6 +108,7 @@ AliMC::AliMC(const char *name, const char *title) :
   fDecayPdg(0),
   fImedia(new TArrayI(1000)),
   fTransParName("\0"),
+  fMonitor(0),
   fHitLists(new TList()),
   fTmpTreeTR(0),
   fTmpFileTR(0),
@@ -114,6 +130,7 @@ AliMC::~AliMC()
   delete fGenerator;
   delete fImedia;
   delete fHitLists;
+  delete fMonitor;
   // Delete track references
 }
 
@@ -147,7 +164,7 @@ void  AliMC::ConstructGeometry()
        return;
       }
     }
-    gMC->SetRootGeometry();
+    TVirtualMC::GetMC()->SetRootGeometry();
   }else{
     // Create modules, materials, geometry
     if (!gGeoManager) new TGeoManager("ALICE", "ALICE geometry");
@@ -203,6 +220,320 @@ void  AliMC::ConstructOpGeometry()
   }  
 }
 
+#include <TPDGCode.h>
+//_______________________________________________________________________
+void  AliMC::AddParticles()
+{
+  //
+  // Add particles (not present in Geant3 or Geant4)
+  //
+  
+  // --------------------------------------------------------------------
+  // An example of adding a particle He5 with defined decay mode
+  // (TO BE REMOVED after a useful code is added)
+  
+  cout << "########## AliMC::AddParticles"  << endl;
+
+  //Hypertriton
+  TVirtualMC::GetMC()->DefineParticle(1010010030, "HyperTriton", kPTHadron, 2.99131 , 1.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 3, kFALSE);
+  //Anti-Hypertriton
+  TVirtualMC::GetMC()->DefineParticle(-1010010030, "AntiHyperTriton", kPTHadron, 2.99131 , 1.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 3, kFALSE);
+
+//Hyper hydrogen 4
+  TVirtualMC::GetMC()->DefineParticle(1010010040, "Hyperhydrog4", kPTHadron, 3.931 , 1.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE);
+  //Anti-Hyper hydrogen 4
+  TVirtualMC::GetMC()->DefineParticle(-1010010040, "AntiHyperhydrog4", kPTHadron, 3.931 , 1.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE);
+
+//Hyper helium 4
+  TVirtualMC::GetMC()->DefineParticle(1010020040, "Hyperhelium4", kPTHadron, 3.929 , 2.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE);
+  //Anti-Hyper helium 4
+  TVirtualMC::GetMC()->DefineParticle(-1010020040, "AntiHyperhelium4", kPTHadron, 3.929 , 2.0, 2.632e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 4, kFALSE);
+
+
+  //Lambda-Neutron 
+  TVirtualMC::GetMC()->DefineParticle(1010000020, "LambdaNeutron", kPTNeutron, 2.054 , 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
+
+  //Anti-Lambda-Neutron
+  TVirtualMC::GetMC()->DefineParticle(-1010000020, "AntiLambdaNeutron", kPTNeutron, 2.054 , 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
+
+  //H-Dibaryon
+  TVirtualMC::GetMC()->DefineParticle(1020000020, "Hdibaryon", kPTNeutron, 2.21, 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
+  //Anti-H-Dibaryon
+  TVirtualMC::GetMC()->DefineParticle(-1020000020, "AntiHdibaryon", kPTNeutron, 2.21  , 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
+
+  //Xi-Proton
+  TVirtualMC::GetMC()->DefineParticle(1030000020, "Xi0Proton", kPTHadron, 2.248 , 1.0, 1.333e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
+
+  //Anti-Xi-Proton
+  TVirtualMC::GetMC()->DefineParticle(-1030000020, "AntiXi0Proton", kPTHadron, 2.248 , 1.0, 1.333e-10,"Ion", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
+  
+  //Lambda-Neutron-Neutron
+  TVirtualMC::GetMC()->DefineParticle(1010000030, "LambdaNeutronNeutron", kPTNeutron, 2.982 , 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
+  
+  //Anti-Lambda-Neutron-Neutron
+  TVirtualMC::GetMC()->DefineParticle(-1010000030, "AntiLambdaNeutronNeutron", kPTNeutron, 2.982 , 0.0, 2.632e-10,"Hadron", 0.0, 0, 1, 0, 0, 0, 0, 0, 2, kFALSE);
+  
+  
+  // Define the 2- and 3-body phase space decay for the Hyper-Triton
+  Int_t mode[6][3];                  
+  Float_t bratio[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+     bratio[kz] = 0.;
+     mode[kz][0] = 0;
+     mode[kz][1] = 0;
+     mode[kz][2] = 0;
+  }
+  bratio[0] = 50.;
+  mode[0][0] = 1000020030; // Helium3 
+  mode[0][1] = -211; // negative pion
+  
+  bratio[1] = 50.;
+  mode[1][0] = 1000010020; // deuteron 
+  mode[1][1] = 2212; // proton
+  mode[1][2] = -211; // negative pion
+
+  TVirtualMC::GetMC()->SetDecayMode(1010010030,bratio,mode);
+
+
+
+  // Define the 2- and 3-body phase space decay for the Anti-Hyper-Triton
+  Int_t amode[6][3];                  
+  Float_t abratio[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+     abratio[kz] = 0.;
+     amode[kz][0] = 0;
+     amode[kz][1] = 0;
+     amode[kz][2] = 0;
+  }
+  abratio[0] = 50.;
+  amode[0][0] = -1000020030; // anti- Helium3 
+  amode[0][1] = 211; // positive pion
+  abratio[1] = 50.;
+  amode[1][0] = -1000010020; // anti-deuteron 
+  amode[1][1] = -2212; // anti-proton
+  amode[1][2] = 211; // positive pion
+
+  TVirtualMC::GetMC()->SetDecayMode(-1010010030,abratio,amode);
+  
+  ////// ----------Hypernuclei with Mass=4 ----------- //////////
+  
+   // Define the 2- and 3-body phase space decay for the Hyper Hydrogen 4
+   
+  Int_t mode3[6][3];                  
+  Float_t bratio3[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+     bratio3[kz] = 0.;
+     mode3[kz][0] = 0;
+     mode3[kz][1] = 0;
+     mode3[kz][2] = 0;
+  }
+  bratio3[0] = 50.;
+  mode3[0][0] = 1000020040; // Helium4 
+  mode3[0][1] = -211; // negative pion
+  
+  bratio3[1] = 50.;
+  mode3[1][0] = 1000010030; // tritium
+  mode3[1][1] = 2212; // proton
+  mode3[1][2] = -211; // negative pion
+
+  TVirtualMC::GetMC()->SetDecayMode(1010010040,bratio3,mode3);
+
+
+  // Define the 2- and 3-body phase space decay for the Hyper Hydrogen 4
+  Int_t amode3[6][3];                  
+  Float_t abratio3[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+     abratio3[kz] = 0.;
+     amode3[kz][0] = 0;
+     amode3[kz][1] = 0;
+     amode3[kz][2] = 0;
+  }
+  abratio3[0] = 50.;
+  amode3[0][0] = -1000020040; // anti- Helium4 
+  amode3[0][1] = 211; // positive pion
+  abratio3[1] = 50.;
+  amode3[1][0] = -1000010030; // anti-tritium
+  amode3[1][1] = -2212; // anti-proton
+  amode3[1][2] = 211; // positive pion
+
+  TVirtualMC::GetMC()->SetDecayMode(-1010010040,abratio3,amode3);
+  
+  
+   // Define the 3-body phase space decay for the Hyper Helium 4
+  Int_t mode4[6][3];                  
+  Float_t bratio4[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+     bratio4[kz] = 0.;
+     mode4[kz][0] = 0;
+     mode4[kz][1] = 0;
+     mode4[kz][2] = 0;
+  }
+  bratio4[0] = 100.;
+  mode4[0][0] = 1000020030; // Helium3 
+  mode4[0][1] = -211; // negative pion
+  mode4[0][2] = 2212; // proton
+  
+  TVirtualMC::GetMC()->SetDecayMode(1010020040,bratio4,mode4);
+
+
+  // Define the 2-body phase space decay for the Anti-Hyper Helium 4
+  Int_t amode4[6][3];                  
+  Float_t abratio4[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+     abratio4[kz] = 0.;
+     amode4[kz][0] = 0;
+     amode4[kz][1] = 0;
+     amode4[kz][2] = 0;
+  }
+  abratio4[0] = 100.;
+  amode4[0][0] = -1000020030; // anti-Helium 3
+  amode4[0][1] = 211; // positive pion
+  amode4[0][2] = -2212; // anti proton
+
+  TVirtualMC::GetMC()->SetDecayMode(-1010020040,abratio4,amode4);
+
+  
+  // Define the 2-body phase space decay for the Lambda-neutron boundstate
+  Int_t mode1[6][3];                  
+  Float_t bratio1[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+     bratio1[kz] = 0.;
+     mode1[kz][0] = 0;
+     mode1[kz][1] = 0;
+     mode1[kz][2] = 0;
+  }
+  bratio1[0] = 100.;
+  mode1[0][0] = 1000010020; // deuteron 
+  mode1[0][1] = -211; // negative pion
+
+  TVirtualMC::GetMC()->SetDecayMode(1010000020,bratio1,mode1);
+
+
+  // Define the 2-body phase space decay for the Anti-Lambda-neutron boundstate
+  Int_t amode1[6][3];                  
+  Float_t abratio1[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+     abratio1[kz] = 0.;
+     amode1[kz][0] = 0;
+     amode1[kz][1] = 0;
+     amode1[kz][2] = 0;
+  }
+  abratio1[0] = 100.;
+  amode1[0][0] = -1000010020; // anti-deuteron 
+  amode1[0][1] = 211; // positive pion
+
+  TVirtualMC::GetMC()->SetDecayMode(-1010000020,abratio1,amode1);
+
+  // Define the 2-body phase space decay for the H-Dibaryon
+  Int_t mode2[6][3];                  
+  Float_t bratio2[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+     bratio2[kz] = 0.;
+     mode2[kz][0] = 0;
+     mode2[kz][1] = 0;
+     mode2[kz][2] = 0;
+  }
+  bratio2[0] = 100.;
+  mode2[0][0] = 3122; // Lambda 
+  mode2[0][1] = 2212; // proton
+  mode2[0][2] = -211; // negative pion
+
+  TVirtualMC::GetMC()->SetDecayMode(1020000020,bratio2,mode2);
+
+  // Define the 2-body phase space decay for the Anti-H-Dibaryon
+  Int_t amode2[6][3];                  
+  Float_t abratio2[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+     abratio2[kz] = 0.;
+     amode2[kz][0] = 0;
+     amode2[kz][1] = 0;
+     amode2[kz][2] = 0;
+  }
+  abratio2[0] = 100.;
+  amode2[0][0] = -3122; // anti-deuteron 
+  amode2[0][1] = -2212; // anti-proton
+  amode2[0][2] = 211; // positive pion
+
+  TVirtualMC::GetMC()->SetDecayMode(-1020000020,abratio2,amode2);
+
+  // Define the 2-body phase space decay for the Xi0P
+  Int_t mode5[6][3];
+  Float_t bratio5[6];
+  
+  for (Int_t kz = 0; kz < 6; kz++) {
+    bratio5[kz] = 0.;
+    mode5[kz][0] = 0;
+    mode5[kz][1] = 0;
+    mode5[kz][2] = 0;
+  }
+  bratio5[0] = 100.;
+  mode5[0][0] = 3122; // Lambda
+  mode5[0][1] = 2212; // proton
+  
+  TVirtualMC::GetMC()->SetDecayMode(1030000020,bratio5,mode5);
+  
+  // Define the 2-body phase space decay for the Anti-Xi0P
+  Int_t amode5[6][3];
+  Float_t abratio5[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+    abratio5[kz] = 0.;
+    amode5[kz][0] = 0;
+    amode5[kz][1] = 0;
+    amode5[kz][2] = 0;
+  }
+  abratio5[0] = 100.;
+  amode5[0][0] = -3122; // anti-Lambda
+  amode5[0][1] = -2212; // anti-proton
+  
+  TVirtualMC::GetMC()->SetDecayMode(-1030000020,abratio5,amode5);
+  
+  // Define the 2-body phase space decay for the Lambda-Neutron-Neutron
+  Int_t mode6[6][3];
+  Float_t bratio6[6];
+
+  for (Int_t kz = 0; kz < 6; kz++) {
+    bratio6[kz] = 0.;
+    mode6[kz][0] = 0;
+    mode6[kz][1] = 0;
+    mode6[kz][2] = 0;
+  }
+  bratio6[0] = 100.;
+  mode6[0][0] = 1000010030; // triton
+  mode6[0][1] = -211; // pion
+  
+  TVirtualMC::GetMC()->SetDecayMode(1010000030,bratio6,mode6);
+  
+  // Define the 2-body phase space decay for the Anti-Lambda-Neutron-Neutron
+  Int_t amode6[6][3];
+  Float_t abratio6[6];
+  
+  for (Int_t kz = 0; kz < 6; kz++) {
+    abratio6[kz] = 0.;
+    amode6[kz][0] = 0;
+    amode6[kz][1] = 0;
+    amode6[kz][2] = 0;
+  }
+  abratio6[0] = 100.;
+  amode6[0][0] = -1000010030; // anti-triton
+  amode6[0][1] = 211; // pion
+  
+  TVirtualMC::GetMC()->SetDecayMode(-1010000030,abratio6,amode6);
+
+  // end of the example
+  // --------------------------------------------------------------------
+}  
+  
 //_______________________________________________________________________
 void  AliMC::InitGeometry()
 { 
@@ -300,6 +631,12 @@ void AliMC::FinishRun()
   // Clean generator information
   AliDebug(1, "fGenerator->FinishRun()");
   fGenerator->FinishRun();
+  
+  // Monitoring information
+  if (fMonitor) {
+    fMonitor->Print();
+    fMonitor->Export("timing.root");
+  }  
 
   //Output energy summary tables
   AliDebug(1, "EnergySummary()");
@@ -323,6 +660,8 @@ void AliMC::PreTrack()
 {
   // Actions before the track's transport
 
+     //verbose.PreTrack();
+
      TObjArray &dets = *gAlice->Modules();
      AliModule *module;
 
@@ -337,21 +676,40 @@ void AliMC::Stepping()
   //
   // Called at every step during transport
   //
-  Int_t id = DetFromMate(gMC->CurrentMedium());
+  //verbose.Stepping();
+
+  Int_t id = DetFromMate(TVirtualMC::GetMC()->CurrentMedium());
   if (id < 0) return;
 
 
-  if ( gMC->IsNewTrack()            && 
-       gMC->TrackTime() == 0.       &&
+  if ( TVirtualMC::GetMC()->IsNewTrack()            && 
+       TVirtualMC::GetMC()->TrackTime() == 0.       &&
        fRDecayMin >= 0.             &&  
        fRDecayMax > fRDecayMin      &&
-       gMC->TrackPid() == fDecayPdg ) 
+       TVirtualMC::GetMC()->TrackPid() == fDecayPdg ) 
   {
       FixParticleDecaytime();
   } 
     
-
-  
+  // --- If monitoring timing was requested, monitor the step
+  if (fUseMonitoring) {
+    if (!fMonitor) {
+      fMonitor = new AliTransportMonitor(TVirtualMC::GetMC()->NofVolumes()+1);
+      fMonitor->Start();
+    }  
+    if (TVirtualMC::GetMC()->IsNewTrack() || TVirtualMC::GetMC()->TrackTime() == 0. || TVirtualMC::GetMC()->TrackStep()<1.1E-10) {
+      fMonitor->DummyStep();
+    } else {
+    // Normal stepping
+      Int_t copy;
+      Int_t volId = TVirtualMC::GetMC()->CurrentVolID(copy);
+      Int_t pdg = TVirtualMC::GetMC()->TrackPid();
+      TLorentzVector xyz, pxpypz;
+      TVirtualMC::GetMC()->TrackPosition(xyz);
+      TVirtualMC::GetMC()->TrackMomentum(pxpypz);
+      fMonitor->StepInfo(volId, pdg, pxpypz.E(), xyz.X(), xyz.Y(), xyz.Z());
+    }  
+  }
   //
   // --- If lego option, do it and leave 
   if (AliSimulation::Instance()->Lego())
@@ -359,12 +717,12 @@ void AliMC::Stepping()
   else {
     Int_t copy;
     //Update energy deposition tables
-    AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
+    AddEnergyDeposit(TVirtualMC::GetMC()->CurrentVolID(copy),TVirtualMC::GetMC()->Edep());
     //
     // write tracke reference for track which is dissapearing - MI
 
-    if (gMC->IsTrackDisappeared() && !(gMC->IsTrackAlive())) {      
-       if (gMC->Etot() > 0.05) AddTrackReference(GetCurrentTrackNumber(), 
+    if (TVirtualMC::GetMC()->IsTrackDisappeared() && !(TVirtualMC::GetMC()->IsTrackAlive())) {      
+       if (TVirtualMC::GetMC()->Etot() > 0.05) AddTrackReference(GetCurrentTrackNumber(), 
                                                AliTrackReference::kDisappeared);
        
 
@@ -417,7 +775,7 @@ void AliMC::EnergySummary()
       for(i=0;i<(3<left?3:left);i++) {
        j=kn*3+i;
         id=Int_t (fEventEnergy[j]+0.1);
-       printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
+       printf(" %s %10.3f +- %10.3f%%;",TVirtualMC::GetMC()->VolName(id),fSummEnergy[j],fSum2Energy[j]);
       }
       printf("\n");
     }
@@ -430,7 +788,7 @@ void AliMC::EnergySummary()
       for(i=0;i<(5<left?5:left);i++) {
        j=kn*5+i;
         id=Int_t (fEventEnergy[j]+0.1);
-       printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
+       printf(" %s %10.3f%%;",TVirtualMC::GetMC()->VolName(id),100*fSummEnergy[j]/edtot);
       }
       printf("\n");
     }
@@ -443,7 +801,7 @@ void AliMC::EnergySummary()
   //  fSummEnergy.Set(0);
   //  fSum2Energy.Set(0);
 }
-
+#include <TFile.h>
 //_____________________________________________________________________________
 void AliMC::BeginEvent()
 {
@@ -457,7 +815,7 @@ void AliMC::BeginEvent()
   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
-    
+
   AliRunLoader *runloader=AliRunLoader::Instance();
 
   /*******************************/    
@@ -479,15 +837,38 @@ void AliMC::BeginEvent()
   else
       runloader->MakeStack();//or make a new one
   
+  // Random engine status
+  //
   
+  if ( fSaveRndmStatus || fSaveRndmEventStatus) {
+    TString fileName="random";
+    if ( fSaveRndmEventStatus ) {
+      fileName += "Evt";
+      fileName += gAlice->GetEventNrInRun();
+    }
+    fileName += ".root";
+       
+    // write ROOT random engine status
+    cout << "Saving random engine status in " << fileName.Data() << endl;
+    TFile f(fileName.Data(),"RECREATE");
+    gRandom->Write(fileName.Data());
+  }     
+
+  if ( fReadRndmStatus ) {
+    //read ROOT random engine status
+    cout << "Reading random engine status from " << fRndmFileName.Data() << endl;
+    TFile f(fRndmFileName.Data());
+    gRandom->Read(fRndmFileName.Data());    
+  }       
+
   if(AliSimulation::Instance()->Lego() == 0x0)
   { 
       AliDebug(1, "fRunLoader->MakeTree(K)");
       runloader->MakeTree("K");
   }
   
-  AliDebug(1, "gMC->SetStack(fRunLoader->Stack())");
-  gMC->SetStack(runloader->Stack());//Was in InitMC - but was moved here 
+  AliDebug(1, "TVirtualMC::GetMC()->SetStack(fRunLoader->Stack())");
+  TVirtualMC::GetMC()->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
@@ -589,7 +970,7 @@ void AliMC::FinishPrimary()
   //  const Int_t times=10;
   // This primary is finished, purify stack
 #if ROOT_VERSION_CODE > 262152
-  if (!(gMC->SecondariesAreOrdered())) {
+  if (!(TVirtualMC::GetMC()->SecondariesAreOrdered())) {
       if (runloader->Stack()->ReorderKine()) RemapHits();
   }
 #endif
@@ -681,6 +1062,7 @@ void AliMC::FinishEvent()
   // Update Header information 
   header->SetNprimary(stack->GetNprimary());
   header->SetNtrack(stack->GetNtrack());  
+  header->SetTimeStamp(AliSimulation::Instance()->GenerateTimeStamp());
 
   // Write out the kinematics
   if (!AliSimulation::Instance()->Lego()) stack->FinishEvent();
@@ -721,8 +1103,9 @@ void AliMC::Init()
 {
   // MC initialization
 
+
    //=================Create Materials and geometry
-   gMC->Init();
+   TVirtualMC::GetMC()->Init();
   // Set alignable volumes for the whole geometry (with old root)
 #if ROOT_VERSION_CODE < 331527
   SetAllAlignableVolumes();
@@ -733,16 +1116,15 @@ void AliMC::Init()
    MediaTable();
 
    //Compute cross-sections
-   gMC->BuildPhysics();
+   TVirtualMC::GetMC()->BuildPhysics();
    
    //Initialise geometry deposition table
-   fEventEnergy.Set(gMC->NofVolumes()+1);
-   fSummEnergy.Set(gMC->NofVolumes()+1);
-   fSum2Energy.Set(gMC->NofVolumes()+1);
+   fEventEnergy.Set(TVirtualMC::GetMC()->NofVolumes()+1);
+   fSummEnergy.Set(TVirtualMC::GetMC()->NofVolumes()+1);
+   fSum2Energy.Set(TVirtualMC::GetMC()->NofVolumes()+1);
 
    // Register MC in configuration 
-   AliConfig::Instance()->Add(gMC);
-
+   AliConfig::Instance()->Add(TVirtualMC::GetMC());
 }
 
 //_______________________________________________________________________
@@ -827,7 +1209,7 @@ void AliMC::ReadTransPar()
   char* filtmp;
   Float_t cut[kncuts];
   Int_t flag[knflags];
-  Int_t i, itmed, iret, ktmed, kz;
+  Int_t i, itmed, iret, jret, ktmed, kz;
   FILE *lun;
   //
   // See whether the file is there
@@ -853,7 +1235,7 @@ void AliMC::ReadTransPar()
       return;
     }
     // Read the end of line
-    fscanf(lun,"%*c");
+    jret = fscanf(lun,"%*c");
     if(!iret) continue;
     if(line[0]=='*') continue;
     // Read the numbers
@@ -884,7 +1266,7 @@ void AliMC::ReadTransPar()
          if(cut[kz]>=0) {
            AliDebug(2, Form("%-6s set to %10.3E for tracking medium code %4d for %s",
                             kpars[kz],cut[kz],itmed,mod->GetName()));
-           gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
+           TVirtualMC::GetMC()->Gstpar(ktmed,kpars[kz],cut[kz]);
          }
        }
        // Set transport mechanisms
@@ -892,7 +1274,7 @@ void AliMC::ReadTransPar()
          if(flag[kz]>=0) {
            AliDebug(2, Form("%-6s set to %10d for tracking medium code %4d for %s",
                             kpars[kncuts+kz],flag[kz],itmed,mod->GetName()));
-           gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
+           TVirtualMC::GetMC()->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
          }
        }
       } else {
@@ -1144,7 +1526,7 @@ void AliMC::FixParticleDecaytime()
     //
 
     TLorentzVector p;
-    gMC->TrackMomentum(p);
+    TVirtualMC::GetMC()->TrackMomentum(p);
     Double_t tmin, tmax;
     Double_t b;
 
@@ -1183,7 +1565,7 @@ void AliMC::FixParticleDecaytime()
     //
     // Force decay time in transport code
     //
-    gMC->ForceDecayTime(t / 2.99792458e10);
+    TVirtualMC::GetMC()->ForceDecayTime(t / 2.99792458e10);
 }
 
 void AliMC::MakeTmpTrackRefsTree()