Possibility to fix the decay time of a primary particle in order to force decays
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Nov 2004 06:16:04 +0000 (06:16 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Nov 2004 06:16:04 +0000 (06:16 +0000)
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.

STEER/AliMC.cxx
STEER/AliMC.h
STEER/libESD.pkg

index defade3..6a8bb1d 100644 (file)
@@ -26,6 +26,8 @@
 #include <TStopwatch.h>
 #include <TSystem.h>
 #include <TVirtualMC.h>
+#include "TGeant3.h"
+
  
 #include "AliLog.h"
 #include "AliDetector.h"
@@ -36,6 +38,7 @@
 #include "AliMCQA.h"
 #include "AliRun.h"
 #include "AliStack.h"
+#include "AliMagF.h"
 #include "AliTrackReference.h"
 
 
@@ -57,6 +60,7 @@ AliMC::AliMC() :
 
 {
   //default constructor
+  DecayLimits();
 }
 
 //_______________________________________________________________________
@@ -77,10 +81,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;
-
 }
 
 //_______________________________________________________________________
@@ -251,10 +254,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())
@@ -623,6 +638,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 +1027,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();
+}
index 238e9b7..0df48b1 100644 (file)
@@ -64,7 +64,10 @@ public:
                                        {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);
@@ -101,8 +104,7 @@ public:
    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;
@@ -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
+   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
index 882b86e..1b8e1dc 100644 (file)
@@ -8,7 +8,7 @@ SRCS = AliESD.cxx \
 
 HDRS:= $(SRCS:.cxx=.h) 
 
-EINCLUDE:=RAW
+EINCLUDE:=$(ALICE)/geant3/TGeant3 RAW
 
 DHDR= ESDLinkDef.h