Updated flags

[u/mrichter/AliRoot.git] / STEER / AliMC.cxx

diff --git a/STEER/AliMC.cxx b/STEER/AliMC.cxx
index defade3d1b2c4e5a31d766a07a847e278ff72c25..16a642ff6fb24f59d6d93427d12ecab0175e5efa 100644 (file)
--- a/STEER/AliMC.cxx
+++ b/STEER/AliMC.cxx
@@ -26,6 +26,9 @@
 #include <TStopwatch.h>
 #include <TSystem.h>
 #include <TVirtualMC.h>
+#include <TGeoManager.h>
+#include "TGeant3.h"
+
  
 #include "AliLog.h"
 #include "AliDetector.h"
@@ -36,6 +39,7 @@
 #include "AliMCQA.h"
 #include "AliRun.h"
 #include "AliStack.h"
+#include "AliMagF.h"
 #include "AliTrackReference.h"
 
 
@@ -49,6 +53,9 @@ AliMC::AliMC() :
   fSum2Energy(0),
   fTrRmax(1.e10),
   fTrZmax(1.e10),
+  fRDecayMax(1.e10),
+  fRDecayMin(0),
+  fDecayPdg(0),
   fImedia(0),
   fTransParName("\0"),
   fMCQA(0),
@@ -57,6 +64,7 @@ AliMC::AliMC() :
 
 {
   //default constructor
+  DecayLimits();
 }
 
 //_______________________________________________________________________
@@ -68,6 +76,9 @@ AliMC::AliMC(const char *name, const char *title) :
   fSum2Energy(0),
   fTrRmax(1.e10),
   fTrZmax(1.e10),
+  fRDecayMax(1.e10),
+  fRDecayMin(0),
+  fDecayPdg(0),
   fImedia(new TArrayI(1000)),
   fTransParName("\0"),
   fMCQA(0),
@@ -77,10 +88,9 @@ AliMC::AliMC(const char *name, const char *title) :
   //constructor
   // Set transport parameters
   SetTransPar();
-
+  DecayLimits();
   // Prepare the tracking medium lists
   for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
-
 }
 
 //_______________________________________________________________________
@@ -92,6 +102,9 @@ AliMC::AliMC(const AliMC &mc) :
   fSum2Energy(0),
   fTrRmax(1.e10),
   fTrZmax(1.e10),
+  fRDecayMax(1.e10),
+  fRDecayMin(0),
+  fDecayPdg(0),
   fImedia(0),
   fTransParName("\0"),
   fMCQA(0),
@@ -99,7 +112,7 @@ AliMC::AliMC(const AliMC &mc) :
   fTrackReferences(0)
 {
   //
-  // Copy constructor for AliRun
+  // Copy constructor for AliMC
   //
   mc.Copy(*this);
 }
@@ -132,9 +145,24 @@ 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
+    const char *geomfilename = gAlice->GetGeometryFileName();
+    if(gSystem->ExpandPathName(geomfilename)){
+      AliInfo(Form("Loading geometry from file:\n %40s\n\n",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,6 +175,8 @@ void  AliMC::ConstructGeometry()
       AliInfo(Form("%10s R:%.2fs C:%.2fs",
                   detector->GetName(),stw.RealTime(),stw.CpuTime()));
     }
+  }
+  
 }
 
 //_______________________________________________________________________
@@ -156,19 +186,20 @@ void  AliMC::InitGeometry()
   // Initialize detectors and display geometry
   //
 
-   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 and display geometry
+    if(gAlice->IsRootGeometry()) detector->CreateMaterials();
+    detector->Init();
+    detector->BuildGeometry();
+    AliInfo(Form("%10s R:%.2fs C:%.2fs",
+                detector->GetName(),stw.RealTime(),stw.CpuTime()));
+  }
+  
 }
 
 //_______________________________________________________________________
@@ -251,10 +282,22 @@ void AliMC::Stepping()
   //
   // Called at every step during transport
   //
-
+    
   Int_t id = DetFromMate(gMC->GetMedium());
   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())
@@ -456,12 +499,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();
 }
 
 //_______________________________________________________________________
@@ -623,6 +666,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;
@@ -1011,3 +1055,52 @@ void AliMC::RemapTrackReferencesIDs(Int_t *map)
   }
   fTrackReferences->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
+    //
+    TGeant3 * geant = (TGeant3*) gMC;
+    geant->Gcphys()->sumlif = t / p.Beta() / p.Gamma();
+}