In decays with AliDecayer put long-lived particles undecayed on the stack.
[u/mrichter/AliRoot.git] / TGeant3 / AliGeant3.cxx
index 0128691..2c8bc63 100644 (file)
 
 /*
 $Log$
+Revision 1.17  2001/06/15 09:31:23  morsch
+In gudcay: write only first generation decay products to stack to respect the possibility of secondary, tertiary, ... vertices during tracking.
+
+Revision 1.16  2001/05/16 14:57:23  alibrary
+New files for folders and Stack
+
+Revision 1.15  2001/03/20 06:28:49  alibrary
+New detector loop split in 2
+
+Revision 1.14  2000/12/20 08:39:39  fca
+Support for Cerenkov and process list in Virtual MC
+
+Revision 1.13  2000/11/30 07:12:54  alibrary
+Introducing new Rndm and QA classes
+
+Revision 1.12  2000/11/06 11:35:46  morsch
+Call BuildGeometry() after Init() to be able to share common detector parameters.
+
+Revision 1.11  2000/10/04 16:30:22  fca
+Add include for exit()
+
+Revision 1.10  2000/10/02 21:28:16  fca
+Removal of useless dependecies via forward declarations
+
+Revision 1.9  2000/09/12 14:36:17  morsch
+In gudcay(): call ForceDecay() before Decay()
+
+Revision 1.8  2000/09/06 14:56:34  morsch
+gudcay() method implemented.
+Decays are performed using  the AliDecayer interface. The pointer to the instance of AliDecayer
+is obtained from geant3 (will be gMC once it has been implemented for Geant4).
+
+Revision 1.7  2000/07/13 16:19:10  fca
+Mainly coding conventions + some small bug fixes
+
+Revision 1.6  2000/07/11 18:24:59  fca
+Coding convention corrections + few minor bug fixes
+
+Revision 1.5  2000/05/20 14:49:48  fca
+Call gdebug at the end of gustep
+
 Revision 1.4  2000/04/26 10:17:32  fca
 Changes in Lego for G4 compatibility
 
@@ -30,25 +71,30 @@ ReadEuclid moved from AliRun to AliModule
 
 */
 
+#include <stdlib.h>
+
 #include <TParticle.h>
+#include <TStopwatch.h>
 
+#include "AliDecayer.h"
 #include "AliGeant3.h"
 #include "AliRun.h"
 #include "TGeant3.h"
 #include "AliCallf77.h" 
+#include "AliModule.h"
+#include "AliMagF.h"
+#include "AliGenerator.h"
 
 #ifndef WIN32
 
 # define rxgtrak rxgtrak_
-# define rxstrak rxstrak_
-# define rxkeep  rxkeep_ 
 # define rxouth  rxouth_
+# define rxinh   rxinh_
 #else
 
 # define rxgtrak RXGTRAK 
-# define rxstrak RXSTRAK 
-# define rxkeep  RXKEEP  
 # define rxouth  RXOUTH
+# define rxinh   RXINH
 #endif
 
 ClassImp(AliGeant3)
@@ -56,8 +102,12 @@ ClassImp(AliGeant3)
 AliGeant3::AliGeant3(const char *title) : 
   TGeant3(title) {}
 
+//____________________________________________________________________________
 void AliGeant3::FinishGeometry()
 {
+  //
+  // Finalise geometry construction
+  //
   TGeant3::FinishGeometry();
   //Create the color table
   SetColors();
@@ -66,26 +116,44 @@ void AliGeant3::FinishGeometry()
 //____________________________________________________________________________
 void AliGeant3::Init()
 {
-  //
-  //=================Create Materials and geometry
-  TObjArray *modules = gAlice->Modules();
-  TIter next(modules);
-  AliModule *detector;
-  while((detector = (AliModule*)next())) {
-    // Initialise detector materials and geometry
-    detector->CreateMaterials();
-    detector->CreateGeometry();
-    detector->BuildGeometry();
-    detector->Init();
-  }
+    //
+    //=================Create Materials and geometry
+    //
+    TStopwatch stw;
+    TObjArray *modules = gAlice->Modules();
+    TIter next(modules);
+    AliModule *detector;
+    printf("Geometry creation:\n");
+    while((detector = (AliModule*)next())) {
+      stw.Start();
+      // Initialise detector materials and geometry
+      detector->CreateMaterials();
+      detector->CreateGeometry();
+      printf("%10s R:%.2fs C:%.2fs\n",
+            detector->GetName(),stw.RealTime(),stw.CpuTime());
+    }
+    //Terminate building of geometry
+    FinishGeometry();
+    
+    printf("Initialisation:\n");
+    next.Reset();
+    while((detector = (AliModule*)next())) {
+      stw.Start();
+      // Initialise detector and display geometry
+      detector->Init();
+      detector->BuildGeometry();
+      printf("%10s R:%.2fs C:%.2fs\n",
+            detector->GetName(),stw.RealTime(),stw.CpuTime());
+    }
 
-  //Terminate building of geometry
-  FinishGeometry();
 }
 
 //____________________________________________________________________________
 void AliGeant3::ProcessRun(Int_t nevent)
 {
+  //
+  // Process the run
+  //
   Int_t todo = TMath::Abs(nevent);
   for (Int_t i=0; i<todo; i++) {
   // Process one run (one run = one event)
@@ -95,8 +163,12 @@ void AliGeant3::ProcessRun(Int_t nevent)
   }
 }
 
+//_____________________________________________________________________________
 void AliGeant3::ProcessEvent()
 {
+  //
+  // Process one event
+  //
   Gtrigi();
   Gtrigc();
   Gtrig();
@@ -166,70 +238,23 @@ extern "C" void type_of_call  rxgtrak (Int_t &mtrack, Int_t &ipart, Float_t *pmo
 }
 
 //_____________________________________________________________________________
-extern "C" void type_of_call 
-#ifndef WIN32
-rxstrak (Int_t &keep, Int_t &parent, Int_t &ipart, Float_t *pmom, 
-              Float_t *vpos, Float_t &tof, const char* cmech, Int_t &ntr, const int cmlen)
-#else
-rxstrak (Int_t &keep, Int_t &parent, Int_t &ipart, Float_t *pmom,
-        Float_t *vpos, Float_t &tof, const char* cmech, const int cmlen,
-        Int_t &ntr)
-#endif
+extern "C" void type_of_call  rxouth ()
 {
   //
-  //     Fetches next track from the ROOT stack for transport. Called by GUKINE
-  //     and GUSTEP.
-  //
-  //              Status of the track. If keep=0 the track is put
-  //      keep    on the ROOT stack but it is not fetched for
-  //              transport.
-  //      parent  Parent track. If parent=0 the track is a primary.
-  //              In GUSTEP the routine is normally called to store
-  //              secondaries generated by the current track whose
-  //              ROOT stack number is MTRACK (common SCKINE.
-  //      ipart   Particle code in the GEANT conventions.
-  //      pmom[3] Particle momentum in GeV/c
-  //      vpos[3] Particle position
-  //      tof     Particle time of flight in seconds
-  //
-  //      cmech   (CHARACTER*10) Particle origin. This field is user
-  //              defined and it is not used inside the GALICE code.
-  //      ntr     Number assigned to the particle in the ROOT stack.
+  // Called by Gtreve at the end of each primary track
   //
-  char mecha[11];
-  Float_t polar[3]={0.,0.,0.};
-  for(int i=0; i<10 && i<cmlen; i++) mecha[i]=cmech[i];
-  mecha[10]=0;
-  Int_t pdg=gMC->PDGFromId(ipart);
-  gAlice->SetTrack(keep, parent-1, pdg, pmom, vpos, polar, tof, mecha, ntr);
-  ntr++;
-}
-
-//_____________________________________________________________________________
-extern "C" void type_of_call  rxkeep(const Int_t &n)
-{
-  if( NULL==gAlice ) exit(1);
-  
-  if( n<=0 || n>gAlice->Particles()->GetEntries() )
-    {
-      printf("  Bad index n=%d must be 0<n<=%d\n",
-            n,gAlice->Particles()->GetEntries());
-      exit(1);
-    }
-  
-  ((TParticle*)(gAlice->Particles()->UncheckedAt(n-1)))->SetBit(Keep_Bit);
+  gAlice->FinishPrimary();
 }
 
 //_____________________________________________________________________________
-extern "C" void type_of_call  rxouth ()
+extern "C" void type_of_call  rxinh ()
 {
   //
-  // Called by Gtreve at the end of each primary track
+  // Called by Gtreve at the beginning of each primary track
   //
-  gAlice->FinishPrimary();
+  gAlice->BeginPrimary();
 }
 
-
 #ifndef WIN32
 #  define gudigi gudigi_
 #  define guhadr guhadr_
@@ -263,7 +288,7 @@ extern "C" void type_of_call  rxouth ()
 #  define ghelix ghelix_
 #  define grkuta grkuta_
 #  define gtrack gtrack_
-#  define gtreve_root gtreve_root_
+#  define gtreveroot gtreveroot_
 #  define glast  glast_
 
 #else
@@ -299,7 +324,7 @@ extern "C" void type_of_call  rxouth ()
 #  define ghelix GHELIX
 #  define grkuta GRKUTA
 #  define gtrack GTRACK
-#  define gtreve_root GTREVE_ROOT
+#  define gtreveroot GTREVEROOT
 #  define glast  GLAST
 
 #endif
@@ -314,7 +339,7 @@ extern "C" type_of_call void ghelx3(Float_t&, Float_t&, Float_t*, Float_t*);
 extern "C" type_of_call void ghelix(Float_t&, Float_t&, Float_t*, Float_t*);
 extern "C" type_of_call void grkuta(Float_t&, Float_t&, Float_t*, Float_t*);
 extern "C" type_of_call void gtrack();
-extern "C" type_of_call void gtreve_root();
+extern "C" type_of_call void gtreveroot();
 extern "C" type_of_call void glast();
 
 extern "C" type_of_call {
@@ -410,6 +435,93 @@ void gudcay()
 //
 //    ------------------------------------------------------------------
 //
+    
+    TGeant3* geant3=(TGeant3*) gMC;
+    // set decay table
+    gMC->Decayer()->ForceDecay();
+
+// Initialize 4-momentum vector    
+    Int_t ipart = geant3->Gckine()->ipart;
+    TLorentzVector p;
+    
+    p[0]=geant3->Gctrak()->vect[3];
+    p[1]=geant3->Gctrak()->vect[4];
+    p[2]=geant3->Gctrak()->vect[5];
+    p[3]=geant3->Gctrak()->vect[6];    
+    
+// Convert from geant to lund particle code
+    Int_t iplund=gMC->PDGFromId(ipart);
+// Particle list
+    static TClonesArray *particles;
+    if(!particles) particles=new TClonesArray("TParticle",1000);
+// Decay
+    gMC->Decayer()->Decay(iplund, &p);
+    
+// Fetch Particles
+    Int_t np = geant3->Decayer()->ImportParticles(particles);
+    if (np <=1) return;
+
+    TParticle *  iparticle = (TParticle *) particles->At(0);
+    Int_t ipF = 0, ipL = 0 ;
+    Int_t i,j;
+
+// Array to flag deselected particles
+    Int_t* pFlag = new Int_t[np];
+    for (i=0; i<np; i++) pFlag[i]=0;
+// Particle loop
+    for (i=1; i < np; i++) 
+    {
+       iparticle = (TParticle *) particles->At(i);
+       ipF = iparticle->GetFirstDaughter();
+       ipL = iparticle->GetLastDaughter();     
+       Int_t kf = iparticle->GetPdgCode();
+       Int_t ks = iparticle->GetStatusCode();
+//
+// Deselect daughters of deselected particles
+// and jump skip the current particle
+       if (pFlag[i] == 1) {
+           if (ipF >= 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
+           continue;
+       } // deselected ??
+// Particles with long life-time are put on the stack for further tracking
+// Decay products are deselected
+//     
+       if (ks != 1) { 
+           Double_t lifeTime = gMC->Decayer()->GetLifetime(kf);
+           if (lifeTime > (Double_t) 1.e-15) {
+               if (ipF >= 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
+           } else{
+               continue;
+           }
+       } // ks==1 ?
+// Skip neutrinos
+       if (kf==12 || kf ==-12) continue;
+       if (kf==14 || kf ==-14) continue;
+       if (kf==16 || kf ==-16) continue;
+       
+       Int_t index=geant3->Gcking()->ngkine;
+// Put particle on geant stack
+// momentum vector
+       
+       (geant3->Gcking()->gkin[index][0]) = iparticle->Px();
+       (geant3->Gcking()->gkin[index][1]) = iparticle->Py();
+       (geant3->Gcking()->gkin[index][2]) = iparticle->Pz();
+       (geant3->Gcking()->gkin[index][3]) = iparticle->Energy();
+       Int_t ilu = gMC->IdFromPDG(kf);
+
+// particle type       
+       (geant3->Gcking()->gkin[index][4]) = Float_t(ilu);
+// position
+       (geant3->Gckin3()->gpos[index][0]) = geant3->Gctrak()->vect[0];
+       (geant3->Gckin3()->gpos[index][1]) = geant3->Gctrak()->vect[1];
+       (geant3->Gckin3()->gpos[index][2]) = geant3->Gctrak()->vect[2];
+// time of flight offset (mm)
+       (geant3->Gcking()->tofd[index])    = 0.;
+// increase stack counter
+       (geant3->Gcking()->ngkine)=index+1;
+    }
+//
+    if (pFlag) delete[] pFlag;
 }
 
 //______________________________________________________________________
@@ -628,20 +740,11 @@ void gutrak()
 //
 //    ------------------------------------------------------------------
 //
-     Int_t ndet = gAlice->Modules()->GetLast();
-     TObjArray &dets = *gAlice->Modules();
-     AliModule *module;
-     Int_t i;
-
-     for(i=0; i<=ndet; i++)
-       if((module = (AliModule*)dets[i]))
-        module->PreTrack();
+     gAlice->PreTrack();
 
      gtrack();
 
-     for(i=0; i<=ndet; i++)
-       if((module = (AliModule*)dets[i]))
-        module->PostTrack();
+     gAlice->PostTrack();
 }
 
 //______________________________________________________________________
@@ -659,7 +762,7 @@ void gutrev()
 //
 //    ------------------------------------------------------------------
 //
-  gtreve_root();
+  gtreveroot();
 }
 
 
@@ -696,10 +799,10 @@ void gustep()
   Int_t ipp, jk, id, nt;
   Float_t polar[3]={0,0,0};
   Float_t mom[3];
-  const char *chproc;
+  AliMCProcess pProc;
+
   
   TGeant3* geant3 = (TGeant3*) gMC;
-
   //     Stop particle if outside user defined tracking region 
   gMC->TrackPosition(x);
   r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
@@ -709,13 +812,13 @@ void gustep()
 
   // --- Add new created particles 
   if (gMC->NSecondaries() > 0) {
-    chproc=gMC->ProdProcess();
+    pProc=gMC->ProdProcess(0);
     for (jk = 0; jk < geant3->Gcking()->ngkine; ++jk) {
       ipp = Int_t (geant3->Gcking()->gkin[jk][4]+0.5);
       // --- Skip neutrinos! 
       if (ipp != 4) {
        gAlice->SetTrack(1,gAlice->CurrentTrack(),gMC->PDGFromId(ipp), geant3->Gcking()->gkin[jk], 
-                        geant3->Gckin3()->gpos[jk], polar,geant3->Gctrak()->tofg, chproc, nt);
+                        geant3->Gckin3()->gpos[jk], polar,geant3->Gctrak()->tofg, pProc, nt);
       }
     }
   }
@@ -730,7 +833,7 @@ void gustep()
                       geant3->Gckin2()->xphot[jk],     //position
                       &geant3->Gckin2()->xphot[jk][7], //polarisation
                       geant3->Gckin2()->xphot[jk][10], //time of flight
-                      "Cherenkov", nt);
+                      kPCerenkov, nt);
       }
   }
   // --- Particle leaving the setup ?