1111//
1212// Author: J. E. Munoz Mendez jesus.munoz@cern.ch
1313
14+ #include <functional>
15+ #include <cmath>
16+ #include <TCanvas.h>
17+ #include <TH2F.h>
18+ #include <TStyle.h>
19+
1420std ::function < void (const double * , double * )> field ()
1521{
1622 return [](const double * x , double * b ) {
17- double Rc ;
18- double R1 ;
19- double R2 ;
20- double B1 ;
21- double B2 ;
22- const double beamStart = 500. ; //[cm]
23- const double tokGauss = 1. / 0.1 ; // conversion from Tesla to kGauss
24-
25- bool isMagAbs = true;
26-
27- // ***********************
28- // LAYOUT 1
29- // ***********************
30-
3123 // RADIUS
32- Rc = 185 . ; //[cm]
33- R1 = 220. ; //[cm]
34- R2 = 290. ; //[cm]
24+ static constexpr double Rc = 170 . ; // [cm] — R_out_coil per Ian DetectorConstruction.cc; confirmed by A. Ortiz definition
25+ static constexpr double R1 = 220. ; // [cm]
26+ static constexpr double R2 = 290. ; // [cm]
3527
36- // To set the B2
37- B1 = 2. ; //[T]
38- B2 = - Rc * Rc / ((R2 * R2 - R1 * R1 ) * B1 ); //[T]
28+ // FIELD
29+ static constexpr double B1 = 2. ; // [T]
30+ static constexpr double B2 = - B1 * Rc * Rc / (R2 * R2 - R1 * R1 ); // [T] — B1 in numerator, confirmed by A. Ortiz Aug 2026
31+ static constexpr double beamStart = 500. ; // [cm]
32+ static constexpr double tokGauss = 1. / 0.1 ; // conversion from Tesla to kGauss
3933
40- if ((abs (x [2 ]) <= beamStart ) && (sqrt (x [0 ] * x [0 ] + x [1 ] * x [1 ]) < Rc )) {
34+ static constexpr bool isMagAbs = true;
35+
36+ const double r = sqrt (x [0 ] * x [0 ] + x [1 ] * x [1 ]);
37+ if ((abs (x [2 ]) <= beamStart ) && (r < Rc )) { // We are inside of the central region
4138 b [0 ] = 0. ;
4239 b [1 ] = 0. ;
4340 b [2 ] = B1 * tokGauss ;
44- } else if ((abs (x [2 ]) <= beamStart ) &&
45- (sqrt (x [0 ] * x [0 ] + x [1 ] * x [1 ]) >= Rc &&
46- sqrt (x [0 ] * x [0 ] + x [1 ] * x [1 ]) < R1 )) {
41+ } else if ((abs (x [2 ]) <= beamStart ) && (r >= Rc && r < R1 )) { // We are in the transition region
4742 b [0 ] = 0. ;
4843 b [1 ] = 0. ;
4944 b [2 ] = 0. ;
50- } else if ((abs (x [2 ]) <= beamStart ) &&
51- (sqrt (x [0 ] * x [0 ] + x [1 ] * x [1 ]) >= R1 &&
52- sqrt (x [0 ] * x [0 ] + x [1 ] * x [1 ]) < R2 )) {
45+ } else if ((abs (x [2 ]) <= beamStart ) && (r >= R1 && r < R2 )) { // We are within the magnet
5346 b [0 ] = 0. ;
5447 b [1 ] = 0. ;
5548 if (isMagAbs ) {
5649 b [2 ] = B2 * tokGauss ;
5750 } else {
5851 b [2 ] = 0. ;
5952 }
60- } else {
53+ } else { // We are outside of the magnet
6154 b [0 ] = 0. ;
6255 b [1 ] = 0. ;
6356 b [2 ] = 0. ;
@@ -67,10 +60,14 @@ std::function<void(const double*, double*)> field()
6760
6861void ALICE3Field ()
6962{
63+ gStyle -> SetPalette (kRainBow );
64+ gStyle -> SetNumberContours (255 );
65+
7066 auto fieldFunc = field ();
7167 // RZ plane visualization
72- TCanvas * cRZ = new TCanvas ("cRZ" , "Field in RZ plane" , 800 , 600 );
73- TH2F * hRZ = new TH2F ("hRZ" , "Magnetic Field B_z in RZ plane;Z [m];R [m]" , 100 , -10 , 10 , 100 , -5 , 5 );
68+ TCanvas * cRZ = new TCanvas ("cRZ" , "Field in RZ plane" , 800 , 800 );
69+ gPad -> SetRightMargin (0.15 );
70+ TH2F * hRZ = new TH2F ("hRZ" , "Magnetic Field B_z in RZ plane;Z [m];R [m];B_{z} [kGauss]" , 100 , -10 , 10 , 100 , -5 , 5 );
7471 hRZ -> SetBit (TH1 ::kNoStats ); // disable stats box
7572 for (int i = 1 ; i <= hRZ -> GetNbinsX (); i ++ ) {
7673 const double Z = hRZ -> GetXaxis ()-> GetBinCenter (i );
@@ -83,12 +80,14 @@ void ALICE3Field()
8380 }
8481 }
8582
83+ hRZ -> GetZaxis ()-> SetRangeUser (-30 , 30 );
8684 hRZ -> Draw ("COLZ" );
8785 cRZ -> Update ();
8886
8987 // XY plane visualization
90- TCanvas * cXY = new TCanvas ("cXY" , "Field in XY plane" , 800 , 600 );
91- TH2F * hXY = new TH2F ("hXY" , "Magnetic Field B_z in XY plane;X [m];Y [m]" , 100 , -5 , 5 , 100 , -5 , 5 );
88+ TCanvas * cXY = new TCanvas ("cXY" , "Field in XY plane" , 800 , 800 );
89+ gPad -> SetRightMargin (0.15 );
90+ TH2F * hXY = new TH2F ("hXY" , "Magnetic Field B_z in XY plane;X [m];Y [m];B_{z} [kGauss]" , 100 , -5 , 5 , 100 , -5 , 5 );
9291 hXY -> SetBit (TH1 ::kNoStats ); // disable stats box
9392
9493 for (int i = 1 ; i <= hXY -> GetNbinsX (); i ++ ) {
@@ -102,6 +101,7 @@ void ALICE3Field()
102101 }
103102 }
104103
104+ hXY -> GetZaxis ()-> SetRangeUser (-30 , 30 );
105105 hXY -> Draw ("COLZ" );
106106 cXY -> Update ();
107- }
107+ }
0 commit comments