Bug fix. Removed delete statement
[u/mrichter/AliRoot.git] / STEER / AliMagF.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
4c039060 16
db83d72f 17#include <TClass.h>
18#include <TFile.h>
19#include <TSystem.h>
fe4da5cc 20
21#include "AliMagF.h"
db83d72f 22#include "AliMagWrapCheb.h"
23#include "AliLog.h"
972ca52f 24
fe4da5cc 25ClassImp(AliMagF)
26
9251fceb 27const Double_t AliMagF::fgkSol2DipZ = -700.;
db83d72f 28
e2afb3b6 29//_______________________________________________________________________
30AliMagF::AliMagF():
db83d72f 31 TVirtualMagField(),
32 fMeasuredMap(0),
33 fMapType(k5kG),
34 fSolenoid(0),
35 fBeamType(kNoBeamField),
36 fBeamEnergy(0),
db83d72f 37 //
e2afb3b6 38 fInteg(0),
db83d72f 39 fPrecInteg(0),
40 fFactorSol(1.),
41 fFactorDip(1.),
42 fMax(15),
43 fDipoleOFF(kFALSE),
e2afb3b6 44 //
db83d72f 45 fQuadGradient(0),
46 fDipoleField(0),
47 fCCorrField(0),
48 fACorr1Field(0),
49 fACorr2Field(0),
50 fParNames("","")
51{
e2afb3b6 52 // Default constructor
53 //
54}
55
56//_______________________________________________________________________
db83d72f 57AliMagF::AliMagF(const char *name, const char* title, Int_t integ,
58 Double_t factorSol, Double_t factorDip,
59 Double_t fmax, BMap_t maptype, const char* path,
9251fceb 60 BeamType_t bt, Double_t be):
db83d72f 61 TVirtualMagField(name),
62 fMeasuredMap(0),
63 fMapType(maptype),
64 fSolenoid(0),
65 fBeamType(bt),
66 fBeamEnergy(be),
db83d72f 67 //
68 fInteg(integ),
604e0531 69 fPrecInteg(1),
db83d72f 70 fFactorSol(1.),
71 fFactorDip(1.),
972ca52f 72 fMax(fmax),
db83d72f 73 fDipoleOFF(factorDip==0.),
74 //
75 fQuadGradient(0),
76 fDipoleField(0),
77 fCCorrField(0),
78 fACorr1Field(0),
79 fACorr2Field(0),
80 fParNames("","")
fe4da5cc 81{
9251fceb 82 // Initialize the field with Geant integration option "integ" and max field "fmax,
83 // Impose scaling of parameterized L3 field by factorSol and of dipole by factorDip.
84 // The "be" is the energy of the beam in GeV/nucleon
aee8290b 85 //
db83d72f 86 SetTitle(title);
87 if(integ<0 || integ > 2) {
88 AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
89 fInteg = 2;
90 }
91 if (fInteg == 0) fPrecInteg = 0;
aee8290b 92 //
db83d72f 93 const char* parname = 0;
94 //
f04e7f5f 95 if (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
96 else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
97 else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
98 else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
db83d72f 99 //
100 SetDataFileName(path);
101 SetParamName(parname);
102 //
db83d72f 103 LoadParameterization();
104 InitMachineField(fBeamType,fBeamEnergy);
f04e7f5f 105 double xyz[3]={0.,0.,0.};
106 fSolenoid = GetBz(xyz);
107 SetFactorSol(factorSol);
108 SetFactorDip(factorDip);
fe4da5cc 109}
110
e2afb3b6 111//_______________________________________________________________________
eeda4611 112AliMagF::AliMagF(const AliMagF &src):
db83d72f 113 TVirtualMagField(src),
114 fMeasuredMap(0),
115 fMapType(src.fMapType),
116 fSolenoid(src.fSolenoid),
117 fBeamType(src.fBeamType),
118 fBeamEnergy(src.fBeamEnergy),
eeda4611 119 fInteg(src.fInteg),
120 fPrecInteg(src.fPrecInteg),
db83d72f 121 fFactorSol(src.fFactorSol),
122 fFactorDip(src.fFactorDip),
eeda4611 123 fMax(src.fMax),
db83d72f 124 fDipoleOFF(src.fDipoleOFF),
125 fQuadGradient(src.fQuadGradient),
126 fDipoleField(src.fDipoleField),
127 fCCorrField(src.fCCorrField),
128 fACorr1Field(src.fACorr1Field),
129 fACorr2Field(src.fACorr2Field),
130 fParNames(src.fParNames)
eeda4611 131{
db83d72f 132 if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
eeda4611 133}
134
135//_______________________________________________________________________
db83d72f 136AliMagF::~AliMagF()
ff66b122 137{
db83d72f 138 delete fMeasuredMap;
139}
140
141//_______________________________________________________________________
142Bool_t AliMagF::LoadParameterization()
143{
144 if (fMeasuredMap) {
145 AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
146 return kTRUE;
147 }
ff66b122 148 //
db83d72f 149 char* fname = gSystem->ExpandPathName(GetDataFileName());
150 TFile* file = TFile::Open(fname);
151 if (!file) {
152 AliError(Form("Failed to open magnetic field data file %s\n",fname));
153 return kFALSE;
154 }
ff66b122 155 //
db83d72f 156 fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
157 if (!fMeasuredMap) {
158 AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
159 return kFALSE;
160 }
161 file->Close();
162 delete file;
163 return kTRUE;
ff66b122 164}
165
db83d72f 166
ff66b122 167//_______________________________________________________________________
db83d72f 168void AliMagF::Field(const Double_t *xyz, Double_t *b)
fe4da5cc 169{
db83d72f 170 // Method to calculate the field at point xyz
aee8290b 171 //
9251fceb 172 // b[0]=b[1]=b[2]=0.0;
173 if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
db83d72f 174 fMeasuredMap->Field(xyz,b);
175 if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
9251fceb 176 else for (int i=3;i--;) b[i] *= fFactorDip;
db83d72f 177 }
9251fceb 178 else MachineField(xyz, b);
aee8290b 179 //
fe4da5cc 180}
eeda4611 181
182//_______________________________________________________________________
db83d72f 183Double_t AliMagF::GetBz(const Double_t *xyz) const
eeda4611 184{
db83d72f 185 // Method to calculate the field at point xyz
186 //
9251fceb 187 if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
188 double bz = fMeasuredMap->GetBz(xyz);
189 return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;
db83d72f 190 }
9251fceb 191 else return 0.;
eeda4611 192}
193
194//_______________________________________________________________________
db83d72f 195AliMagF& AliMagF::operator=(const AliMagF& src)
eeda4611 196{
db83d72f 197 if (this != &src && src.fMeasuredMap) {
198 if (fMeasuredMap) delete fMeasuredMap;
199 fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
200 SetName(src.GetName());
201 fSolenoid = src.fSolenoid;
202 fBeamType = src.fBeamType;
203 fBeamEnergy = src.fBeamEnergy;
db83d72f 204 fInteg = src.fInteg;
205 fPrecInteg = src.fPrecInteg;
206 fFactorSol = src.fFactorSol;
207 fFactorDip = src.fFactorDip;
208 fMax = src.fMax;
209 fDipoleOFF = src.fDipoleOFF;
210 fParNames = src.fParNames;
211 }
212 return *this;
eeda4611 213}
214
215//_______________________________________________________________________
db83d72f 216void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
eeda4611 217{
9251fceb 218 if (btype==kNoBeamField || benergy<1.) {
db83d72f 219 fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
9251fceb 220 return;
db83d72f 221 }
222 //
9251fceb 223 double rigScale = benergy/7000.; // scale according to ratio of E/Enominal
224 // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
225 if (btype == kBeamTypeAA) rigScale *= 208./82.;
226 //
227 fQuadGradient = 22.0002*rigScale;
228 fDipoleField = 37.8781*rigScale;
229 //
230 // SIDE C
231 fCCorrField = -9.6980;
232 // SIDE A
233 fACorr1Field = -13.2247;
234 fACorr2Field = 11.7905;
db83d72f 235 //
eeda4611 236}
eed8a1a2 237
db83d72f 238//_______________________________________________________________________
239void AliMagF::MachineField(const Double_t *x, Double_t *b) const
eed8a1a2 240{
db83d72f 241 // ---- This is the ZDC part
9251fceb 242 // Compansators for Alice Muon Arm Dipole
243 const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
244 const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
245 //
246 const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
90ae20c9 247 const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
9251fceb 248 const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
249 const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
db83d72f 250 //
9251fceb 251 const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
252 const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
253 const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
db83d72f 254 //
255 double rad2 = x[0] * x[0] + x[1] * x[1];
256 //
9251fceb 257 b[0] = b[1] = b[2] = 0;
258 //
db83d72f 259 // SIDE C **************************************************
260 if(x[2]<0.){
9251fceb 261 if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
262 b[0] = fCCorrField*fFactorDip;
db83d72f 263 }
9251fceb 264 else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f 265 b[0] = fQuadGradient*x[1];
266 b[1] = fQuadGradient*x[0];
db83d72f 267 }
9251fceb 268 else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
db83d72f 269 b[0] = -fQuadGradient*x[1];
270 b[1] = -fQuadGradient*x[0];
db83d72f 271 }
9251fceb 272 else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
db83d72f 273 b[0] = -fQuadGradient*x[1];
274 b[1] = -fQuadGradient*x[0];
db83d72f 275 }
9251fceb 276 else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f 277 b[0] = fQuadGradient*x[1];
278 b[1] = fQuadGradient*x[0];
db83d72f 279 }
9251fceb 280 else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
db83d72f 281 b[1] = fDipoleField;
db83d72f 282 }
9251fceb 283 else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
284 double dxabs = TMath::Abs(x[0])-kDip2DXC;
285 if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
db83d72f 286 b[1] = -fDipoleField;
db83d72f 287 }
288 }
289 }
290 //
291 // SIDE A **************************************************
292 else{
9251fceb 293 if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
db83d72f 294 // Compensator magnet at z = 1075 m
9251fceb 295 b[0] = fACorr1Field*fFactorDip;
db83d72f 296 }
297 //
9251fceb 298 if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
299 b[0] = fACorr2Field*fFactorDip;
300 }
301 else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f 302 b[0] = -fQuadGradient*x[1];
303 b[1] = -fQuadGradient*x[0];
eed8a1a2 304 }
9251fceb 305 else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
306 b[0] = fQuadGradient*x[1];
307 b[1] = fQuadGradient*x[0];
eed8a1a2 308 }
9251fceb 309 else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
310 b[0] = fQuadGradient*x[1];
311 b[1] = fQuadGradient*x[0];
db83d72f 312 }
9251fceb 313 else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f 314 b[0] = -fQuadGradient*x[1];
315 b[1] = -fQuadGradient*x[0];
db83d72f 316 }
9251fceb 317 else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
db83d72f 318 b[1] = -fDipoleField;
db83d72f 319 }
9251fceb 320 else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
321 double dxabs = TMath::Abs(x[0])-kDip2DXA;
322 if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
db83d72f 323 b[1] = fDipoleField;
324 }
325 }
326 }
9251fceb 327 //
db83d72f 328}
329
330//_______________________________________________________________________
331void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
332{
333 // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
334 b[0]=b[1]=b[2]=0.0;
335 if (fMeasuredMap) {
336 fMeasuredMap->GetTPCInt(xyz,b);
337 for (int i=3;i--;) b[i] *= fFactorSol;
338 }
339}
340
341//_______________________________________________________________________
342void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
343{
344 // Method to calculate the integral of magnetic integral from point to nearest cathode plane
345 // in cylindrical coordiates ( -pi<phi<pi convention )
346 b[0]=b[1]=b[2]=0.0;
347 if (fMeasuredMap) {
348 fMeasuredMap->GetTPCIntCyl(rphiz,b);
349 for (int i=3;i--;) b[i] *= fFactorSol;
350 }
eed8a1a2 351}