TFluka/stuprf.cxx

   1 #include <Riostream.h>
   2 #include "AliRun.h"
   3 #include "AliStack.h"
   4 #ifndef WIN32
   5 # define stuprf stuprf_
   6 #else
   7 # define stuprf STUPRF
   8 #endif
   9 //
  10 // Fluka include
  11 #include "Fdimpar.h"  //(DIMPAR) fluka include
  12 // Fluka commons
  13 #include "Fdblprc.h"  //(DBLPRC) fluka common
  14 #include "Fevtflg.h"  //(EVTFLG) fluka common
  15 #include "Fpaprop.h"  //(PAPROP) fluka common
  16 #include "Fstack.h"   //(STACK)  fluka common
  17 #include "Ftrackr.h"  //(TRACKR) fluka common
  18 #include "Ffinuc.h"   //(FINUC)  fluka common
  19
  20 //Virtual MC
  21 #include "TFluka.h"
  22
  23 #include "TVirtualMCStack.h"
  24 #include "TVirtualMCApplication.h"
  25 #include "TParticle.h"
  26 #include "TVector3.h"
  27
  28 extern "C" {
  29     void stuprf(Int_t& /*ij*/, Int_t& /*mreg*/,
  30                 Double_t& xx, Double_t& yy, Double_t& zz,
  31                 Int_t& numsec, Int_t& npprmr)
  32 {
  33 //*----------------------------------------------------------------------*
  34 //*                                                                      *
  35 //*  SeT User PRoperties for Fluka particles                             *
  36 //*                                                                      *
  37 //*----------------------------------------------------------------------*
  38
  39 // STACK.lstack  = stack pointer
  40 // STACK.louse   = user flag
  41 // TRACKR.llouse = user defined flag for the current particle
  42   STACK.louse[STACK.lstack] = TRACKR.llouse;
  43
  44 // mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
  45 // mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
  46 // STACK.sparek  = spare real variables available for k.w.burn
  47 // STACK.ispark  = spare integer variables available for k.w.burn
  48 // TRACKR.spausr = user defined spare variables for the current particle
  49 // TRACKR.ispusr = user defined spare flags for the current particle
  50   Int_t ispr;
  51   for (ispr = 0; ispr <= mkbmx1 - 1; ispr++) {
  52     STACK.sparek[STACK.lstack][ispr] = TRACKR.spausr[ispr];
  53   }
  54   for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
  55     STACK.ispark[STACK.lstack][ispr] = TRACKR.ispusr[ispr];
  56   }
  57
  58 // Get the pointer to the VMC
  59   TFluka* fluka =  (TFluka*) gMC;
  60   Int_t verbosityLevel = fluka->GetVerbosityLevel();
  61   Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
  62   fluka->SetTrackIsNew(kTRUE);
  63 //  TVirtualMC* fluka = TFluka::GetMC();
  64 // Get the stack produced from the generator
  65   TVirtualMCStack* cppstack = fluka->GetStack();
  66
  67 // EVTFLG.ntrcks = track number
  68 // Increment the track number and put it into the last flag
  69 // was numsec -1
  70 // clarify with Alberto
  71   if (numsec > npprmr) {
  72 // Now call the PushTrack(...)
  73     Int_t done = 0;
  74
  75     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
  76     Int_t kpart  = FINUC.kpart[numsec-1];
  77     if (kpart < -6) return;
  78
  79     Int_t pdg = fluka->PDGFromId(kpart);
  80
  81
  82     Double_t px = FINUC.plr[numsec-1] * FINUC.cxr[numsec-1];
  83     Double_t py = FINUC.plr[numsec-1] * FINUC.cyr[numsec-1];
  84     Double_t pz = FINUC.plr[numsec-1] * FINUC.czr[numsec-1];
  85     Double_t e  = FINUC.tki[numsec-1] + PAPROP.am[FINUC.kpart[numsec-1]+6];
  86
  87     Double_t vx = xx;
  88     Double_t vy = yy;
  89     Double_t vz = zz;
  90
  91     Double_t tof  = TRACKR.atrack;
  92     Double_t polx = FINUC.cxrpol[numsec-1];
  93     Double_t poly = FINUC.cyrpol[numsec-1];
  94     Double_t polz = FINUC.czrpol[numsec-1];
  95
  96
  97     TMCProcess mech = kPHadronic;
  98
  99     if (EVTFLG.ldecay == 1) {
 100         mech = kPDecay;
 101         if (debug) cout << endl << "Decay" << endl;
 102
 103     } else if (EVTFLG.ldltry == 1) {
 104         mech = kPDeltaRay;
 105         if (debug) cout << endl << "Delta Ray" << endl;
 106
 107     } else if (EVTFLG.lpairp == 1) {
 108         mech = kPPair;
 109         if (debug) cout << endl << "Pair Production" << endl;
 110
 111     } else if (EVTFLG.lbrmsp == 1) {
 112         mech = kPBrem;
 113         if (debug) cout << endl << "Bremsstrahlung" << endl;
 114
 115     }
 116
 117
 118     Double_t weight = FINUC.wei[numsec-1];
 119     Int_t is = 0;
 120     Int_t ntr;
 121     //
 122     // Save particle in VMC stack
 123     cppstack->PushTrack(done, parent, pdg,
 124                        px, py, pz, e,
 125                        vx, vy, vz, tof,
 126                        polx, poly, polz,
 127                        mech, ntr, weight, is);
 128     if (debug) cout << endl << " !!! stuprf: ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "numsec "
 129          << numsec << "npprmr " << npprmr << endl;
 130 //
 131 //  Save current track number
 132     STACK.ispark[STACK.lstack][mkbmx2-1] = ntr;
 133     STACK.ispark[STACK.lstack][mkbmx2-2] = 0;
 134   } // end of if (numsec-1 > npprmr)
 135 } // end of stuprf
 136 } // end of extern "C"
 137