]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliMC.cxx
LUTs mapping symnames and original global matrices removed from AliGeomManager, which...
[u/mrichter/AliRoot.git] / STEER / AliMC.cxx
index defade3d1b2c4e5a31d766a07a847e278ff72c25..c7d6fe88915573b62d192fd4bef7a3710939f315 100644 (file)
 
 #include <RVersion.h>
 #include <TBrowser.h>
+#include <TClonesArray.h>
+#include <TGeoManager.h>
 #include <TStopwatch.h>
 #include <TSystem.h>
 #include <TVirtualMC.h>
+#include <TParticle.h>
+#include <TROOT.h>
+#include <TFile.h>
  
 #include "AliLog.h"
 #include "AliDetector.h"
 #include "AliMC.h"
 #include "AliMCQA.h"
 #include "AliRun.h"
+#include "AliHit.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)
@@ -49,14 +61,21 @@ AliMC::AliMC() :
   fSum2Energy(0),
   fTrRmax(1.e10),
   fTrZmax(1.e10),
+  fRDecayMax(1.e10),
+  fRDecayMin(-1.),
+  fDecayPdg(0),
   fImedia(0),
   fTransParName("\0"),
   fMCQA(0),
   fHitLists(0),
-  fTrackReferences(0)
+  fTmpTreeTR(0),
+  fTmpFileTR(0),
+  fTrackReferences(),
+  fTmpTrackReferences()
 
 {
   //default constructor
+  DecayLimits();
 }
 
 //_______________________________________________________________________
@@ -68,19 +87,24 @@ AliMC::AliMC(const char *name, const char *title) :
   fSum2Energy(0),
   fTrRmax(1.e10),
   fTrZmax(1.e10),
+  fRDecayMax(1.e10),
+  fRDecayMin(-1.),
+  fDecayPdg(0),
   fImedia(new TArrayI(1000)),
   fTransParName("\0"),
   fMCQA(0),
   fHitLists(new TList()),
-  fTrackReferences(new TClonesArray("AliTrackReference", 100))
+  fTmpTreeTR(0),
+  fTmpFileTR(0),
+  fTrackReferences("AliTrackReference", 100),
+  fTmpTrackReferences("AliTrackReference", 100)
 {
   //constructor
   // Set transport parameters
   SetTransPar();
-
+  DecayLimits();
   // Prepare the tracking medium lists
   for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
-
 }
 
 //_______________________________________________________________________
@@ -92,14 +116,20 @@ AliMC::AliMC(const AliMC &mc) :
   fSum2Energy(0),
   fTrRmax(1.e10),
   fTrZmax(1.e10),
+  fRDecayMax(1.e10),
+  fRDecayMin(-1.),
+  fDecayPdg(0),
   fImedia(0),
   fTransParName("\0"),
   fMCQA(0),
   fHitLists(0),
-  fTrackReferences(0)
+  fTmpTreeTR(0),
+  fTmpFileTR(0),
+  fTrackReferences(),
+  fTmpTrackReferences()
 {
   //
-  // Copy constructor for AliRun
+  // Copy constructor for AliMC
   //
   mc.Copy(*this);
 }
@@ -113,12 +143,6 @@ AliMC::~AliMC()
   delete fMCQA;
   delete fHitLists;
   // Delete track references
-  if (fTrackReferences) {
-    fTrackReferences->Delete();
-    delete fTrackReferences;
-    fTrackReferences     = 0;
-  }
-
 }
 
 //_______________________________________________________________________
@@ -132,9 +156,34 @@ void AliMC::Copy(TObject &) const
 void  AliMC::ConstructGeometry() 
 {
   //
-  // Create modules, materials, geometry
-  //
-
+  // Either load geometry from file or create it through usual
+  // loop on detectors. In the first case the method
+  // AliModule::CreateMaterials() only builds fIdtmed and is postponed
+  // at InitGeometry().
+  //
+
+  if(gAlice->IsRootGeometry()){ //load geometry either from CDB or from file
+    if(gAlice->IsGeomFromCDB()){
+      AliInfo("Loading geometry from CDB default storage");
+      AliCDBPath path("GRP","Geometry","Data");
+      AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
+      if(!entry) AliFatal("Unable to load geometry from CDB!");
+      entry->SetOwner(0);
+      gGeoManager = (TGeoManager*) entry->GetObject();
+      if (!gGeoManager) AliFatal("TGeoManager object not found in the specified CDB entry!");
+    }else{
+      // Load geometry
+      const char *geomfilename = gAlice->GetGeometryFileName();
+      if(gSystem->ExpandPathName(geomfilename)){
+       AliInfo(Form("Loading geometry from file:\n %40s",geomfilename));
+       TGeoManager::Import(geomfilename);
+      }else{
+       AliInfo(Form("Geometry file %40s not found!\n",geomfilename));
+       return;
+      }
+    }
+  }else{
+    // Create modules, materials, geometry
     TStopwatch stw;
     TIter next(gAlice->Modules());
     AliModule *detector;
@@ -147,28 +196,79 @@ void  AliMC::ConstructGeometry()
       AliInfo(Form("%10s R:%.2fs C:%.2fs",
                   detector->GetName(),stw.RealTime(),stw.CpuTime()));
     }
+  }
+  
+}
+
+//_______________________________________________________________________
+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);
+   if(!AliGeomManager::CheckSymNamesLUT())
+    AliFatal("Current loaded geometry differs in the definition of symbolic names!");
+
+   return AliSimulation::GetInstance()->MisalignGeometry(gAlice->GetRunLoader());
+}   
+
+//_______________________________________________________________________
+void  AliMC::ConstructOpGeometry() 
+{
+  //
+  // Loop all detector modules and call DefineOpticalProperties() method 
+  //
+
+  TIter next(gAlice->Modules());
+  AliModule *detector;
+  AliInfo("Optical properties definition");
+  while((detector = dynamic_cast<AliModule*>(next()))) {
+    // Initialise detector optical properties
+    detector->DefineOpticalProperties();
+  }  
 }
 
 //_______________________________________________________________________
 void  AliMC::InitGeometry()
 { 
   //
-  // Initialize detectors and display geometry
+  // Initialize detectors
   //
 
-   AliInfo("Initialisation:");
-    TStopwatch stw;
-    TIter next(gAlice->Modules());
-    AliModule *detector;
-    while((detector = dynamic_cast<AliModule*>(next()))) {
-      stw.Start();
-      // Initialise detector and display geometry
-      detector->Init();
-      detector->BuildGeometry();
-      AliInfo(Form("%10s R:%.2fs C:%.2fs",
-                  detector->GetName(),stw.RealTime(),stw.CpuTime()));
-    }
+  AliInfo("Initialisation:");
+  TStopwatch stw;
+  TIter next(gAlice->Modules());
+  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::SetAllAlignableVolumes()
+{ 
+  //
+  // Add alignable volumes (TGeoPNEntries) looping on all
+  // active modules
+  //
+
+  AliInfo(Form("Setting entries for all alignable volumes of active detectors"));
+  AliModule *detector;
+  TIter next(gAlice->Modules());
+  while((detector = dynamic_cast<AliModule*>(next()))) {
+    detector->AddAlignableVolumes();
+  }
 }
 
 //_______________________________________________________________________
@@ -196,13 +296,16 @@ void AliMC::ResetGenerator(AliGenerator *generator)
   //
   // Load the event generator
   //
-  if(fGenerator)
-    if(generator)
+  if(fGenerator) {
+    if(generator) {
       AliWarning(Form("Replacing generator %s with %s",
-                     fGenerator->GetName(),generator->GetName()))
-    else
+                     fGenerator->GetName(),generator->GetName()));
+    }
+    else {
       AliWarning(Form("Replacing generator %s with NULL",
                      fGenerator->GetName()));
+    }
+  }
   fGenerator = generator;
 }
 
@@ -228,7 +331,6 @@ void AliMC::BeginPrimary()
   // Reset Hits info
   ResetHits();
   ResetTrackReferences();
-
 }
 
 //_______________________________________________________________________
@@ -251,10 +353,21 @@ void AliMC::Stepping()
   //
   // Called at every step during transport
   //
-
-  Int_t id = DetFromMate(gMC->GetMedium());
+  Int_t id = DetFromMate(gMC->CurrentMedium());
   if (id < 0) return;
 
+
+  if ( gMC->IsNewTrack()            && 
+       gMC->TrackTime() == 0.       &&
+       fRDecayMin >= 0.             &&  
+       fRDecayMax > fRDecayMin      &&
+       gMC->TrackPid() == fDecayPdg ) 
+  {
+      FixParticleDecaytime();
+  } 
+    
+
+  
   //
   // --- If lego option, do it and leave 
   if (gAlice->Lego())
@@ -265,14 +378,18 @@ void AliMC::Stepping()
     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()) {
-      fMCQA->StepManager(id);
+      if(AliLog::GetGlobalDebugLevel()>0) fMCQA->StepManager(id);
       det->StepManager();
     }
   }
@@ -281,7 +398,7 @@ void AliMC::Stepping()
 //_______________________________________________________________________
 void AliMC::EnergySummary()
 {
-  //
+  //e
   // Print summary of deposited energy
   //
 
@@ -358,7 +475,7 @@ void AliMC::BeginEvent()
   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
     
-    AliRunLoader *runloader=gAlice->GetRunLoader();
+  AliRunLoader *runloader=gAlice->GetRunLoader();
 
   /*******************************/    
   /*   Clean after eventual      */
@@ -375,16 +492,17 @@ void AliMC::BeginEvent()
     // Clean detector information
   
   if (runloader->Stack())
-    runloader->Stack()->Reset();//clean stack -> tree is unloaded
+      runloader->Stack()->Reset();//clean stack -> tree is unloaded
   else
-    runloader->MakeStack();//or make a new one
+      runloader->MakeStack();//or make a new one
+  
   
   if(gAlice->Lego() == 0x0)
-   
-     AliDebug(1, "fRunLoader->MakeTree(K)");
-     runloader->MakeTree("K");
-   }
-   
+  { 
+      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 
                                      //because we don't have guarantee that 
@@ -394,15 +512,15 @@ void AliMC::BeginEvent()
   //  Reset all Detectors & kinematics & make/reset trees
   //
     
-  runloader->GetHeader()->Reset(gAlice->GetRunNumber(),gAlice->GetEvNumber(),
-    gAlice->GetEventNrInRun());
+  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()) 
-   {
-    gAlice->Lego()->BeginEvent();
-    return;
-   }
+  {
+      gAlice->Lego()->BeginEvent();
+      return;
+  }
   
 
   AliDebug(1, "ResetHits()");
@@ -411,32 +529,17 @@ void AliMC::BeginEvent()
   AliDebug(1, "fRunLoader->MakeTree(H)");
   runloader->MakeTree("H");
   
-  AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
-  runloader->MakeTrackRefsContainer();//for insurance
 
 
+  MakeTmpTrackRefsTree();
   //create new branches and SetAdresses
   TIter next(gAlice->Modules());
   AliModule *detector;
   while((detector = (AliModule*)next()))
    {
-    AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName()));
-    detector->MakeBranch("H"); 
-    AliDebug(2, Form("%s->MakeBranchTR()",detector->GetName()));
-    detector->MakeBranchTR();
-    AliDebug(2, Form("%s->SetTreeAddress()",detector->GetName()));
-    detector->SetTreeAddress();
+       AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName()));
+       detector->MakeBranch("H"); 
    }
-  // make branch for AliRun track References
-  TTree * treeTR = runloader->TreeTR();
-  if (treeTR){
-    // make branch for central track references
-    if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference",0);
-    TBranch *branch;
-    branch = treeTR->Branch("AliRun",&fTrackReferences);
-    branch->SetAddress(&fTrackReferences);
-  }
-  //
 }
 
 //_______________________________________________________________________
@@ -456,12 +559,12 @@ void AliMC::ResetHits()
 void AliMC::PostTrack()
 {
   // Posts tracks for each module
-     TObjArray &dets = *gAlice->Modules();
-     AliModule *module;
-
-     for(Int_t i=0; i<=gAlice->GetNdets(); i++)
-       if((module = dynamic_cast<AliModule*>(dets[i])))
-        module->PostTrack();
+  TObjArray &dets = *gAlice->Modules();
+  AliModule *module;
+  
+  for(Int_t i=0; i<=gAlice->GetNdets(); i++)
+    if((module = dynamic_cast<AliModule*>(dets[i])))
+      module->PostTrack();
 }
 
 //_______________________________________________________________________
@@ -475,10 +578,11 @@ void AliMC::FinishPrimary()
   //  const Int_t times=10;
   // This primary is finished, purify stack
 #if ROOT_VERSION_CODE > 262152
-  if (!(gMC->SecondariesAreOrdered()))
-       runloader->Stack()->ReorderKine();
+  if (!(gMC->SecondariesAreOrdered())) {
+      if (runloader->Stack()->ReorderKine()) RemapHits();
+  }
 #endif
-  runloader->Stack()->PurifyKine();
+  if (runloader->Stack()->PurifyKine()) RemapHits();
   
   TIter next(gAlice->Modules());
   AliModule *detector;
@@ -493,7 +597,42 @@ void AliMC::FinishPrimary()
   }
 
   // Write out track references if any
-  if (runloader->TreeTR()) runloader->TreeTR()->Fill();
+  if (fTmpTreeTR) fTmpTreeTR->Fill();
+}
+
+void AliMC::RemapHits()
+{
+//    
+// Remaps the track labels of the hits
+    AliRunLoader *runloader=gAlice->GetRunLoader();
+    AliStack* stack = runloader->Stack();
+    TList* hitLists = GetHitLists();
+    TIter next(hitLists);
+    TCollection *hitList;
+    
+    while((hitList = dynamic_cast<TCollection*>(next()))) {
+       TIter nexthit(hitList);
+       AliHit *hit;
+       while((hit = dynamic_cast<AliHit*>(nexthit()))) {
+           hit->SetTrack(stack->TrackLabel(hit->GetTrack()));
+       }
+    }
+    
+    // 
+    // This for detectors which have a special mapping mechanism
+    // for hits, such as TPC and TRD
+    //
+
+    
+    TObjArray* modules = gAlice->Modules();
+    TIter nextmod(modules);
+    AliModule *module;
+    while((module = (AliModule*) nextmod())) {
+       AliDetector* det = dynamic_cast<AliDetector*> (module);
+       if (det) det->RemapTrackHitIDs(stack->TrackLabelMap());
+    }
+    //
+    RemapTrackReferencesIDs(stack->TrackLabelMap());
 }
 
 //_______________________________________________________________________
@@ -504,6 +643,7 @@ void AliMC::FinishEvent()
   //
   
   //
+    
   if(gAlice->Lego()) gAlice->Lego()->FinishEvent();
 
   TIter next(gAlice->Modules());
@@ -533,9 +673,11 @@ void AliMC::FinishEvent()
   header->SetNprimary(stack->GetNprimary());
   header->SetNtrack(stack->GetNtrack());  
 
-  
   // Write out the kinematics
-  stack->FinishEvent();
+  if (!gAlice->Lego()) stack->FinishEvent();
+
+  // Synchronize the TreeTR with TreeK
+  if (fTmpTreeTR) ReorderAndExpandTreeTR();
    
   // Write out the event Header information
   TTree* treeE = runloader->TreeE();
@@ -579,6 +721,10 @@ void AliMC::Init()
 
    //=================Create Materials and geometry
    gMC->Init();
+  // Set alignable volumes for the whole geometry (with old root)
+#if ROOT_VERSION_CODE < 331527
+  SetAllAlignableVolumes();
+#endif
    //Read the cuts for all materials
    ReadTransPar();
    //Build the special IMEDIA table
@@ -587,9 +733,6 @@ void AliMC::Init()
    //Compute cross-sections
    gMC->BuildPhysics();
    
-   //Write Geometry object to current file.
-   gAlice->GetRunLoader()->WriteGeometry();
-
    //Initialise geometry deposition table
    fEventEnergy.Set(gMC->NofVolumes()+1);
    fSummEnergy.Set(gMC->NofVolumes()+1);
@@ -623,6 +766,7 @@ void AliMC::MediaTable()
     if((det=dynamic_cast<AliModule*>(dets[kz]))) {
         TArrayI &idtmed = *(det->GetIdtmed()); 
         for(nz=0;nz<100;nz++) {
+           
        // Find max and min material number
        if((idt=idtmed[nz])) {
          det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
@@ -784,7 +928,6 @@ void AliMC::Browse(TBrowser *b)
   b->Add(fMCQA,"AliMCQA");
 }
 
-//PH
 //_______________________________________________________________________
 void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
 {
@@ -873,7 +1016,7 @@ TParticle* AliMC::Particle(Int_t i) const
 }
 
 //_______________________________________________________________________
-TObjArray* AliMC::Particles() const {
+const TObjArray* AliMC::Particles() const {
   //
   // Returns pointer to Particles array
   //
@@ -961,16 +1104,15 @@ void AliMC::SetCurrentTrack(Int_t track) const
 }
 
 //_______________________________________________________________________
-void  AliMC::AddTrackReference(Int_t label){
+AliTrackReference*  AliMC::AddTrackReference(Int_t label, Int_t id) 
+{
   //
   // add a trackrefernce to the list
-  if (!fTrackReferences) {
-    AliError("Container trackrefernce not active");
-    return;
-  }
-  Int_t nref = fTrackReferences->GetEntriesFast();
-  TClonesArray &lref = *fTrackReferences;
-  new(lref[nref]) AliTrackReference(label);
+  Int_t primary = GetPrimary(label);
+  Particle(primary)->SetBit(kKeepBit);
+
+  Int_t nref = fTmpTrackReferences.GetEntriesFast();
+  return new(fTmpTrackReferences[nref]) AliTrackReference(label, id);
 }
 
 
@@ -981,13 +1123,7 @@ void AliMC::ResetTrackReferences()
   //
   //  Reset all  references
   //
-  if (fTrackReferences)   fTrackReferences->Clear();
-
-  TIter next(gAlice->Modules());
-  AliModule *detector;
-  while((detector = dynamic_cast<AliModule*>(next()))) {
-     detector->ResetTrackReferences();
-  }
+    fTmpTrackReferences.Clear();
 }
 
 void AliMC::RemapTrackReferencesIDs(Int_t *map)
@@ -996,18 +1132,193 @@ void AliMC::RemapTrackReferencesIDs(Int_t *map)
   // Remapping track reference
   // Called at finish primary
   //
-  if (!fTrackReferences) return;
-  for (Int_t i=0;i<fTrackReferences->GetEntries();i++){
-    AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTrackReferences->UncheckedAt(i));
-    if (ref) {
-      Int_t newID = map[ref->GetTrack()];
-      if (newID>=0) ref->SetTrack(newID);
-      else {
-        //ref->SetTrack(-1);
-        ref->SetBit(kNotDeleted,kFALSE);
-        fTrackReferences->RemoveAt(i);  
-      }      
-    }
+    
+  Int_t nEntries = fTmpTrackReferences.GetEntries();
+  for (Int_t i=0; i < nEntries; 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);  
+         }      
+      } // if ref
   }
-  fTrackReferences->Compress();
+  fTmpTrackReferences.Compress();
 }
+
+void AliMC::FixParticleDecaytime()
+{
+    //
+    // Fix the particle decay time according to rmin and rmax for decays
+    //
+
+    TLorentzVector p;
+    gMC->TrackMomentum(p);
+    Double_t tmin, tmax;
+    Double_t b;
+
+    // Transverse velocity 
+    Double_t vt    = p.Pt() / p.E();
+    
+    if ((b = gAlice->Field()->SolenoidField()) > 0.) {     // [kG]
+
+       // Radius of helix
+       
+       Double_t rho   = p.Pt() / 0.0003 / b; // [cm]
+       
+       // Revolution frequency
+       
+       Double_t omega = vt / rho;
+       
+       // Maximum and minimum decay time
+       //
+       // Check for curlers first
+       if (fRDecayMax * fRDecayMax / rho / rho / 2. > 1.) return;
+       
+       //
+       tmax  = TMath::ACos(1. - fRDecayMax * fRDecayMax / rho / rho / 2.) / omega;   // [ct]
+       tmin  = TMath::ACos(1. - fRDecayMin * fRDecayMin / rho / rho / 2.) / omega;   // [ct]
+    } else {
+       tmax =  fRDecayMax / vt;                                                      // [ct] 
+       tmin =  fRDecayMin / vt;                                                      // [ct]
+    }
+    
+    //
+    // Dial t using the two limits
+    Double_t t = tmin + (tmax - tmin) * gRandom->Rndm();                              // [ct]
+    //
+    //
+    // Force decay time in transport code
+    //
+    gMC->ForceDecayTime(t / 2.99792458e10);
+}
+
+void AliMC::MakeTmpTrackRefsTree()
+{
+    // Make the temporary track reference tree
+    fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
+    fTmpTreeTR = new TTree("TreeTR", "Track References");
+    TClonesArray* pRef = &fTmpTrackReferences;
+    fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &pRef, 4000);
+}
+
+void AliMC::ReorderAndExpandTreeTR()
+{
+//
+//  Reorder and expand the temporary track reference tree in order to match the kinematics tree
+//
+
+    AliRunLoader *rl = gAlice->GetRunLoader();
+//
+//  TreeTR
+    AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
+    rl->MakeTrackRefsContainer(); 
+    TTree * treeTR = rl->TreeTR();
+    if (treeTR){
+       // make branch for central track references
+       TBranch *branch;
+       TClonesArray* pRef = &fTrackReferences;
+       branch = treeTR->Branch("TrackReferences", &pRef);
+       branch->SetAddress(&pRef);
+    }
+
+    AliStack* stack  = rl->Stack();
+    Int_t np = stack->GetNprimary();
+    Int_t nt = fTmpTreeTR->GetEntries();
+    //
+    // Loop over tracks and find the secondaries with the help of the kine tree
+    Int_t ifills = 0;
+    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());
+       
+       // Skip primaries that have not been transported
+       Int_t dau1  = part->GetFirstDaughter();
+       Int_t dau2  = -1;
+       if (!part->TestBit(kTransportBit)) continue;
+       //
+       fTmpTreeTR->GetEntry(it++);
+       Int_t nh = fTmpTrackReferences.GetEntries();
+       // Determine range of secondaries produced by this primary
+       if (dau1 > -1) {
+           Int_t inext = ip - 1;
+           while (dau2 < 0) {
+               if (inext >= 0) {
+                   part = stack->Particle(inext);
+                   dau2 =  part->GetFirstDaughter();
+                   if (!(part->TestBit(kTransportBit)) || dau2 == -1 || dau2 < np) {
+//                 if (dau2 == -1 || dau2 < np) {
+                       dau2 = -1;
+                   } else {
+                       dau2--;
+                   }
+               } else {
+                   dau2 = stack->GetNtrack() - 1;
+               }
+               inext--;
+           } // 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);
+               Int_t label = tr->Label();
+               // Skip primaries
+               if (label == ip) continue;
+               if (label > dau2 || label < dau1) 
+                   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();
+                   new(fTrackReferences[nref]) AliTrackReference(*tr);
+               }
+           } // hits
+           treeTR->Fill();
+           fTrackReferences.Clear();
+           ifills++;
+       } // daughters
+    } // tracks
+    //
+    // Now loop again and write the primaries
+    it = nt - 1;
+    for (Int_t ip = 0; ip < np; ip++) {
+       TParticle* part = stack->Particle(ip);
+//     if ((part->GetFirstDaughter() == -1 && part->GetStatusCode() <= 1) || part->GetFirstDaughter() >= np) 
+       if (part->TestBit(kTransportBit))
+       {
+           // Skip particles that have not been transported
+           fTmpTreeTR->GetEntry(it--);
+           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);
+               Int_t label = tr->Label();
+               if (label == ip) {
+                   Int_t nref = fTrackReferences.GetEntriesFast();
+                   new(fTrackReferences[nref]) AliTrackReference(*tr);
+               }
+           } 
+       }
+       treeTR->Fill();
+       fTrackReferences.Clear();
+       ifills++;
+    } // tracks
+    // Check
+    if (ifills != stack->GetNtrack()) 
+       AliWarning(Form("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, stack->GetNtrack()));
+//
+//  Clean-up
+    delete fTmpTreeTR;
+    fTmpFileTR->Close();
+    delete fTmpFileTR;
+    fTmpTrackReferences.Clear();
+    gSystem->Exec("rm -rf TrackRefsTmp.root");
+}
+