TFluka/magfld.cxx

   1 #include "TVirtualMCApplication.h"
   2 #include "Fdblprc.h"  //(DBLPRC) fluka common
   3 //
   4 // #include "TCallf77.h"
   5
   6 #ifndef WIN32
   7 #define magfld magfld_
   8 #define type_of_call
   9 #else
  10 #define magfld MAGFLD
  11 #define type_of_call  _stdcall
  12 #endif
  13
  14 extern "C" void type_of_call magfld(double& x,   double& y,   double& z,
  15                                     double& btx, double& bty, double& btz, double& b,
  16                                     int&    nreg,int& idisc)
  17 {
  18
  19 /*
  20 *----------------------------------------------------------------------*
  21 *                                                                      *
  22 *                                                                      *
  23 *     Input variables:                                                 *
  24 *            x,y,z = current position                                  *
  25 *            nreg  = current region                                    *
  26 *     Output variables:                                                *
  27 *            btx,bty,btz = cosines of the magn. field vector           *
  28 *            B = magnetic field intensity (Tesla)                      *
  29 *            idisc = set to 1 if the particle has to be discarded      *
  30 *                                                                      *
  31 *----------------------------------------------------------------------*
  32 */
  33
  34     printf("\n ========> magfld\n");
  35
  36
  37     idisc = 0;
  38
  39     Double_t bc[3];
  40     Double_t xc[3];
  41
  42     xc[1] = x;
  43     xc[0] = y;
  44     xc[2] = z;
  45
  46
  47
  48     (TVirtualMCApplication::Instance())->Field(xc, bc);
  49
  50     b = sqrt(bc[0] * bc[0] + bc[1] * bc[1] + bc[2] * bc[2]);
  51     if (b) {
  52         btx = bc[1]/b;
  53         bty = bc[0]/b;
  54         Double_t btt = btx * btx + bty * bty;
  55         if (btt >= (Double_t) 1.) {
  56             btx /= TMath::Sqrt(btt);
  57             bty /= TMath::Sqrt(btt);
  58             b   /= TMath::Sqrt(btt);
  59             btz =  (Double_t) 0.;
  60         } else {
  61             btz = TMath::Sqrt((Double_t) 1. -  btt);
  62         }
  63     } else {
  64         btx = 0.;
  65         bty = 0.;
  66         btz = 1.;
  67     }
  68
  69     // from kG to T
  70     b /= (Double_t) 10.;
  71     printf("\n #magfld: x, y, z, btx, bty, btz, b %f %f %f %f %f %f %f \n",
  72            x, y, z, btx, bty, btz, b);
  73 }