TFluka/stuprf.cxx

   1 #include <Riostream.h>
   2 #ifndef WIN32
   3 # define stuprf stuprf_
   4 #else
   5 # define stuprf STUPRF
   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 "Fevtflg.h"  //(EVTFLG) fluka common
  13 #include "Fpaprop.h"  //(PAPROP) fluka common
  14 #include "Fflkstk.h"  //(FLKSTK)  fluka common
  15 #include "Ftrackr.h"  //(TRACKR) fluka common
  16 #include "Fgenstk.h"  //(GENSTK)  fluka common
  17
  18
  19 //Virtual MC
  20 #include "TFluka.h"
  21
  22 #include "TVirtualMCStack.h"
  23 #include "TVirtualMCApplication.h"
  24 #include "TParticle.h"
  25 #include "TVector3.h"
  26
  27 extern "C" {
  28     void stuprf(Int_t& /*ij*/, Int_t& /*mreg*/,
  29                 Double_t& xx, Double_t& yy, Double_t& zz,
  30                 Int_t& numsec, Int_t& npprmr)
  31 {
  32 //*----------------------------------------------------------------------*
  33 //*                                                                      *
  34 //*  SeT User PRoperties for Fluka particles                             *
  35 //*                                                                      *
  36 //*----------------------------------------------------------------------*
  37
  38 // FLKSTK.npflka  = stack pointer
  39 // FLKSTK.louse   = user flag
  40 // TRACKR.llouse = user defined flag for the current particle
  41
  42    FLKSTK.louse[FLKSTK.npflka] = 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 // FLKSTK.sparek  = spare real variables available for k.w.burn
  47 // FLKSTK.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   if (numsec <= npprmr) {
  52       for (ispr = 0; ispr <= mkbmx1 - 1; ispr++) {
  53           FLKSTK.sparek[FLKSTK.npflka][ispr] = TRACKR.spausr[ispr];
  54           TRACKR.spausr[ispr] = 0.;
  55       }
  56       for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
  57           FLKSTK.ispark[FLKSTK.npflka][ispr] = TRACKR.ispusr[ispr];
  58           TRACKR.ispusr[ispr] = 0;
  59       }
  60   }
  61
  62   // save parent info
  63   FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 3] = TRACKR.jtrack;              // fluka particle id
  64   FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 4] = TRACKR.ispusr[mkbmx2 - 1];  // current track number
  65   FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 5] = npprmr;                     // flag npprmr>0
  66
  67 // Get the pointer to the VMC
  68   TFluka* fluka =  (TFluka*) gMC;
  69   Int_t verbosityLevel = fluka->GetVerbosityLevel();
  70   Bool_t debug = (verbosityLevel>=3)? kTRUE:kFALSE;
  71
  72   fluka->SetTrackIsNew(kTRUE);
  73 //  TVirtualMC* fluka = TFluka::GetMC();
  74 // Get the stack produced from the generator
  75   TVirtualMCStack* cppstack = fluka->GetStack();
  76
  77 // EVTFLG.ntrcks = track number
  78 // Increment the track number and put it into the last flag
  79 // was numsec -1
  80 // clarify with Alberto
  81
  82 //  npprmr > 0, the secondary being loaded is actually still the interacting
  83 //  particle (it can happen in some biasing situations)
  84
  85   if (numsec > npprmr) {
  86 // Now call the PushTrack(...)
  87     Int_t done = 0;
  88
  89     Int_t parent  =  TRACKR.ispusr[mkbmx2-1];
  90     Int_t kpart   = GENSTK.kpart[numsec-1];
  91     if (kpart < -6) return; // -7 to -12 = "heavy" fragment
  92     Int_t  pdg  = fluka->PDGFromId(kpart);
  93
  94     Double_t px = GENSTK.plr[numsec-1] * GENSTK.cxr[numsec-1];
  95     Double_t py = GENSTK.plr[numsec-1] * GENSTK.cyr[numsec-1];
  96     Double_t pz = GENSTK.plr[numsec-1] * GENSTK.czr[numsec-1];
  97     Double_t e  = GENSTK.tki[numsec-1] + PAPROP.am[GENSTK.kpart[numsec-1]+6];
  98
  99     Double_t vx = xx;
 100     Double_t vy = yy;
 101     Double_t vz = zz;
 102
 103     Double_t tof  = TRACKR.atrack;
 104     Double_t polx = GENSTK.cxrpol[numsec-1];
 105     Double_t poly = GENSTK.cyrpol[numsec-1];
 106     Double_t polz = GENSTK.czrpol[numsec-1];
 107
 108
 109
 110     TMCProcess mech = kPHadronic;
 111     if (EVTFLG.ldecay == 1) {
 112         mech = kPDecay;
 113         if (debug) cout << endl << "Decay" << endl;
 114         FLKSTK.nlattc[FLKSTK.npflka] = TRACKR.lt1trk;
 115     } else if (EVTFLG.ldltry == 1) {
 116         mech = kPDeltaRay;
 117         if( fluka->GetIcode() == kKASHEA ) {
 118            //  For all interactions secondaries are put on GENSTK common (kp=1,np)
 119            //  but for KASHEA delta ray generation where only the secondary elec-
 120            //  tron is present and stacked on FLKSTK common for kp=lstack
 121            pdg  = fluka->PDGFromId( FLKSTK.iloflk[FLKSTK.npflka] );
 122            px   = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.txflk[FLKSTK.npflka];
 123            py   = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.tyflk[FLKSTK.npflka];
 124            pz   = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.tzflk[FLKSTK.npflka];
 125            e    = FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.am[FLKSTK.iloflk[FLKSTK.npflka]+6];
 126            polx = FLKSTK.txpol[FLKSTK.npflka];
 127            poly = FLKSTK.typol[FLKSTK.npflka];
 128            polz = FLKSTK.tzpol[FLKSTK.npflka];
 129            if (debug) cout << endl << "Delta Ray from KASHEA...." << " pdg from FLKSTK=" << pdg << endl;
 130         } else {
 131             if (debug) cout << endl << "Delta Ray" << endl;
 132         }
 133     } else if (EVTFLG.lpairp == 1) {
 134         mech = kPPair;
 135         if (debug) cout << endl << "Pair Production" << endl;
 136     } else if (EVTFLG.lbrmsp == 1) {
 137         mech = kPBrem;
 138         if (debug) cout << endl << "Bremsstrahlung" << endl;
 139     }
 140
 141     Double_t weight = GENSTK.wei[numsec-1];
 142     Int_t is = 0;
 143     Int_t ntr;
 144     //
 145     // Save particle in VMC stack
 146     cppstack->PushTrack(done, parent, pdg,
 147                         px, py, pz, e,
 148                         vx, vy, vz, tof,
 149                         polx, poly, polz,
 150                         mech, ntr, weight, is);
 151     if (debug)
 152         cout << endl << " !!! stuprf: ntr=" << ntr << " pdg " << pdg << " parent=" << parent
 153              << " parent_pdg="<< fluka->PDGFromId(TRACKR.jtrack) << " numsec "
 154              << numsec << " npprmr " << npprmr << " icode=" << fluka->GetIcode()
 155              << "kin. energy [keV] = " <<  GENSTK.tki[numsec-1] * 1.e6 << endl << endl;
 156
 157
 158 //
 159 //  Save current track number
 160     FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = ntr;
 161     FLKSTK.ispark[FLKSTK.npflka][mkbmx2-2] = 0;
 162   } // end of if (numsec > npprmr)
 163 //  else {
 164 //     if(debug) {
 165 //        cout << endl << " !!! stuprf: skipping pushtrack   track=" << TRACKR.ispusr[mkbmx2-1]
 166 //              << " pdg " << fluka->PDGFromId(TRACKR.jtrack) << " numsec=" << numsec<< " npprmr=" << npprmr
 167 //              << " GENSTK pdg=" << fluka->PDGFromId(GENSTK.kpart[numsec-1]) << endl;
 168 //     }
 169 //  }
 170 } // end of stuprf
 171 } // end of extern "C"
 172