No debug message.
[u/mrichter/AliRoot.git] / 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