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 // Fluka commons
  11 #include "Fdblprc.h"  //(DBLPRC) fluka common
  12 #include "Femfstk.h"  //(EMFSTK) fluka common
  13 #include "Fevtflg.h"  //(EVTFLG) fluka common
  14 #include "Fpaprop.h"  //(PAPROP) fluka common
  15 #include "Ftrackr.h"  //(TRACKR) fluka common
  16
  17 //Virtual MC
  18
  19 #include "TFluka.h"
  20
  21 #include "TVirtualMCStack.h"
  22 #include "TVirtualMCApplication.h"
  23 #include "TParticle.h"
  24 #include "TVector3.h"
  25
  26 extern "C" {
  27 void stupre()
  28 {
  29 //*----------------------------------------------------------------------*
  30 //*                                                                      *
  31 //*  SeT User PRoperties for Emf particles                               *
  32 //*                                                                      *
  33 //*----------------------------------------------------------------------*
  34   Int_t lbhabh = 0;
  35   if (EVTFLG.ldltry == 1) {
  36     if (EMFSTK.ichemf[EMFSTK.npemf-1] * EMFSTK.ichemf[EMFSTK.npemf-2] < 0) lbhabh = 1;
  37   }
  38
  39 // mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
  40 // mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
  41 // EMFSTK.espark  = spare real variables available for
  42 // EMFSTK.iespak  = spare integer variables available for
  43 // TRACKR.spausr = user defined spare variables for the current particle
  44 // TRACKR.ispusr = user defined spare flags for the current particle
  45 // EMFSTK.louemf = user flag
  46 // TRACKR.llouse = user defined flag for the current particle
  47
  48   Int_t npnw, ispr;
  49   for (npnw=EMFSTK.npstrt-1; npnw<=EMFSTK.npemf-1; npnw++) {
  50
  51     for (ispr=0; ispr<=mkbmx1-1; ispr++)
  52       EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
  53
  54     for (ispr=0; ispr<=mkbmx2-1; ispr++)
  55       EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
  56
  57     EMFSTK.louemf[npnw] = TRACKR.llouse;
  58   }
  59
  60 // Get the pointer to the VMC
  61   TFluka* fluka =  (TFluka*) gMC;
  62   Int_t verbosityLevel = fluka->GetVerbosityLevel();
  63   Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
  64   fluka->SetTrackIsNew(kTRUE);
  65 //  TVirtualMC* fluka = TFluka::GetMC();
  66 // Get the stack produced from the generator
  67   TVirtualMCStack* cppstack = fluka->GetStack();
  68
  69 // EVTFLG.ntrcks = track number
  70 // Increment the track number and put it into the last flag
  71
  72   Int_t kp;
  73   for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
  74
  75 //* save the parent track number and reset it at each loop
  76     Int_t done = 0;
  77
  78     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
  79
  80     Int_t flukaid = 0;
  81
  82     if (EMFSTK.ichemf[kp] == -1) flukaid = 3;
  83     else if (EMFSTK.ichemf[kp] == 0)  flukaid = 7;
  84     else if (EMFSTK.ichemf[kp] == 1)  flukaid = 4;
  85
  86     Int_t pdg       = fluka->PDGFromId(flukaid);
  87     Double_t e      = EMFSTK.etemf[kp] * emvgev;
  88     if (flukaid + 6 < 0) printf("stupre: Calling PDGFromId for %10d %10d  \n", kp, flukaid);
  89     Double_t p      = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
  90     Double_t px     = p * EMFSTK.u[kp];
  91     Double_t pz     = p * EMFSTK.v[kp];
  92     Double_t py     = p * EMFSTK.w[kp];
  93     Double_t tof    = EMFSTK.agemf[kp];
  94     Double_t polx   = EMFSTK.upol[kp];
  95     Double_t poly   = EMFSTK.vpol[kp];
  96     Double_t polz   = EMFSTK.wpol[kp];
  97     Double_t vx     = EMFSTK.x[kp];
  98     Double_t vy     = EMFSTK.y[kp];
  99     Double_t vz     = EMFSTK.z[kp];
 100     Double_t weight = EMFSTK.wtemf[kp];
 101
 102     Int_t ntr;
 103     TMCProcess mech;
 104     Int_t is = 0;
 105
 106 //* case of no parent left (pair, photoelectric, annihilation):
 107 //* all secondaries are true
 108     if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
 109         (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
 110
 111         if (EVTFLG.lpairp == 1) mech = kPPair;
 112         else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
 113         else mech = kPAnnihilation;
 114         cppstack->PushTrack(done, parent, pdg,
 115                            px, py, pz, e, vx, vy, vz, tof,
 116                            polx, poly, polz, mech, ntr, weight, is);
 117         if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
 118
 119         EMFSTK.iespak[kp][mkbmx2-1] = ntr;
 120         EMFSTK.iespak[kp][mkbmx2-2] = 0;
 121     } // end of lpairp, lphoel, lannfl, lannrs
 122
 123 //* Compton: secondary is true only if charged (e+, e-)
 124     else if ((EVTFLG.lcmptn == 1)) {
 125
 126         if (EMFSTK.ichemf[kp] != 0) {
 127             mech = kPCompton;
 128             cppstack->PushTrack(done, parent, pdg,
 129                                px, py, pz, e, vx, vy, vz, tof,
 130                                polx, poly, polz, mech, ntr, weight, is);
 131             if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
 132             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
 133             EMFSTK.iespak[kp][mkbmx2-2] = 0;
 134         }
 135     } // end of lcmptn
 136
 137 //* Bremsstrahlung: true secondary only if charge = 0 (photon)
 138     else if ((EVTFLG.lbrmsp == 1)) {
 139         if (EMFSTK.ichemf[kp] == 0) {
 140             mech = kPBrem;
 141             cppstack->PushTrack(done, parent, pdg,
 142                                px, py, pz, e, vx, vy, vz, tof,
 143                                polx, poly, polz, mech, ntr, weight, is);
 144             if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
 145             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
 146             EMFSTK.iespak[kp][mkbmx2-2] = 0;
 147         }
 148     } // end of lbrmsp
 149
 150 //* Delta ray: If Bhabha, true secondary only if negative (electron)
 151     else if ((EVTFLG.ldltry == 1)) {
 152         if (lbhabh == 1) {
 153             if (EMFSTK.ichemf[kp] == -1) {
 154                 mech = kPDeltaRay;
 155                 cppstack->PushTrack(done, parent, pdg,
 156                                    px, py, pz, e, vx, vy, vz, tof,
 157                                    polx, poly, polz, mech, ntr, weight, is);
 158                 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
 159                 EMFSTK.iespak[kp][mkbmx2-2] = 0;
 160            if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
 161             } // end of Bhabha
 162         } // lbhabh == 1
 163
 164 //* Delta ray: Otherwise Moller: true secondary is the electron with
 165 //*            lower energy, which has been put higher in the stack
 166         else if (kp == EMFSTK.npemf-1) {
 167             mech = kPDeltaRay;
 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 (Moller) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
 172             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
 173             EMFSTK.iespak[kp][mkbmx2-2] = 0;
 174         } // end of Delta ray
 175     } // end of ldltry
 176
 177   } // end of loop
 178
 179 // !!! TO BE CONFIRMED !!!
 180 } // end of stupre
 181 } // end of extern "C"
 182