In decays with AliDecayer put long-lived particles undecayed on the stack.
[u/mrichter/AliRoot.git] / TGeant3 / AliGeant3.cxx
index 3f02c2f..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.
 
@@ -59,6 +74,7 @@ ReadEuclid moved from AliRun to AliModule
 #include <stdlib.h>
 
 #include <TParticle.h>
+#include <TStopwatch.h>
 
 #include "AliDecayer.h"
 #include "AliGeant3.h"
@@ -73,10 +89,12 @@ ReadEuclid moved from AliRun to AliModule
 
 # define rxgtrak rxgtrak_
 # define rxouth  rxouth_
+# define rxinh   rxinh_
 #else
 
 # define rxgtrak RXGTRAK 
 # define rxouth  RXOUTH
+# define rxinh   RXINH
 #endif
 
 ClassImp(AliGeant3)
@@ -98,22 +116,36 @@ 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->Init();
-    detector->BuildGeometry();
-  }
+    //
+    //=================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();
 }
 
 //____________________________________________________________________________
@@ -214,6 +246,15 @@ extern "C" void type_of_call  rxouth ()
   gAlice->FinishPrimary();
 }
 
+//_____________________________________________________________________________
+extern "C" void type_of_call  rxinh ()
+{
+  //
+  // Called by Gtreve at the beginning of each primary track
+  //
+  gAlice->BeginPrimary();
+}
+
 #ifndef WIN32
 #  define gudigi gudigi_
 #  define guhadr guhadr_
@@ -396,8 +437,8 @@ void gudcay()
 //
     
     TGeant3* geant3=(TGeant3*) gMC;
-  // set decay table
-  gMC->Decayer()->ForceDecay();
+    // set decay table
+    gMC->Decayer()->ForceDecay();
 
 // Initialize 4-momentum vector    
     Int_t ipart = geant3->Gckine()->ipart;
@@ -419,15 +460,41 @@ void gudcay()
 // Fetch Particles
     Int_t np = geant3->Decayer()->ImportParticles(particles);
     if (np <=1) return;
-    
-    for (Int_t i=0; i< np; i++) 
+
+    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++) 
     {
-       TParticle *  iparticle = (TParticle *) particles->At(i);
+       iparticle = (TParticle *) particles->At(i);
+       ipF = iparticle->GetFirstDaughter();
+       ipL = iparticle->GetLastDaughter();     
+       Int_t kf = iparticle->GetPdgCode();
        Int_t ks = iparticle->GetStatusCode();
-// Final state particles only
-       if (ks < 1 || ks >  10) continue;
+//
+// 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
-       Int_t  kf=iparticle->GetPdgCode();
        if (kf==12 || kf ==-12) continue;
        if (kf==14 || kf ==-14) continue;
        if (kf==16 || kf ==-16) continue;
@@ -435,7 +502,7 @@ void gudcay()
        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();
@@ -450,10 +517,11 @@ void gudcay()
        (geant3->Gckin3()->gpos[index][2]) = geant3->Gctrak()->vect[2];
 // time of flight offset (mm)
        (geant3->Gcking()->tofd[index])    = 0.;
-//     (geant3->Gcking()->tofd[index])    = iparticle->T()/(10*kSpeedOfLight);
 // increase stack counter
        (geant3->Gcking()->ngkine)=index+1;
     }
+//
+    if (pFlag) delete[] pFlag;
 }
 
 //______________________________________________________________________
@@ -744,7 +812,7 @@ void gustep()
 
   // --- Add new created particles 
   if (gMC->NSecondaries() > 0) {
-    pProc=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!