]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TFluka/stupre.cxx
Example macros for using AliAODv0 class
[u/mrichter/AliRoot.git] / TFluka / stupre.cxx
index 9f4bc4a48db36114940d60c900c51d47485a1adb..027497f5585350381c855533c1e2263117d79f31 100644 (file)
@@ -1,5 +1,4 @@
 #include <Riostream.h>
-#include "AliRun.h"
 #ifndef WIN32
 # define stupre stupre_
 #else
@@ -8,21 +7,18 @@
 //
 // Fluka include
 #include "Fdimpar.h"  //(DIMPAR) fluka include
+
 // Fluka commons
 #include "Fdblprc.h"  //(DBLPRC) fluka common
 #include "Femfstk.h"  //(EMFSTK) fluka common
 #include "Fevtflg.h"  //(EVTFLG) fluka common
 #include "Fpaprop.h"  //(PAPROP) fluka common
 #include "Ftrackr.h"  //(TRACKR) fluka common
+#include "Femfrgn.h"  //(EFMRGN) fluka common
 
 //Virtual MC
-
-#ifndef WITH_ROOT
 #include "TFluka.h"
-#else
-#include "TFlukaGeo.h"
-#endif
-
+#include "TFlukaCodes.h"
 #include "TVirtualMCStack.h"
 #include "TVirtualMCApplication.h"
 #include "TParticle.h"
@@ -36,39 +32,41 @@ void stupre()
 //*  SeT User PRoperties for Emf particles                               *
 //*                                                                      *
 //*----------------------------------------------------------------------*
-
+// Get the pointer to the VMC
+  TFluka* fluka =  (TFluka*) gMC;
+  static Double_t emassmev = PAPROP.am[kFLUKAelectron + 6] * 1000.;
+    
   Int_t lbhabh = 0;
   if (EVTFLG.ldltry == 1) {
     if (EMFSTK.ichemf[EMFSTK.npemf-1] * EMFSTK.ichemf[EMFSTK.npemf-2] < 0) lbhabh = 1;
   }
 
-// mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
-// mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
+// mkbmx1         = dimension for kwb real spare array in fluka stack in DIMPAR
+// mkbmx2         = dimension for kwb int. spare array in fluka stack in DIMPAR
 // EMFSTK.espark  = spare real variables available for 
 // EMFSTK.iespak  = spare integer variables available for
-// TRACKR.spausr = user defined spare variables for the current particle
-// TRACKR.ispusr = user defined spare flags for the current particle
-// EMFSTK.louemf = user flag
-// TRACKR.llouse = user defined flag for the current particle
+// TRACKR.spausr  = user defined spare variables for the current particle
+// TRACKR.ispusr  = user defined spare flags for the current particle
+// EMFSTK.louemf  = user flag
+// TRACKR.llouse  = user defined flag for the current particle
 
   Int_t npnw, ispr;
-  for (npnw=EMFSTK.npstrt-1; npnw<=EMFSTK.npemf-1; npnw++) {
-
-    for (ispr=0; ispr<=mkbmx1-1; ispr++) 
-      EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
-
-    for (ispr=0; ispr<=mkbmx2-1; ispr++) 
-      EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
-
-    EMFSTK.louemf[npnw] = TRACKR.llouse;
+  for (npnw = EMFSTK.npstrt-1; npnw <= EMFSTK.npemf-1; npnw++) {
+      
+      for (ispr = 0; ispr <= mkbmx1-1; ispr++) 
+         EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
+      
+      for (ispr = 0; ispr <= mkbmx2-1; ispr++) 
+         EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
+    
+      EMFSTK.louemf[npnw] = TRACKR.llouse;
   }
 
-// Get the pointer to the VMC
-  TFluka* fluka =  (TFluka*) gMC;
+
   Int_t verbosityLevel = fluka->GetVerbosityLevel();
   Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
   fluka->SetTrackIsNew(kTRUE);
-//  TVirtualMC* fluka = TFluka::GetMC();
+
 // Get the stack produced from the generator
   TVirtualMCStack* cppstack = fluka->GetStack();
   
@@ -76,31 +74,52 @@ void stupre()
 // Increment the track number and put it into the last flag
 
   Int_t kp;
+  
   for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
-
-//* save the parent track number and reset it at each loop
+    
+// Ckeck transport cut first
+    Int_t    ireg   = EMFSTK.iremf[kp];
+    Double_t cut    = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EMFRGN.elethr[ireg-1] :  EMFRGN.phothr[ireg-1];
+    Double_t e      = EMFSTK.etemf[kp];
+
+    if ((e < cut) 
+       && ( 
+           (EMFSTK.ichemf[kp] ==  0) ||
+           (EMFSTK.ichemf[kp] == -1) ||
+           (EMFSTK.ichemf[kp] ==  1 &&  EMFRGN.phothr[ireg-1] > emassmev)
+           )
+       )
+    {
+       EMFSTK.iespak[kp][mkbmx2-1] = TRACKR.ispusr[mkbmx2-1];
+       EMFSTK.iespak[kp][mkbmx2-2] =  0;
+       EMFSTK.iespak[kp][mkbmx2-3] = TRACKR.jtrack;
+       continue;
+    }
+
+// Save the parent track number and reset it at each loop
     Int_t done = 0;
-
     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
-    
     Int_t flukaid = 0;
+    
+// Identify particle type
+    if      (EMFSTK.ichemf[kp] == -1)  flukaid = kFLUKAelectron;
+    else if (EMFSTK.ichemf[kp] ==  0)  flukaid = kFLUKAphoton;
+    else if (EMFSTK.ichemf[kp] ==  1)  flukaid = kFLUKApositron;
+    
 
-    if (EMFSTK.ichemf[kp] == -1) flukaid = 3;
-    else if (EMFSTK.ichemf[kp] == 0)  flukaid = 7;
-    else if (EMFSTK.ichemf[kp] == 0)  flukaid = 4;
-    Int_t pdg       = fluka->PDGFromId(flukaid);
-    Double_t e      = EMFSTK.etemf[kp] * emvgev;
+    e *= emvgev;
+    Int_t    pdg    = fluka->PDGFromId(flukaid);
     Double_t p      = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
-    Double_t px     = p * EMFSTK.u[kp];
-    Double_t pz     = p * EMFSTK.v[kp];
-    Double_t py     = p * EMFSTK.w[kp];
+    Double_t px     = p * EMFSTK.uemf[kp];
+    Double_t py     = p * EMFSTK.vemf[kp];
+    Double_t pz     = p * EMFSTK.wemf[kp];
     Double_t tof    = EMFSTK.agemf[kp];
     Double_t polx   = EMFSTK.upol[kp];
     Double_t poly   = EMFSTK.vpol[kp];
     Double_t polz   = EMFSTK.wpol[kp];
-    Double_t vx     = EMFSTK.x[kp];
-    Double_t vy     = EMFSTK.y[kp];
-    Double_t vz     = EMFSTK.z[kp];
+    Double_t vx     = EMFSTK.xemf[kp];
+    Double_t vy     = EMFSTK.yemf[kp];
+    Double_t vz     = EMFSTK.zemf[kp];
     Double_t weight = EMFSTK.wtemf[kp];
 
     Int_t ntr;
@@ -118,9 +137,10 @@ void stupre()
         cppstack->PushTrack(done, parent, pdg,
                           px, py, pz, e, vx, vy, vz, tof,
                           polx, poly, polz, mech, ntr, weight, is);
-       if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
+       if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "energy " << e-PAPROP.am[flukaid+6] << endl;
 
        EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+       EMFSTK.iespak[kp][mkbmx2-2] = 0;
     } // end of lpairp, lphoel, lannfl, lannrs
     
 //* Compton: secondary is true only if charged (e+, e-)
@@ -133,6 +153,7 @@ void stupre()
                               polx, poly, polz, mech, ntr, weight, is);
            if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
            EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+           EMFSTK.iespak[kp][mkbmx2-2] = 0;
        }
     } // end of lcmptn
     
@@ -145,6 +166,7 @@ void stupre()
                               polx, poly, polz, mech, ntr, weight, is);
            if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
            EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+           EMFSTK.iespak[kp][mkbmx2-2] = 0;
        }
     } // end of lbrmsp
     
@@ -157,6 +179,7 @@ void stupre()
                                   px, py, pz, e, vx, vy, vz, tof,
                                   polx, poly, polz, mech, ntr, weight, is);
                EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+               EMFSTK.iespak[kp][mkbmx2-2] = 0;
           if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
            } // end of Bhabha
        } // lbhabh == 1
@@ -170,12 +193,11 @@ void stupre()
                               polx, poly, polz, mech, ntr, weight, is);
            if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
            EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+           EMFSTK.iespak[kp][mkbmx2-2] = 0;
        } // end of Delta ray
     } // end of ldltry
     
   } // end of loop
-  
-// !!! TO BE CONFIRMED !!!
 } // end of stupre
 } // end of extern "C"