HYDRA_development_version
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
hmdctrackgfield.cc
Go to the documentation of this file.
1 //*--Author A.Rustamov
2 //
3 //Input position vector,output field vector
4 //In the Sector coordinate system
5 //
6 // |Y
7 // |
8 // ***************
9 // * | *
10 // * | *
11 // * | *
12 // * | *
13 // * | *
14 // ***
15 // |
16 // -------------
17 // X
18 //
19 
20 using namespace std;
21 #include <iostream>
22 #include <iomanip>
23 #include <fstream>
24 #include "hmdctrackgfield.h"
25 #include "TTimer.h"
26 #include "TBuffer.h"
27 #include <math.h>
28 #include <unistd.h>
29 
30 ClassImp(HMdcTrackGField)
31 
32 HMdcTrackGField::HMdcTrackGField(const Char_t* name, const Char_t* title) :
33 TNamed(name,title)
34 {
35  dconv=0.0001;
36 
37  p_tzfl=new Double_t[708400]; //161*176*25
38  p_tpfl=new Double_t[708400];
39  p_trfl=new Double_t[708400];
40  nfz=161;
41  nfr=176;
42  nfp=25;
43  nfz_nfr=nfz*nfr;
44 
45  one_sixtyth=1.0/60.0;
46 
47 // printf("building acos table of size %d points...",ACOS_TABLE_SIZE);
48  acos_table=new Double_t[ACOS_TABLE_SIZE*2+1];
49  Double_t temp;
50  Double_t t1,t2;
51  t2=Double_t(ACOS_TABLE_SIZE);
52  Double_t radd=57.29577951;
53 //fill negative arguments:
54  for(Int_t i=0;i<ACOS_TABLE_SIZE;i++)
55  {
56  t1=Double_t(i-ACOS_TABLE_SIZE);
57  temp=t1/t2;
58  acos_table[i]=acos(temp)*radd;
59  }
60 
61 //fill zero argument:
62  acos_table[ACOS_TABLE_SIZE]=acos(0.0)*radd;
63 
64 //fill positive arguments:
65  for(Int_t i=ACOS_TABLE_SIZE+1;i<(ACOS_TABLE_SIZE*2+1);i++)
66  {
67  t1=Double_t(i-ACOS_TABLE_SIZE);
68  temp=t1/t2;
69  acos_table[i]=acos(temp)*radd;
70  }
71 // printf("Ok\n");
72 
73 }
74 
75 HMdcTrackGField::~HMdcTrackGField()
76 {
77  delete [] p_tzfl;
78  delete [] p_trfl;
79  delete [] p_tpfl;
80  delete [] acos_table;
81 }
82 void HMdcTrackGField::clear()
83 {
84 }
85 
86 void HMdcTrackGField::init(TString infile)
87 {
88  ifstream input;
89  input.open(infile.Data(),ios::in);
90  while(!input.eof())
91  {
92 
93  input>>zflmin ;
94  input>>zflmax ;
95  input>>zfldel;
96  input>>rflmin ;
97  input>>rflmax ;
98  input>>rfldel ;
99  input>>pflmin ;
100  input>>pflmax ;
101  input>>pfldel ;
102  for (Int_t i1 = 0; i1 < nfp; i1 ++){
103  for (Int_t i2 = 0; i2 < nfr; i2 ++){
104  for (Int_t i3 = 0;i3 < nfz; i3 ++){input>>p_tzfl[i3+i2*nfz+i1*nfz*nfr];}
105  for (Int_t i3 = 0;i3 < nfz; i3 ++){input>>p_trfl[i3+i2*nfz+i1*nfz*nfr];}
106  for (Int_t i3 = 0;i3 < nfz; i3 ++){input>>p_tpfl[i3+i2*nfz+i1*nfz*nfr];}
107  }
108  }
109  }
110  cout<<"READING FIELDMAP FINISHED!"<<endl;
111  input.close();
112  step1z=1.0/zfldel;
113  step1r=1.0/rfldel;
114  step1p=1.0/pfldel;
115 }
116 
117 void HMdcTrackGField::calcField(Double_t *xv,Double_t *btos,Double_t fpol )
118 {
119  Double_t xloc,yloc,zloc,rhog,phigd,phil,bz,br,bp;
120  Double_t m1,m2,m3;
121  Double_t dc,ds,a,b,c,ac,ab,bc,abc,w1,w2,w3,w4,w5,w6,w7,w8;
122  Double_t spol;
123  Double_t mult;
124  Int_t ifz,ifr,ifp,ifr1,ifp1;
125  Int_t i1,i2,i3,i4,i5,i6,i7,i8;
126 
127  zloc=xv[2];
128 
129  if(zloc<zflmin||zloc>zflmax) {
130  btos[0]=0.; btos[1]=0.;btos[2]=0.;
131 // printf("WARN!! z returned,xv=%f %f %f, zloc=%f \n",xv[0],xv[1],xv[2]);
132  return;
133  }
134  xloc=xv[0];
135  yloc=xv[1];
136  rhog=sqrt(xloc*xloc+yloc*yloc);
137  if(rhog>rflmax){
138  btos[0]=0.;
139  btos[1]=0.;
140  btos[2]=0.;
141 // printf("WARN!!! r returned xv=%f %f %f rhog=%f\n",xv[0],xv[1],xv[2]);
142  return;
143  }
144  if(rhog>0.01){
145  dc=xloc/rhog;
146  ds=yloc/rhog;
147  phigd=acos_table[Int_t(dc*ACOS_TABLE_SIZE+ACOS_TABLE_SIZE)];
148 
149  }else{
150  rhog=0.;
151  phigd=0.;
152  dc=0.;
153  ds=0.;
154 // printf("WARN!!!! rhog<=eps!!!\n");
155  }
156 
157  phil=phigd-(Int_t)(phigd*one_sixtyth)*60.0;
158 //phil=phigd-60.0;
159 
160  if(phil>30.0){
161  spol=-1.0;
162  phil=60.0-phil;
163  }else {
164  spol=1.0;
165  }
166 
167 
168  m1=(zloc-zflmin)*step1z;
169  m2=(rhog-rflmin)*step1r;
170  m3=(phil-pflmin)*step1p;
171 
172  ifz=(Int_t)(m1);
173  ifr=(Int_t)(m2);
174  ifp=(Int_t)(m3);
175 
176  a=m1-ifz;
177  b=m2-ifr;
178  c=m3-ifp;
179 
180  ifz++;
181  ifr=(ifr+1)*nfz;
182  ifp=(ifp+1)*nfz_nfr;
183  ifr1=ifr-nfz;
184  ifp1=ifp-nfz_nfr;
185 
186  ab=a*b;
187  ac=a*c;
188  bc=b*c;
189  abc=a*bc;
190  w1=1.0-a-b+ab-c+ac+bc-abc;
191  w2=b-ab-bc+abc;
192  w3=bc-abc;
193  w4=c-ac-bc+abc;
194  w5=a-ab-ac+abc;
195  w6=ab-abc;
196  w7=abc;
197  w8=ac-abc;
198 
199  i1=ifz-1+ifr1+ifp1;
200  i2=ifz-1+ifr+ifp1;
201  i3=ifz-1+ifr+ifp;
202  i4=ifz-1+ifr1+ifp;
203  i5=ifz+ifr1+ifp1;
204  i6=ifz+ifr+ifp1;
205  i7=ifz+ifr+ifp;
206  i8=ifz+ifr1+ifp;
207 
208  bz=w1*p_tzfl[i1]+w2*p_tzfl[i2]
209  +w3*p_tzfl[i3]+w4*p_tzfl[i4]
210  +w5*p_tzfl[i5]+w6*p_tzfl[i6]
211  +w7*p_tzfl[i7]+w8*p_tzfl[i8];
212 
213  br=w1*p_trfl[i1]+w2*p_trfl[i2]
214  +w3*p_trfl[i3]+w4*p_trfl[i4]
215  +w5*p_trfl[i5]+w6*p_trfl[i6]
216  +w7*p_trfl[i7]+w8*p_trfl[i8];
217 
218  bp=w1*p_tpfl[i1]+w2*p_tpfl[i2]
219  +w3*p_tpfl[i3]+w4*p_tpfl[i4]
220  +w5*p_tpfl[i5]+w6*p_tpfl[i6]
221  +w7*p_tpfl[i7]+w8*p_tpfl[i8];
222 
223 
224 
225  mult=fpol*dconv;
226  bz=-bz*spol*mult;
227  br=-br*spol*mult;
228  bp=-bp*mult;
229 
230  btos[0]=(dc*br-ds*bp);
231  btos[1]=(ds*br+dc*bp);
232  btos[2]=bz;
233 // printf("xv=%f %f %f, btos=%f %f %f\n",xv[0],xv[1],xv[2],btos[0],btos[1],btos[2]);
234 }
235 
236 
237 void HMdcTrackGField::Streamer(TBuffer &R__b) {
238  // Stream an object of class HMdcTrackGField.
239  // The data in the 1-dimensional arrays are streamed into/from the old
240  // 3-dimensional arrays (no change of class version)
241  UInt_t R__s, R__c;
242  if (R__b.IsReading()) {
243  Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { }
244  TNamed::Streamer(R__b);
245  R__b >> nfz;
246  R__b >> nfr;
247  R__b >> nfp;
248  R__b >> zflmin;
249  R__b >> zflmax;
250  R__b >> zfldel;
251  R__b >> rflmin;
252  R__b >> rflmax;
253  R__b >> rfldel;
254  R__b >> pflmin;
255  R__b >> pflmax;
256  R__b >> pfldel;
257 
258  Int_t n;
259  R__b >> n;
260  for (Int_t i3 = 0; i3 < nfz; i3 ++)
261  for (Int_t i2 = 0; i2 < nfr; i2 ++)
262  for (Int_t i1 = 0; i1 < nfp; i1 ++)
263  R__b >> p_tzfl[i3+i2*nfz+i1*nfz*nfr];
264 
265  R__b >> n;
266  for (Int_t i3 = 0; i3 < nfz; i3 ++)
267  for (Int_t i2 = 0; i2 < nfr; i2 ++)
268  for (Int_t i1 = 0; i1 < nfp; i1 ++)
269  R__b >> p_trfl[i3+i2*nfz+i1*nfz*nfr];
270 
271  R__b >> n;
272  for (Int_t i3 = 0; i3 < nfz; i3 ++)
273  for (Int_t i2 = 0; i2 < nfr; i2 ++)
274  for (Int_t i1 = 0; i1 < nfp; i1 ++)
275  R__b >> p_tpfl[i3+i2*nfz+i1*nfz*nfr];
276 
277  step1z=1.0/zfldel;
278  step1r=1.0/rfldel;
279  step1p=1.0/pfldel;
280 
281  R__b.ReadStaticArray((Double_t*)Pvector);
282  R__b.ReadStaticArray((Double_t*)Fvector);
283 
284  R__b.CheckByteCount(R__s, R__c, HMdcTrackGField::IsA());
285 
286  } else {
287  R__c = R__b.WriteVersion(HMdcTrackGField::IsA(), kTRUE);
288  TNamed::Streamer(R__b);
289  R__b << nfz;
290  R__b << nfr;
291  R__b << nfp;
292  R__b << zflmin;
293  R__b << zflmax;
294  R__b << zfldel;
295  R__b << rflmin;
296  R__b << rflmax;
297  R__b << rfldel;
298  R__b << pflmin;
299  R__b << pflmax;
300  R__b << pfldel;
301 
302  Int_t n=nfz*nfr*nfp;
303  R__b << n;
304  for (Int_t i3 = 0; i3 < nfz; i3 ++)
305  for (Int_t i2 = 0; i2 < nfr; i2 ++)
306  for (Int_t i1 = 0; i1 < nfp; i1 ++)
307  R__b << p_tzfl[i3+i2*nfz+i1*nfz*nfr];
308 
309  R__b << n;
310  for (Int_t i3 = 0; i3 < nfz; i3 ++)
311  for (Int_t i2 = 0; i2 < nfr; i2 ++)
312  for (Int_t i1 = 0; i1 < nfp; i1 ++)
313  R__b << p_trfl[i3+i2*nfz+i1*nfz*nfr];
314 
315  R__b << n;
316  for (Int_t i3 = 0; i3 < nfz; i3 ++)
317  for (Int_t i2 = 0; i2 < nfr; i2 ++)
318  for (Int_t i1 = 0; i1 < nfp; i1 ++)
319  R__b << p_tpfl[i3+i2*nfz+i1*nfz*nfr];
320 
321  R__b.WriteArray(Pvector, 3);
322  R__b.WriteArray(Fvector, 3);
323  R__b.SetByteCount(R__c, kTRUE);
324  }
325 }
326 
HMdcTrackGField::calcField
void calcField(Double_t *xv, Double_t *btos, Double_t fpol)
Definition: hmdctrackgfield.cc:117
n
Int_t n
Definition: drawResolution.C:34
m2
Int_t m2
Definition: drawAccCuts.C:11
m1
Int_t m1
Definition: drawAccCuts.C:11
HMdcTrackGField::init
void init(TString infile)
Definition: hmdctrackgfield.cc:86
input
TString input
Definition: GarReader.C:5
ClassImp
ClassImp(HMdcTrackGField) HMdcTrackGField
Definition: hmdctrackgfield.cc:30
HMdcTrackGField::~HMdcTrackGField
virtual ~HMdcTrackGField()
Definition: hmdctrackgfield.cc:75
ACOS_TABLE_SIZE
#define ACOS_TABLE_SIZE
Definition: hmdctrackgfield.h:26
m3
Int_t m3
Definition: drawAccCuts.C:11