TFluka/magfld.cxx

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