1+ // Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+ // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+ // All rights not expressly granted are reserved.
4+ //
5+ // This software is distributed under the terms of the GNU General Public
6+ // License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+ //
8+ // In applying this license CERN does not waive the privileges and immunities
9+ // granted to it by virtue of its status as an Intergovernmental Organization
10+ // or submit itself to any jurisdiction.
11+ //
12+ // Author: J. E. Munoz Mendez jesus.munoz@cern.ch
13+
14+ std ::function < void (const double * double * )> field ()
15+ {
16+ return [](const double * x , double * b ) {
17+ const double Rc = 185. ; //[cm]
18+ const double R1 = 220. ; //[cm]
19+ const double R2 = 290. ; //[cm]
20+ const double B1 = 2. ; //[T]
21+ const double B2 = - Rc * Rc / ((R2 * R2 - R1 * R1 ) * B1 ); //[T]
22+ const double beamStart = 370. ; //[cm]
23+ const double tokGauss = 1. / 0.1 ; // conversion from Tesla to kGauss
24+
25+ const bool isMagAbs = true;
26+
27+ const double r = sqrt (x [0 ] * x [0 ] + x [1 ] * x [1 ]);
28+ if ((abs (x [2 ]) <= beamStart ) && (r < Rc )) { // We are inside of the central region
29+ b [0 ] = 0. ;
30+ b [1 ] = 0. ;
31+ b [2 ] = B1 * tokGauss ;
32+ } else if ((abs (x [2 ]) <= beamStart ) && (r >= Rc && r < R1 )) { // We are in the transition region
33+ b [0 ] = 0. ;
34+ b [1 ] = 0. ;
35+ b [2 ] = 0. ;
36+ } else if ((abs (x [2 ]) <= beamStart ) && (r >= R1 && r < R2 )) { // We are within the magnet
37+ b [0 ] = 0. ;
38+ b [1 ] = 0. ;
39+ if (isMagAbs ) {
40+ b [2 ] = B2 * tokGauss ;
41+ } else {
42+ b [2 ] = 0. ;
43+ }
44+ } else { // We are outside of the magnet
45+ b [0 ] = 0. ;
46+ b [1 ] = 0. ;
47+ b [2 ] = 0. ;
48+ }
49+ };
50+ }
51+
52+ void ALICE3V3Magnet ()
53+ {
54+ auto fieldFunc = field ();
55+ // RZ plane visualization
56+ TCanvas * cRZ = new TCanvas ("cRZ" , "Field in RZ plane" , 800 , 600 );
57+ TH2F * hRZ = new TH2F ("hRZ" , "Magnetic Field B_z in RZ plane;Z [m];R [m]" , 100 , -10 , 10 , 100 , -5 , 5 );
58+ hRZ -> SetBit (TH1 ::kNoStats ); // disable stats box
59+ for (int i = 1 ; i <= hRZ -> GetNbinsX (); i ++ ) {
60+ const double Z = hRZ -> GetXaxis ()-> GetBinCenter (i );
61+ for (int j = 1 ; j <= hRZ -> GetNbinsY (); j ++ ) {
62+ const double R = hRZ -> GetYaxis ()-> GetBinCenter (j );
63+ const double pos [3 ] = {R * 100 , 0 , Z * 100 }; // convert to cm
64+ double b [3 ] = {0 , 0 , 0 };
65+ fieldFunc (pos , b );
66+ hRZ -> SetBinContent (i , j , b [2 ]);
67+ }
68+ }
69+
70+ hRZ -> Draw ("COLZ" );
71+ cRZ -> Update ();
72+
73+ // XY plane visualization
74+ TCanvas * cXY = new TCanvas ("cXY" , "Field in XY plane" , 800 , 600 );
75+ TH2F * hXY = new TH2F ("hXY" , "Magnetic Field B_z in XY plane;X [m];Y [m]" , 100 , -5 , 5 , 100 , -5 , 5 );
76+ hXY -> SetBit (TH1 ::kNoStats ); // disable stats box
77+
78+ for (int i = 1 ; i <= hXY -> GetNbinsX (); i ++ ) {
79+ const double X = hXY -> GetXaxis ()-> GetBinCenter (i );
80+ for (int j = 1 ; j <= hXY -> GetNbinsY (); j ++ ) {
81+ const double Y = hXY -> GetYaxis ()-> GetBinCenter (j );
82+ const double pos [3 ] = {X * 100 , Y * 100 , 0 }; // convert to cm
83+ double b [3 ] = {0 , 0 , 0 };
84+ fieldFunc (pos , b );
85+ hXY -> SetBinContent (i , j , b [2 ]);
86+ }
87+ }
88+
89+ hXY -> Draw ("COLZ" );
90+ cXY -> Update ();
91+ }
0 commit comments