EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PndMultiField.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PndMultiField.cxx
1 // -------------------------------------------------------------------------
2 // ----- PndMultiField source file -----
3 // ----- Created 29/01/07 by M. Al/Turany -----
4 // -------------------------------------------------------------------------
5 
6 
7 
8 #include "PndMultiField.h"
9 #include "PndRegion.h"
10 #include "PndConstField.h"
11 #include "PndFieldMap.h"
12 #include "PndMapPar.h"
13 #include "PndMultiFieldPar.h"
14 #include "PndSolenoidMap.h"
15 #include "PndTransMap.h"
16 #include "PndDipoleMap.h"
17 
18 
19 #include "FairRun.h"
20 #include "FairRuntimeDb.h"
21 
22 
23 #include "TObjArray.h"
24 
25 
26 
27 #include <iomanip>
28 #include <iostream>
29 #include <fstream>
30 
31 
32 
33 using namespace std;
34 
35 // ------------- Default constructor ----------------------------------
37  : fMaps(new TObjArray(10)), fNoOfMaps(0), fFieldMaps(), fMapIter(), fBeamMom(0.)
38 {
39  fType = 5;
40 }
41 
42 // ------------- Default constructor ----------------------------------
43 PndMultiField::PndMultiField(TString Map, Double_t BeamMom)
44  : fMaps(new TObjArray(10)), fNoOfMaps(0), fFieldMaps(), fMapIter(), fBeamMom(BeamMom)
45 {
46  fType = 5;
47 
48  Map.ToUpper();
49 
50  if (Map=="FULL") {
51  PndTransMap *map_t= new PndTransMap("TransMap", "R", fBeamMom);
52  PndDipoleMap *map_d1= new PndDipoleMap("DipoleMap1", "R", fBeamMom);
53  PndDipoleMap *map_d2= new PndDipoleMap("DipoleMap2", "R", fBeamMom);
54  PndSolenoidMap *map_s1= new PndSolenoidMap("SolenoidMap1", "R");
55  PndSolenoidMap *map_s2= new PndSolenoidMap("SolenoidMap2", "R");
56  PndSolenoidMap *map_s3= new PndSolenoidMap("SolenoidMap3", "R");
57  PndSolenoidMap *map_s4= new PndSolenoidMap("SolenoidMap4", "R");
58 
59  AddField(map_t);
60  AddField(map_d1);
61  AddField(map_d2);
62  AddField(map_s1);
63  AddField(map_s2);
64  AddField(map_s3);
65  AddField(map_s4);
66 
67  }else if (Map=="DIPOLE") {
68  PndDipoleMap *map_d1= new PndDipoleMap("DipoleMap1", "R", fBeamMom);
69  PndDipoleMap *map_d2= new PndDipoleMap("DipoleMap2", "R", fBeamMom);
70 
71  AddField(map_d1);
72  AddField(map_d2);
73 
74  }else if (Map=="SOLENOID") {
75 
76  PndSolenoidMap *map_s1= new PndSolenoidMap("SolenoidMap1", "R");
77  PndSolenoidMap *map_s2= new PndSolenoidMap("SolenoidMap2", "R");
78  PndSolenoidMap *map_s3= new PndSolenoidMap("SolenoidMap3", "R");
79  PndSolenoidMap *map_s4= new PndSolenoidMap("SolenoidMap4", "R");
80 
81  AddField(map_s1);
82  AddField(map_s2);
83  AddField(map_s3);
84  AddField(map_s4);
85  }
86 
87 
88 
89 }
90 
91 
92 
93 
94 // ------------------------------------------------------------------------
95 
96 // ------------ Constructor from PndFieldPar --------------------------
98  : fMaps( new TObjArray(10)), fNoOfMaps(0), fFieldMaps(), fMapIter(), fBeamMom(0.)
99 {
100  fType = 5;
101  TObjArray *fArray= fieldPar->GetParArray();
102  if(fArray->IsEmpty()) fType=-1;
103 
104 }
105 
106 // ------------ Destructor --------------------------------------------
108 
109  //printf("PndMultiField::~PndMultiField() \n");
110  fMaps->Delete();
111  delete fMaps;
112 
113 }
114 
115 // ----------- Adding fields ------------------------------------------
116 void PndMultiField::AddField(FairField *field){
117 
118  if(field){
119  fMaps->AddLast(field);
120  fNoOfMaps++;
121  }
122 
123 }
124 
125 // ----------- Intialisation ------------------------------------------
127  PndConstField *field=0;
128  PndFieldMap *fieldMap=0;
129  for (Int_t n=0; n<=fNoOfMaps; n++){
130  fieldMap = dynamic_cast<PndFieldMap *>(fMaps->At(n));
131  field = dynamic_cast<PndConstField *>(fMaps->At(n));
132  if(fieldMap){
133  fieldMap->Init();
134  fFieldMaps.insert( pair<PndRegion*, FairField*>(new PndRegion(fieldMap->GetZmin(),fieldMap->GetZmax()) , fieldMap ));
135  }else if(field){
136  field->Init();
137  fFieldMaps.insert( pair<PndRegion*, FairField*>(new PndRegion(field->GetZmin(),field->GetZmax() ), field ));
138  }
139  }
140 
141 }
142 
143 // --------- Screen output --------------------------------------------
145  for (Int_t n=0; n<=fNoOfMaps; n++){
146  FairField *fieldMap = dynamic_cast<FairField *>(fMaps->At(n));
147  if(fieldMap) fieldMap->Print();
148  }
149 }
150 
151 
152 
153 
154 
155 // --------- Fill the parameters --------------------------------------------
157 // for (Int_t n=0; n<=fNoOfMaps; n++){
158 // FairField *fieldMap = dynamic_cast<FairField *>(fMaps->At(n));
159 // if(fieldMap) fieldMap->FillParContainer();
160 // }
161  FairRun *fRun=FairRun::Instance();
162  FairRuntimeDb *rtdb=fRun->GetRuntimeDb();
163  Bool_t kParameterMerged=kTRUE;
164  PndMultiFieldPar* Par = (PndMultiFieldPar*) rtdb->getContainer("PndMultiFieldPar");
165  Par->SetParameters(this);
166  Par->setInputVersion(fRun->GetRunId(),1);
167  Par->setChanged();
168 
169 
170 
171 }
172 
173 // -------------------------------------------------------------------------
174 void PndMultiField::GetFieldValue(const Double_t point[3], Double_t* bField)
175 {
176  PndRegion *fReg=0;
177  FairField *fField=0;
178  for (fMapIter=fFieldMaps.begin(); fMapIter!= fFieldMaps.end();fMapIter++ ){
179  fReg=fMapIter->first;
180  if(fReg->IsInside(point[2])){
181  fField=fMapIter->second;
182  break;
183  }
184  }
185  if(fField){
186  //printf("%s\n", fField->GetName());
187  fField->GetBxyz(point, bField);
188  }else{
189  bField[0] = 0;
190  bField[1] = 0;
191  bField[2] = 0;
192 
193  }
194 
195  // 2015/08/28 -> NB: EicMagneticField does this job; this code has never been used;
196 #if 0
197  // Well, I rather want a superposition of all fields contributiong to
198  // this Z-range; in fact this crap is needed for DID test studies only;
199  for(unsigned iq=0; iq<3; iq++)
200  bField[iq] = 0.0;
201 
202  for (fMapIter=fFieldMaps.begin(); fMapIter!= fFieldMaps.end();fMapIter++ ){
203  PndRegion *fReg = fMapIter->first;
204  if(fReg->IsInside(point[2])){
205  FairField *fField = fMapIter->second;
206 
207  double buffer[3];
208  fField->GetBxyz(point, buffer);
209 
210  for(unsigned iq=0; iq<3; iq++)
211  bField[iq] += buffer[iq];
212  } //if
213  } //for
214 #endif
215 }
216 
217