TFluka/stuprf.cxx

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