TFluka/stupre.cxx

   1 #include <Riostream.h>
   2 #ifndef WIN32
   3 # define stupre stupre_
   4 #else
   5 # define stupre STUPRE
   6 #endif
   7 //
   8 // Fluka include
   9 #include "Fdimpar.h"  //(DIMPAR) fluka include
  10
  11 // Fluka commons
  12 #include "Fdblprc.h"  //(DBLPRC) fluka common
  13 #include "Femfstk.h"  //(EMFSTK) fluka common
  14 #include "Fevtflg.h"  //(EVTFLG) fluka common
  15 #include "Fpaprop.h"  //(PAPROP) fluka common
  16 #include "Ftrackr.h"  //(TRACKR) fluka common
  17 #include "Femfrgn.h"  //(EFMRGN) fluka common
  18
  19 //Virtual MC
  20 #include "TFluka.h"
  21 #include "TFlukaCodes.h"
  22 #include "TVirtualMCStack.h"
  23 #include "TVirtualMCApplication.h"
  24 #include "TParticle.h"
  25 #include "TVector3.h"
  26
  27 extern "C" {
  28 void stupre()
  29 {
  30 //*----------------------------------------------------------------------*
  31 //*                                                                      *
  32 //*  SeT User PRoperties for Emf particles                               *
  33 //*                                                                      *
  34 //*----------------------------------------------------------------------*
  35 // Get the pointer to the VMC
  36   TFluka* fluka =  (TFluka*) gMC;
  37   static Double_t emassmev = PAPROP.am[kFLUKAelectron + 6] * 1000.;
  38
  39   Int_t lbhabh = 0;
  40   if (EVTFLG.ldltry == 1) {
  41     if (EMFSTK.ichemf[EMFSTK.npemf-1] * EMFSTK.ichemf[EMFSTK.npemf-2] < 0) lbhabh = 1;
  42   }
  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 // EMFSTK.espark  = spare real variables available for
  47 // EMFSTK.iespak  = spare integer variables available for
  48 // TRACKR.spausr  = user defined spare variables for the current particle
  49 // TRACKR.ispusr  = user defined spare flags for the current particle
  50 // EMFSTK.louemf  = user flag
  51 // TRACKR.llouse  = user defined flag for the current particle
  52
  53   Int_t npnw, ispr;
  54   for (npnw = EMFSTK.npstrt-1; npnw <= EMFSTK.npemf-1; npnw++) {
  55
  56       for (ispr = 0; ispr <= mkbmx1-1; ispr++)
  57           EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
  58
  59       for (ispr = 0; ispr <= mkbmx2-1; ispr++)
  60           EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
  61
  62       EMFSTK.louemf[npnw] = TRACKR.llouse;
  63   }
  64
  65
  66   Int_t verbosityLevel = fluka->GetVerbosityLevel();
  67   Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
  68
  69   fluka->SetTrackIsNew(kTRUE);
  70
  71 // Get the stack produced from the generator
  72   TVirtualMCStack* cppstack = fluka->GetStack();
  73
  74 // EVTFLG.ntrcks = track number
  75 // Increment the track number and put it into the last flag
  76
  77   Int_t kp;
  78
  79   for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
  80
  81 // Check transport cut first
  82     Int_t    ireg   = EMFSTK.iremf[kp];
  83     Double_t cut    = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EMFRGN.elethr[ireg-1] :  EMFRGN.phothr[ireg-1];
  84     Double_t e      = EMFSTK.etemf[kp];
  85
  86     if ((e < cut)
  87         && (
  88             (EMFSTK.ichemf[kp] ==  0) ||
  89             (EMFSTK.ichemf[kp] == -1) ||
  90             (EMFSTK.ichemf[kp] ==  1 &&  EMFRGN.phothr[ireg-1] > emassmev)
  91             )
  92         )
  93     {
  94         EMFSTK.iespak[kp][mkbmx2-1] = TRACKR.ispusr[mkbmx2-1];
  95         EMFSTK.iespak[kp][mkbmx2-2] =  0;
  96         EMFSTK.iespak[kp][mkbmx2-3] = TRACKR.jtrack;
  97         EMFSTK.irlatt[kp] = TRACKR.lt1trk;
  98         continue;
  99     }
 100
 101 // Save the parent track number and reset it at each loop
 102     Int_t done = 0;
 103     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
 104     Int_t flukaid = 0;
 105
 106 // Identify particle type
 107     if      (EMFSTK.ichemf[kp] == -1)  flukaid = kFLUKAelectron;
 108     else if (EMFSTK.ichemf[kp] ==  0)  flukaid = kFLUKAphoton;
 109     else if (EMFSTK.ichemf[kp] ==  1)  flukaid = kFLUKApositron;
 110
 111
 112     e *= emvgev;
 113     Int_t    pdg    = fluka->PDGFromId(flukaid);
 114     Double_t p      = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
 115     Double_t px     = p * EMFSTK.uemf[kp];
 116     Double_t py     = p * EMFSTK.vemf[kp];
 117     Double_t pz     = p * EMFSTK.wemf[kp];
 118     Double_t tof    = EMFSTK.agemf[kp];
 119     Double_t polx   = EMFSTK.upol[kp];
 120     Double_t poly   = EMFSTK.vpol[kp];
 121     Double_t polz   = EMFSTK.wpol[kp];
 122     Double_t vx     = EMFSTK.xemf[kp];
 123     Double_t vy     = EMFSTK.yemf[kp];
 124     Double_t vz     = EMFSTK.zemf[kp];
 125     Double_t weight = EMFSTK.wtemf[kp];
 126
 127     Int_t ntr;
 128     TMCProcess mech;
 129     Int_t is = 0;
 130
 131 //* case of no parent left (pair, photoelectric, annihilation):
 132 //* all secondaries are true
 133     if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
 134         (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
 135
 136         if (EVTFLG.lpairp == 1) mech = kPPair;
 137         else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
 138         else mech = kPAnnihilation;
 139         cppstack->PushTrack(done, parent, pdg,
 140                            px, py, pz, e, vx, vy, vz, tof,
 141                            polx, poly, polz, mech, ntr, weight, is);
 142         if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << " energy " << e-PAPROP.am[flukaid+6] << endl;
 143
 144         EMFSTK.iespak[kp][mkbmx2-1] = ntr;
 145         EMFSTK.iespak[kp][mkbmx2-2] = 0;
 146     } // end of lpairp, lphoel, lannfl, lannrs
 147
 148 //* Compton: secondary is true only if charged (e+, e-)
 149     else if ((EVTFLG.lcmptn == 1)) {
 150
 151         if (EMFSTK.ichemf[kp] != 0) {
 152             mech = kPCompton;
 153             cppstack->PushTrack(done, parent, pdg,
 154                                px, py, pz, e, vx, vy, vz, tof,
 155                                polx, poly, polz, mech, ntr, weight, is);
 156             if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
 157             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
 158             EMFSTK.iespak[kp][mkbmx2-2] = 0;
 159         } else {
 160             fluka->SetPint(px, py, pz, e);
 161         }
 162     } // end of lcmptn
 163
 164 //* Bremsstrahlung: true secondary only if charge = 0 (photon)
 165     else if ((EVTFLG.lbrmsp == 1)) {
 166         if (EMFSTK.ichemf[kp] == 0) {
 167             mech = kPBrem;
 168             cppstack->PushTrack(done, parent, pdg,
 169                                px, py, pz, e, vx, vy, vz, tof,
 170                                polx, poly, polz, mech, ntr, weight, is);
 171             if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
 172             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
 173             EMFSTK.iespak[kp][mkbmx2-2] = 0;
 174         } else {
 175             fluka->SetPint(px, py, pz, e);
 176         }
 177
 178     } // end of lbrmsp
 179
 180 //* Delta ray: If Bhabha, true secondary only if negative (electron)
 181     else if ((EVTFLG.ldltry == 1)) {
 182         if (lbhabh == 1) {
 183             if (EMFSTK.ichemf[kp] == -1) {
 184                 mech = kPDeltaRay;
 185                 cppstack->PushTrack(done, parent, pdg,
 186                                     px, py, pz, e, vx, vy, vz, tof,
 187                                     polx, poly, polz, mech, ntr, weight, is);
 188                 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
 189                 EMFSTK.iespak[kp][mkbmx2-2] = 0;
 190                 if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
 191             } else {
 192                 fluka->SetPint(px, py, pz, e);
 193             }
 194         } // lbhabh == 1
 195
 196 //* Delta ray: Otherwise Moller: true secondary is the electron with
 197 //*            lower energy, which has been put higher in the stack
 198         else if (kp == EMFSTK.npemf-1) {
 199             mech = kPDeltaRay;
 200             cppstack->PushTrack(done, parent, pdg,
 201                                px, py, pz, e, vx, vy, vz, tof,
 202                                polx, poly, polz, mech, ntr, weight, is);
 203             if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
 204             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
 205             EMFSTK.iespak[kp][mkbmx2-2] = 0;
 206         } else {
 207             fluka->SetPint(px, py, pz, e);
 208         }
 209     } // end of ldltry
 210
 211   } // end of loop
 212 } // end of stupre
 213 } // end of extern "C"
 214