Possibility to fix the decay time of a primary particle in order to force decays

within a user chosen minumum and maximum radius.
The current implementation is TGeant3 dependent. This is a temporary solution
to be replaced by corresponding TVirtualMC method calls.

diff --git a/STEER/AliMC.cxx b/STEER/AliMC.cxx
index defade3d1b2c4e5a31d766a07a847e278ff72c25..6a8bb1d484b1cf66455fbac494fe81a93522d953 100644 (file)
--- a/STEER/AliMC.cxx
+++ b/STEER/AliMC.cxx
@@ -26,6 +26,8 @@
 #include <TStopwatch.h>
 #include <TSystem.h>
 #include <TVirtualMC.h>
 #include <TStopwatch.h>
 #include <TSystem.h>
 #include <TVirtualMC.h>

+#include "TGeant3.h"
+

  
 #include "AliLog.h"
 #include "AliDetector.h"
  
 #include "AliLog.h"
 #include "AliDetector.h"


@@ -36,6 +38,7 @@
 #include "AliMCQA.h"
 #include "AliRun.h"
 #include "AliStack.h"
 #include "AliMCQA.h"
 #include "AliRun.h"
 #include "AliStack.h"

+#include "AliMagF.h"

 #include "AliTrackReference.h"
 
 
 #include "AliTrackReference.h"
 
 


@@ -57,6 +60,7 @@ AliMC::AliMC() :
 
 {
   //default constructor
 
 {
   //default constructor

+  DecayLimits();

 }
 
 //_______________________________________________________________________
 }
 
 //_______________________________________________________________________


@@ -77,10 +81,9 @@ AliMC::AliMC(const char *name, const char *title) :
   //constructor
   // Set transport parameters
   SetTransPar();
   //constructor
   // Set transport parameters
   SetTransPar();

-
+  DecayLimits();

   // Prepare the tracking medium lists
   for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
   // Prepare the tracking medium lists
   for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;

-

 }
 
 //_______________________________________________________________________
 }
 
 //_______________________________________________________________________


@@ -251,10 +254,22 @@ void AliMC::Stepping()
   //
   // Called at every step during transport
   //
   //
   // Called at every step during transport
   //

-
+    

   Int_t id = DetFromMate(gMC->GetMedium());
   if (id < 0) return;
 
   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())
   //
   // --- If lego option, do it and leave 
   if (gAlice->Lego())


@@ -623,6 +638,7 @@ void AliMC::MediaTable()
     if((det=dynamic_cast<AliModule*>(dets[kz]))) {
         TArrayI &idtmed = *(det->GetIdtmed()); 
         for(nz=0;nz<100;nz++) {
     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;
        // Find max and min material number
        if((idt=idtmed[nz])) {
          det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;


@@ -1011,3 +1027,52 @@ void AliMC::RemapTrackReferencesIDs(Int_t *map)
   }
   fTrackReferences->Compress();
 }
   }
   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();
+}

diff --git a/STEER/AliMC.h b/STEER/AliMC.h
index 238e9b76c00945222a5f23447a1352527573b510..0df48b1455eed3dee6e1ff22c2264605ad85d12f 100644 (file)
--- a/STEER/AliMC.h
+++ b/STEER/AliMC.h
@@ -64,7 +64,10 @@ public:
                                        {fEventEnergy[id]+=edep;}
    virtual  void  ResetHits();
    virtual  void  TrackingLimits( Float_t rmax=1.e10, Float_t zmax=1.e10)
                                        {fEventEnergy[id]+=edep;}
    virtual  void  ResetHits();
    virtual  void  TrackingLimits( Float_t rmax=1.e10, Float_t zmax=1.e10)

-    {fTrRmax=rmax; fTrZmax=zmax;}
+       {fTrRmax=rmax; fTrZmax=zmax;}
+   virtual  void  DecayLimits( Float_t rmin = -1., Float_t rmax = -1., Int_t pdg = 0)
+       {fRDecayMin = rmin; fRDecayMax = rmax; fDecayPdg = pdg;}
+   

    virtual  void  Init();
    virtual  void  SetTransPar(const char *filename="$(ALICE_ROOT)/data/galice.cuts");
    virtual  void  Browse(TBrowser *b);
    virtual  void  Init();
    virtual  void  SetTransPar(const char *filename="$(ALICE_ROOT)/data/galice.cuts");
    virtual  void  Browse(TBrowser *b);


@@ -101,8 +104,7 @@ public:
    TClonesArray   *TrackReferences()   const {return fTrackReferences;}
    virtual void   RemapTrackReferencesIDs(Int_t *map); //remaping track references MI
    virtual void   ResetTrackReferences();
    TClonesArray   *TrackReferences()   const {return fTrackReferences;}
    virtual void   RemapTrackReferencesIDs(Int_t *map); //remaping track references MI
    virtual void   ResetTrackReferences();

-   
-
+   virtual void   FixParticleDecaytime();

 
 private:
    void Copy (TObject &mc) const;
 
 private:
    void Copy (TObject &mc) const;


@@ -112,6 +114,9 @@ private:
    TArrayF        fSum2Energy;        //! Energy squared per event in each volume
    Float_t        fTrRmax;            //  Maximum radius for tracking
    Float_t        fTrZmax;            //  Maximu z for tracking
    TArrayF        fSum2Energy;        //! Energy squared per event in each volume
    Float_t        fTrRmax;            //  Maximum radius for tracking
    Float_t        fTrZmax;            //  Maximu z for tracking

+   Float_t        fRDecayMax;         //  Maximum radius for decay
+   Float_t        fRDecayMin;         //  Minimum radius for decay
+   Int_t          fDecayPdg;          //  PDG code of particle with forced decay length

    TArrayI       *fImedia;            //! Array of correspondence between media and detectors
    TString        fTransParName;      //  Name of the transport parameters file
    AliMCQA       *fMCQA;              //  Pointer to MC Quality assurance class
    TArrayI       *fImedia;            //! Array of correspondence between media and detectors
    TString        fTransParName;      //  Name of the transport parameters file
    AliMCQA       *fMCQA;              //  Pointer to MC Quality assurance class

diff --git a/STEER/libESD.pkg b/STEER/libESD.pkg
index 882b86e4c5ec538403947961ac9f37801bfb0305..1b8e1dc2d5059fa4a63506a2407aa3c0dcb4a6f8 100644 (file)
--- a/STEER/libESD.pkg
+++ b/STEER/libESD.pkg
@@ -8,7 +8,7 @@ SRCS = AliESD.cxx \
 
 HDRS:= $(SRCS:.cxx=.h) 
 
 
 HDRS:= $(SRCS:.cxx=.h) 
 

-EINCLUDE:=RAW
+EINCLUDE:=$(ALICE)/geant3/TGeant3 RAW

 
 DHDR= ESDLinkDef.h
 
 
 DHDR= ESDLinkDef.h