]>
Commit | Line | Data |
---|---|---|
d92ed900 | 1 | /////////////////////////////////////////////////////////////////////////// |
2 | // // | |
3 | // AliFemtoCoulomb: This is a Coulomb correction class which // | |
4 | // 1. Reads in the dat from a file // | |
5 | // 2. Performs a linear interpolation in R and creates any array of // | |
6 | // interpolations // | |
7 | // 3. Interpolates in eta and returns the Coulomb correction to user // | |
8 | // // | |
9 | /////////////////////////////////////////////////////////////////////////// | |
67427ff7 | 10 | |
11 | #include "AliFemtoCoulomb.h" | |
12 | //#include "Stiostream.h" | |
13 | #include <stdio.h> | |
14 | #include <cassert> | |
556650dc | 15 | //#include "PhysicalConstants.h" |
16 | #define fine_structure_const 0.00729735 | |
67427ff7 | 17 | |
18 | #ifdef __ROOT__ | |
19 | ClassImp(AliFemtoCoulomb) | |
20 | #endif | |
21 | ||
0215f606 | 22 | AliFemtoCoulomb::AliFemtoCoulomb() : |
23 | fFile(""), | |
24 | fRadius(-1.0), | |
25 | fZ1Z2(1.0), | |
26 | fNLines(0) | |
27 | { | |
d92ed900 | 28 | // Default constructor |
67427ff7 | 29 | fFile = "/afs/rhic/star/hbt/coul/AliFemtoCorrectionFiles/correctionpp.dat"; |
30 | if (!fFile) { | |
31 | cout << " No file, dummy!" << endl; | |
32 | assert(0); | |
33 | } | |
67427ff7 | 34 | cout << "You have 1 default Coulomb correction!" << endl; |
14c25543 | 35 | |
36 | for (int ie=0; ie<1000; ie++) { | |
37 | fEta[ie] = 0.0; | |
38 | fCoulomb[ie] = 0.0; | |
39 | } | |
67427ff7 | 40 | } |
41 | ||
0215f606 | 42 | AliFemtoCoulomb::AliFemtoCoulomb(const AliFemtoCoulomb& aCoul) : |
43 | fFile(aCoul.fFile), | |
44 | fRadius(aCoul.fRadius), | |
45 | fZ1Z2(aCoul.fZ1Z2), | |
46 | fNLines(0) | |
47 | { | |
d92ed900 | 48 | // copy constructor |
0215f606 | 49 | CreateLookupTable(fRadius); |
50 | } | |
51 | ||
52 | AliFemtoCoulomb::AliFemtoCoulomb(const char* readFile, const double& radius, const double& charge) : | |
53 | fFile(readFile), | |
54 | fRadius(radius), | |
55 | fZ1Z2(0), | |
56 | fNLines(0) | |
57 | { | |
d92ed900 | 58 | // constructor with explicit filename |
67427ff7 | 59 | fFile = readFile; |
60 | fRadius = radius; | |
61 | CreateLookupTable(fRadius); | |
62 | fZ1Z2 = charge; | |
63 | cout << "You have 1 Coulomb correction!" << endl; | |
64 | } | |
65 | ||
66 | AliFemtoCoulomb::~AliFemtoCoulomb() { | |
d92ed900 | 67 | // destructor |
67427ff7 | 68 | } |
69 | ||
0215f606 | 70 | AliFemtoCoulomb& AliFemtoCoulomb::operator=(const AliFemtoCoulomb& aCoul) |
71 | { | |
d92ed900 | 72 | // assignment operator |
0215f606 | 73 | if (this == &aCoul) |
74 | return *this; | |
75 | ||
76 | fFile = aCoul.fFile; | |
77 | fRadius = aCoul.fRadius; | |
78 | fZ1Z2 = aCoul.fZ1Z2; | |
79 | ||
80 | CreateLookupTable(fRadius); | |
81 | ||
82 | return *this; | |
83 | } | |
84 | ||
85 | ||
67427ff7 | 86 | void AliFemtoCoulomb::SetRadius(const double& radius) { |
d92ed900 | 87 | // set the coulomb radius |
67427ff7 | 88 | cout << " AliFemtoCoulomb::setRadius() " << endl; |
89 | fRadius = radius; | |
90 | CreateLookupTable(fRadius); | |
91 | } | |
92 | ||
d92ed900 | 93 | double AliFemtoCoulomb::GetRadius() const { |
94 | // return coulomb radius | |
67427ff7 | 95 | return (fRadius); |
96 | } | |
97 | ||
98 | void AliFemtoCoulomb::SetFile(const char* readFile) { | |
d92ed900 | 99 | // set the filename with coulomb calculations |
67427ff7 | 100 | cout << " AliFemtoCoulomb::SetFile() " << endl; |
101 | fFile = readFile; | |
102 | // Create new lookup table since file has changed | |
103 | if (fRadius>0.0) { | |
104 | CreateLookupTable(fRadius); | |
105 | } | |
106 | } | |
107 | ||
108 | void AliFemtoCoulomb::SetChargeProduct(const double& charge) { | |
d92ed900 | 109 | // set pair charge |
67427ff7 | 110 | cout << " AliFemtoCoulomb::SetChargeProduct() " << endl; |
111 | if ( fZ1Z2!=charge ) { | |
112 | fZ1Z2 = charge; | |
113 | if ( fZ1Z2>0 ) { | |
114 | fFile = "/afs/rhic/star/hbt/coul/AliFemtoCorrectionFiles/correctionpp.dat"; | |
115 | } | |
116 | else { | |
117 | fFile = "/afs/rhic/star/hbt/coul/AliFemtoCorrectionFiles/correctionpm.dat"; | |
118 | } | |
119 | CreateLookupTable(fRadius); | |
120 | } | |
121 | } | |
122 | ||
123 | void AliFemtoCoulomb::CreateLookupTable(const double& radius) { | |
67427ff7 | 124 | // Read radii from fFile |
125 | // Create array(pair) of linear interpolation between radii | |
d92ed900 | 126 | cout << " AliFemtoCoulomb::CreateLookupTable() " << endl; |
67427ff7 | 127 | |
128 | if (radius<0.0) { | |
129 | cout << " AliFemtoCoulomb::CreateLookupTable -> NEGATIVE RADIUS " << endl; | |
130 | cout << " call AliFemtoCoulomb::SetRadius(r) with positive r " << endl; | |
131 | cerr << " AliFemtoCoulomb::CreateLookupTable -> NEGATIVE RADIUS " << endl; | |
132 | cerr << " call AliFemtoCoulomb::SetRadius(r) with positive r " << endl; | |
133 | assert(0); | |
134 | } | |
135 | ifstream mystream(fFile); | |
136 | if (!mystream) { | |
137 | cout << "Could not open file" << endl; | |
138 | assert(0); | |
139 | } | |
140 | else { | |
141 | cout << "Input correction file opened" << endl; | |
142 | } | |
143 | ||
144 | static char tempstring[2001]; | |
145 | static float radii[2000]; | |
d92ed900 | 146 | static int tNRadii = 0; |
147 | tNRadii = 0; | |
67427ff7 | 148 | if (!mystream.getline(tempstring,2000)) { |
149 | cout << "Could not read radii from file" << endl; | |
150 | assert(0); | |
151 | } | |
152 | for (unsigned int ii=0; ii<strlen(tempstring); ii++) { | |
153 | while (tempstring[ii]==' ') ii++; | |
d92ed900 | 154 | sscanf(&tempstring[ii++],"%f",&radii[++tNRadii]); |
67427ff7 | 155 | while ( tempstring[ii]!=' ' && (ii)<strlen(tempstring) )ii++; |
156 | } | |
d92ed900 | 157 | cout << " Read " << tNRadii << " radii from file" << endl; |
158 | ||
159 | static double tLowRadius = -1.0; | |
160 | static double tHighRadius = -1.0; | |
161 | static int tLowIndex = 0; | |
162 | tLowRadius = -1.0; | |
163 | tHighRadius = -1.0; | |
164 | tLowIndex = 0; | |
165 | for(int iii=1; iii<=tNRadii-1; iii++) { // Loop to one less than #radii | |
67427ff7 | 166 | if ( radius >= radii[iii] && radius <= radii[iii+1] ) { |
d92ed900 | 167 | tLowRadius = radii[iii]; |
168 | tHighRadius = radii[iii+1]; | |
169 | tLowIndex = iii; | |
67427ff7 | 170 | } |
171 | } | |
d92ed900 | 172 | if ( (tLowRadius < 0.0) || (tHighRadius < 0.0) ) { |
67427ff7 | 173 | cout << "AliFemtoCoulomb::CreateLookupTable --> Problem interpolating radius" << endl; |
174 | cout << " Check range of radii in lookup file...." << endl; | |
175 | cerr << "AliFemtoCoulomb::CreateLookupTable --> Problem interpolating radius" << endl; | |
176 | cerr << " Check range of radii in lookup file...." << endl; | |
177 | assert(0); | |
178 | } | |
179 | ||
d92ed900 | 180 | static double corr[100]; // array of corrections ... must be > tNRadii |
67427ff7 | 181 | fNLines = 0; |
182 | static double tempEta = 0; | |
183 | tempEta = 0; | |
184 | while (mystream >> tempEta) { | |
d92ed900 | 185 | for (int i=1; i<=tNRadii; i++) { |
67427ff7 | 186 | mystream >> corr[i]; |
187 | } | |
d92ed900 | 188 | static double tLowCoulomb = 0; |
189 | static double tHighCoulomb = 0; | |
67427ff7 | 190 | static double nCorr = 0; |
d92ed900 | 191 | tLowCoulomb = corr[tLowIndex]; |
192 | tHighCoulomb = corr[tLowIndex+1]; | |
193 | nCorr = ( (radius-tLowRadius)*tHighCoulomb+(tHighRadius-radius)*tLowCoulomb )/(tHighRadius-tLowRadius); | |
67427ff7 | 194 | fEta[fNLines] = tempEta; // Eta |
195 | fCoulomb[fNLines] = nCorr; // Interpolated Coulomb correction for radius | |
196 | fNLines++; | |
197 | } | |
198 | mystream.close(); | |
199 | cout << "Lookup Table is created with " << fNLines << " points" << endl; | |
200 | } | |
201 | ||
202 | double AliFemtoCoulomb::CoulombCorrect(const double& eta) { | |
203 | // Interpolates in eta | |
204 | if (fRadius < 0.0) { | |
205 | cout << "AliFemtoCoulomb::CoulombCorrect(eta) --> Trying to correct for negative radius!" << endl; | |
206 | cerr << "AliFemtoCoulomb::CoulombCorrect(eta) --> Trying to correct for negative radius!" << endl; | |
207 | assert(0); | |
208 | } | |
209 | static int middle=0; | |
210 | middle=int( (fNLines-1)/2 ); | |
211 | if (eta*fEta[middle]<0.0) { | |
212 | cout << "AliFemtoCoulomb::CoulombCorrect(eta) --> eta: " << eta << " has wrong sign for data file! " << endl; | |
213 | cerr << "AliFemtoCoulomb::CoulombCorrect(eta) --> eta: " << eta << " has wrong sign for data file! " << endl; | |
214 | assert(0); | |
215 | } | |
216 | ||
d92ed900 | 217 | static double tCorr = 0; |
218 | tCorr = -1.0; | |
67427ff7 | 219 | |
220 | if ( (eta>fEta[0]) && (fEta[0]>0.0) ) { | |
d92ed900 | 221 | tCorr = fCoulomb[0]; |
222 | return (tCorr); | |
67427ff7 | 223 | } |
224 | if ( (eta<fEta[fNLines-1]) && (fEta[fNLines-1]<0.0) ) { | |
d92ed900 | 225 | tCorr = fCoulomb[fNLines-1]; |
226 | return (tCorr); | |
67427ff7 | 227 | } |
228 | // This is a binary search for the bracketing pair of data points | |
229 | static int high = 0; | |
230 | static int low = 0; | |
231 | static int width = 0; | |
232 | high = fNLines-1; | |
233 | low = 0; | |
234 | width = high-low; | |
235 | middle = int(width/2.0); // Was instantiated above | |
236 | while (middle > 0) { | |
237 | if (fEta[low+middle] < eta) { | |
238 | // eta is in the 1st half | |
239 | high-=middle; | |
240 | width = high-low; | |
241 | middle = int(width/2.0); | |
242 | } | |
243 | else { | |
244 | // eta is in the 2nd half | |
245 | low+=middle; | |
246 | width = high-low; | |
247 | middle = int(width/2.0); | |
248 | } | |
249 | } | |
250 | // Make sure we found the right one | |
251 | if ( (fEta[low] >= eta) && (eta >= fEta[low+1]) ) { | |
d92ed900 | 252 | static double tLowEta = 0; |
253 | static double tHighEta = 0; | |
254 | static double tLowCoulomb = 0; | |
255 | static double tHighCoulomb = 0; | |
256 | tLowEta = fEta[low]; | |
257 | tHighEta = fEta[low+1]; | |
258 | tLowCoulomb = fCoulomb[low]; | |
259 | tHighCoulomb = fCoulomb[low+1]; | |
260 | // cout << tLowEta << " *** Eta *** " << tHighEta << endl; | |
261 | // cout << tLowCoulomb << " *** Coulomb *** " << tHighCoulomb << endl; | |
262 | tCorr = ( (eta-tLowEta)*tHighCoulomb+(tHighEta-eta)*tLowCoulomb )/(tHighEta-tLowEta); | |
67427ff7 | 263 | } |
d92ed900 | 264 | if (tCorr<0.0) { |
67427ff7 | 265 | cout << "AliFemtoCoulomb::CoulombCorrect(eta) --> No correction" << endl; |
266 | cout << " Check range of eta in file: Input eta " << eta << endl; | |
267 | cerr << "AliFemtoCoulomb::CoulombCorrect(eta) --> No correction" << endl; | |
268 | cerr << " Check range of eta in file: Input eta " << eta << endl; | |
269 | assert(0); | |
270 | } | |
d92ed900 | 271 | return (tCorr); |
67427ff7 | 272 | |
273 | } | |
274 | ||
275 | double AliFemtoCoulomb::CoulombCorrect(const double& eta, | |
276 | const double& radius) { | |
277 | // Checks radii ... input radius and fRadius | |
278 | // Calls createLookupTable if neccessary | |
279 | // Interpolate(linear) between etas in the created lookup table | |
280 | ||
281 | if (radius < 0.0) { | |
282 | if (fRadius < 0.0) { | |
283 | // Both radii are negative | |
284 | cout << "AliFemtoCoulomb::CoulombCorrect(eta,r) --> input and member radii are negative!" << endl; | |
285 | cerr << "AliFemtoCoulomb::CoulombCorrect(eta,r) --> input and member radii are negative!" << endl; | |
286 | assert(0); | |
287 | } | |
288 | } | |
289 | else { | |
290 | // radius > 0.0 | |
291 | if (radius == fRadius) { | |
292 | // Both radii are positive and equal | |
293 | // cout << "Radii are the same!!!" << endl; | |
294 | } | |
295 | else { | |
296 | // Both radii are positive but not equal | |
297 | fRadius = radius; | |
298 | CreateLookupTable(fRadius); | |
299 | } | |
300 | } | |
301 | ||
302 | // Interpolate in eta | |
303 | return ( CoulombCorrect(eta) ); | |
304 | } | |
305 | ||
306 | double AliFemtoCoulomb::CoulombCorrect(const AliFemtoPair* pair) { | |
307 | return ( CoulombCorrect( Eta(pair) ) );; | |
308 | } | |
309 | ||
310 | double AliFemtoCoulomb::CoulombCorrect(const AliFemtoPair* pair, const double& radius) { | |
311 | return ( CoulombCorrect( Eta(pair),radius ) ); | |
312 | } | |
313 | ||
314 | double AliFemtoCoulomb::Eta(const AliFemtoPair* pair) { | |
d92ed900 | 315 | // calculate eta |
67427ff7 | 316 | static double px1,py1,pz1,px2,py2,pz2; |
317 | static double px1new,py1new,pz1new; | |
318 | static double px2new,py2new,pz2new; | |
319 | static double vx1cms,vy1cms,vz1cms; | |
320 | static double vx2cms,vy2cms,vz2cms; | |
d92ed900 | 321 | static double tVcmsX,tVcmsY,tVcmsZ; |
67427ff7 | 322 | static double dv = 0.0; |
323 | static double e1,e2,e1new,e2new; | |
324 | static double psi,theta; | |
325 | static double beta,gamma; | |
d92ed900 | 326 | static double tVcmsXnew; |
67427ff7 | 327 | |
d0e92d9a | 328 | px1 = pair->Track1()->FourMomentum().px(); |
329 | py1 = pair->Track1()->FourMomentum().py(); | |
330 | pz1 = pair->Track1()->FourMomentum().pz(); | |
331 | e1 = pair->Track1()->FourMomentum().e(); | |
332 | px2 = pair->Track2()->FourMomentum().px(); | |
333 | py2 = pair->Track2()->FourMomentum().py(); | |
334 | pz2 = pair->Track2()->FourMomentum().pz(); | |
335 | e2 = pair->Track2()->FourMomentum().e(); | |
67427ff7 | 336 | |
d92ed900 | 337 | tVcmsX = ( px1+px2 )/( e1+e2 ); |
338 | tVcmsY = ( py1+py2 )/( e1+e2 ); | |
339 | tVcmsZ = ( pz1+pz2 )/( e1+e2 ); | |
340 | // Rotate tVcms to x-direction | |
341 | psi = atan(tVcmsY/tVcmsX); | |
342 | tVcmsXnew = tVcmsX*cos(psi)+tVcmsY*sin(psi); | |
343 | tVcmsX = tVcmsXnew; | |
344 | theta = atan(tVcmsZ/tVcmsX); | |
345 | tVcmsXnew = tVcmsX*cos(theta)+tVcmsZ*sin(theta); | |
346 | tVcmsX = tVcmsXnew; | |
67427ff7 | 347 | // Gamma and Beta |
d92ed900 | 348 | beta = tVcmsX; |
67427ff7 | 349 | gamma = 1.0/::sqrt( 1.0-beta*beta ); |
350 | ||
351 | // Rotate p1 and p2 to new frame | |
352 | px1new = px1*cos(psi)+py1*sin(psi); | |
353 | py1new = -px1*sin(psi)+py1*cos(psi); | |
354 | px1 = px1new; | |
355 | px1new = px1*cos(theta)+pz1*sin(theta); | |
356 | pz1new = -px1*sin(theta)+pz1*cos(theta); | |
357 | px1 = px1new; | |
358 | py1 = py1new; | |
359 | pz1 = pz1new; | |
360 | ||
361 | px2new = px2*cos(psi)+py2*sin(psi); | |
362 | py2new = -px2*sin(psi)+py2*cos(psi); | |
363 | px2 = px2new; | |
364 | px2new = px2*cos(theta)+pz2*sin(theta); | |
365 | pz2new = -px2*sin(theta)+pz2*cos(theta); | |
366 | px2 = px2new; | |
367 | py2 = py2new; | |
368 | pz2 = pz2new; | |
369 | ||
370 | // Lorentz transform the x component and energy | |
371 | e1new = gamma*e1 - gamma*beta*px1; | |
372 | px1new = -gamma*beta*e1 + gamma*px1; | |
373 | e2new = gamma*e2 - gamma*beta*px2; | |
374 | px2new = -gamma*beta*e2 + gamma*px2; | |
375 | px1 = px1new; | |
376 | px2 = px2new; | |
377 | ||
378 | // New velocities | |
379 | vx1cms = px1/e1new; | |
380 | vy1cms = py1/e1new; | |
381 | vz1cms = pz1/e1new; | |
382 | vx2cms = px2/e2new; | |
383 | vy2cms = py2/e2new; | |
384 | vz2cms = pz2/e2new; | |
385 | ||
386 | // Velocity difference in CMS frame | |
387 | dv = ::sqrt( (vx1cms-vx2cms)*(vx1cms-vx2cms) + | |
388 | (vy1cms-vy2cms)*(vy1cms-vy2cms) + | |
389 | (vz1cms-vz2cms)*(vz1cms-vz2cms) ); | |
390 | ||
391 | return ( fZ1Z2*fine_structure_const/(dv) ); | |
392 | } | |
393 | ||
394 | TH1D* AliFemtoCoulomb::CorrectionHistogram(const double& mass1, const double& mass2, const int& nBins, | |
395 | const double& low, const double& high) { | |
d92ed900 | 396 | // return correction histogram |
397 | ||
67427ff7 | 398 | if ( mass1!=mass2 ) { |
399 | cout << "Masses not equal ... try again. No histogram created." << endl; | |
400 | assert(0); | |
401 | } | |
402 | TH1D* correction = new TH1D("correction","Coulomb correction",nBins,low,high); | |
d92ed900 | 403 | const double kReducedMass = mass1*mass2/(mass1+mass2); |
67427ff7 | 404 | double qInv = low; |
405 | //double dQinv = (high-low)/( (double)nBins ); | |
406 | double eta; | |
407 | for (int ii=1; ii<=nBins; ii++) | |
408 | { | |
409 | qInv = correction->GetBinCenter(ii); | |
d92ed900 | 410 | eta = 2.0*fZ1Z2*kReducedMass*fine_structure_const/( qInv ); |
67427ff7 | 411 | CoulombCorrect( eta ); |
412 | correction->Fill( qInv, CoulombCorrect(eta,fRadius) ); | |
413 | } | |
414 | ||
415 | return (correction); | |
416 | } | |
417 | ||
418 | #ifdef __ROOT__ | |
419 | TH1D* AliFemtoCoulomb::CorrectionHistogram(const TH1D* histo, const double mass) { | |
d92ed900 | 420 | // return correction histogram - 1D case |
67427ff7 | 421 | TH1D* correction = (TH1D*) ((TH1D*)histo)->Clone(); |
422 | correction->Reset(); | |
423 | correction->SetDirectory(0); | |
424 | int nBins = correction->GetXaxis()->GetNbins(); | |
d92ed900 | 425 | const double kReducedMass = 0.5*mass; |
67427ff7 | 426 | double qInv; |
427 | double eta; | |
428 | for (int ii=1; ii<=nBins; ii++) | |
429 | { | |
430 | qInv = correction->GetBinCenter(ii); | |
d92ed900 | 431 | eta = 2.0*fZ1Z2*kReducedMass*fine_structure_const/( qInv ); |
67427ff7 | 432 | correction->Fill( qInv, CoulombCorrect(eta,fRadius) ); |
433 | } | |
434 | ||
435 | return (correction); | |
436 | } | |
437 | ||
438 | TH3D* AliFemtoCoulomb::CorrectionHistogram(const TH3D* histo, const double mass) { | |
d92ed900 | 439 | // return correction histogram - 3D case |
67427ff7 | 440 | TH3D* correction = (TH3D*) ((TH3D*)histo)->Clone(); |
441 | correction->Reset(); | |
442 | correction->SetDirectory(0); | |
443 | int nBinsX = correction->GetXaxis()->GetNbins(); | |
444 | int nBinsY = correction->GetYaxis()->GetNbins(); | |
445 | int nBinsZ = correction->GetZaxis()->GetNbins(); | |
d92ed900 | 446 | const double kReducedMass = 0.5*mass; |
67427ff7 | 447 | double eta; |
448 | double qInv; | |
449 | int binNumber; | |
450 | for (int ii=1; ii<=nBinsX; ii++) { | |
451 | for (int iii=1; iii<=nBinsY; iii++) { | |
452 | for (int iv=1; iv<=nBinsZ; iv++) { | |
453 | binNumber = histo->GetBin(ii,iii,iv); | |
454 | qInv = histo->GetBinContent(binNumber); | |
d92ed900 | 455 | eta = 2.0*fZ1Z2*kReducedMass*fine_structure_const/( qInv ); |
67427ff7 | 456 | correction->SetBinContent(binNumber, CoulombCorrect(eta,fRadius) ); |
457 | } | |
458 | } | |
459 | } | |
460 | return (correction); | |
461 | } | |
462 | #endif | |
463 | ||
464 | double AliFemtoCoulomb::CoulombCorrect(const double& mass, const double& charge, | |
465 | const double& radius, const double& qInv) { | |
d92ed900 | 466 | // return correction factor |
67427ff7 | 467 | fRadius = radius; |
468 | fZ1Z2 = charge; | |
d92ed900 | 469 | const double kReducedMass = 0.5*mass; // must be same mass particles |
470 | double eta = 2.0*fZ1Z2*kReducedMass*fine_structure_const/( qInv ); | |
67427ff7 | 471 | return ( CoulombCorrect(eta,fRadius) ); |
472 | } |