]>
Commit | Line | Data |
---|---|---|
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 | ||
16 | /* $Id$ */ | |
17 | ||
18 | //---------------------------------------------------------------------------- | |
19 | // Implementation of the class to calculate the parton energy loss | |
20 | // Based on the "BDMPS" quenching weights by C.A.Salgado and U.A.Wiedemann | |
21 | // | |
22 | // References: | |
23 | // C.A.Salgado and U.A.Wiedemann, Phys.Rev.D68 (2003) 014008 [hep-ph/0302184] | |
24 | // A.Dainese, Eur.Phys.J.C, in press, [nucl-ex/0312005] | |
25 | // | |
26 | // | |
27 | // Origin: C. Loizides constantinos.loizides@cern.ch | |
28 | // A. Dainese andrea.dainese@pd.infn.it | |
29 | // | |
30 | //---------------------------------------------------------------------------- | |
31 | ||
32 | #include <Riostream.h> | |
33 | #include <TH1F.h> | |
34 | #include <TH2F.h> | |
35 | #include <TCanvas.h> | |
36 | #include <TGraph.h> | |
37 | #include <TROOT.h> | |
38 | #include <TSystem.h> | |
39 | #include <TLegend.h> | |
40 | #include "AliQuenchingWeights.h" | |
41 | ||
42 | ClassImp(AliQuenchingWeights) | |
43 | ||
44 | // conversion from fm to GeV^-1: 1 fm = fmGeV GeV^-1 | |
45 | const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197; | |
46 | ||
47 | // maximum value of R | |
48 | const Double_t AliQuenchingWeights::fgkRMax = 1.e6; | |
49 | ||
50 | // counter for histogram labels | |
51 | Int_t AliQuenchingWeights::fgCounter = 0; | |
52 | ||
53 | AliQuenchingWeights::AliQuenchingWeights() | |
54 | : TObject() | |
55 | { | |
56 | //default constructor | |
57 | ||
58 | fTablesLoaded=kFALSE; | |
59 | fMultSoft=kTRUE; | |
60 | fHistos=0; | |
61 | SetMu(); | |
62 | SetQTransport(); | |
63 | SetECMethod(); | |
64 | SetLengthMax(); | |
65 | fLengthMaxOld=0; | |
66 | fInstanceNumber=fgCounter++; | |
67 | } | |
68 | ||
69 | AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a) | |
70 | : TObject() | |
71 | { | |
72 | //copy constructor | |
73 | ||
74 | fTablesLoaded=kFALSE; | |
75 | fHistos=0; | |
76 | fLengthMaxOld=0; | |
77 | fMultSoft=a.GetMultSoft();; | |
78 | fMu=a.GetMu(); | |
79 | fQTransport=a.GetQTransport(); | |
80 | fECMethod=(kECMethod)a.GetECMethod(); | |
81 | fLengthMax=a.GetLengthMax(); | |
82 | fInstanceNumber=fgCounter++; | |
83 | //Missing in the class is the pathname | |
84 | //to the tables, can be added if needed | |
85 | } | |
86 | ||
87 | AliQuenchingWeights::~AliQuenchingWeights() | |
88 | { | |
89 | Reset(); | |
90 | } | |
91 | ||
92 | void AliQuenchingWeights::Reset() | |
93 | { | |
94 | //reset tables if there were used | |
95 | ||
96 | if(!fHistos) return; | |
97 | for(Int_t l=0;l<fLengthMaxOld;l++){ | |
98 | delete fHistos[0][l]; | |
99 | delete fHistos[1][l]; | |
100 | } | |
101 | delete[] fHistos; | |
102 | fHistos=0; | |
103 | fLengthMaxOld=0; | |
104 | } | |
105 | ||
106 | void AliQuenchingWeights::SetECMethod(kECMethod type) | |
107 | { | |
108 | //set energy constraint method | |
109 | ||
110 | fECMethod=type; | |
111 | if(fECMethod==kDefault) | |
112 | Info("SetECMethod","Energy Constraint Method set to DEFAULT:\nIf (sampled energy loss > parton energy) then sampled energy loss = parton energy."); | |
113 | else | |
114 | Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy."); | |
115 | } | |
116 | ||
117 | Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall) | |
118 | { | |
119 | // read in tables for multiple scattering approximation | |
120 | // path to continuum and to discrete part | |
121 | ||
122 | fTablesLoaded = kFALSE; | |
123 | fMultSoft=kTRUE; | |
124 | ||
125 | Char_t fname[1024]; | |
126 | sprintf(fname,"%s",gSystem->ExpandPathName(contall)); | |
127 | //PH ifstream fincont(fname); | |
128 | fstream fincont(fname,ios::in); | |
129 | #if defined(__HP_aCC) || defined(__DECCXX) | |
130 | if(!fincont.rdbuf()->is_open()) return -1; | |
131 | #else | |
132 | if(!fincont.is_open()) return -1; | |
133 | #endif | |
134 | ||
135 | Int_t nn=0; //quarks | |
136 | while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>> | |
137 | fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>> | |
138 | fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>fcaq[13][nn]>> | |
139 | fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>fcaq[18][nn]>> | |
140 | fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>fcaq[23][nn]>> | |
141 | fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>fcaq[28][nn]>> | |
142 | fcaq[29][nn]>>fcaq[30][nn]>>fcaq[31][nn]>>fcaq[32][nn]>>fcaq[33][nn]) | |
143 | { | |
144 | nn++; | |
145 | if(nn==261) break; | |
146 | } | |
147 | ||
148 | nn=0; //gluons | |
149 | while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>> | |
150 | fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>> | |
151 | fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>fcag[13][nn]>> | |
152 | fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>fcag[18][nn]>> | |
153 | fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>fcag[23][nn]>> | |
154 | fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>fcag[28][nn]>> | |
155 | fcag[29][nn]>>fcag[30][nn]>>fcag[31][nn]>>fcag[32][nn]>>fcag[33][nn]) | |
156 | { | |
157 | nn++; | |
158 | if(nn==261) break; | |
159 | } | |
160 | fincont.close(); | |
161 | ||
162 | sprintf(fname,"%s",gSystem->ExpandPathName(discall)); | |
163 | //PH ifstream findisc(fname); | |
164 | fstream findisc(fname,ios::in); | |
165 | #if defined(__HP_aCC) || defined(__DECCXX) | |
166 | if(!findisc.rdbuf()->is_open()) return -1; | |
167 | #else | |
168 | if(!findisc.is_open()) return -1; | |
169 | #endif | |
170 | ||
171 | nn=0; //quarks | |
172 | while(findisc>>frrr[nn]>>fdaq[nn]) { | |
173 | nn++; | |
174 | if(nn==34) break; | |
175 | } | |
176 | nn=0; //gluons | |
177 | while(findisc>>frrrg[nn]>>fdag[nn]) { | |
178 | nn++; | |
179 | if(nn==34) break; | |
180 | } | |
181 | findisc.close(); | |
182 | fTablesLoaded = kTRUE; | |
183 | return 0; | |
184 | } | |
185 | ||
186 | /* | |
187 | C*************************************************************************** | |
188 | C Quenching Weights for Multiple Soft Scattering | |
189 | C February 10, 2003 | |
190 | C | |
191 | C Refs: | |
192 | C | |
193 | C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184. | |
194 | C | |
195 | C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002. | |
196 | C | |
197 | C | |
198 | C This package contains quenching weights for gluon radiation in the | |
199 | C multiple soft scattering approximation. | |
200 | C | |
201 | C swqmult returns the quenching weight for a quark (ipart=1) or | |
202 | C a gluon (ipart=2) traversing a medium with transport coeficient q and | |
203 | C length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where | |
204 | C wc=0.5*q*L^2 and w is the energy radiated. The output values are | |
205 | C the continuous and discrete (prefactor of the delta function) parts | |
206 | C of the quenching weights. | |
207 | C | |
208 | C In order to use this routine, the files cont.all and disc.all need to be | |
209 | C in the working directory. | |
210 | C | |
211 | C An initialization of the tables is needed by doing call initmult before | |
212 | C using swqmult. | |
213 | C | |
214 | C Please, send us any comment: | |
215 | C | |
216 | C urs.wiedemann@cern.ch | |
217 | C carlos.salgado@cern.ch | |
218 | C | |
219 | C | |
220 | C------------------------------------------------------------------- | |
221 | ||
222 | SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete) | |
223 | * | |
224 | REAL*8 xx(400), daq(34), caq(34,261), rrr(34) | |
225 | COMMON /dataqua/ xx, daq, caq, rrr | |
226 | * | |
227 | REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34) | |
228 | COMMON /dataglu/ xxg, dag, cag, rrrg | |
229 | ||
230 | REAL*8 rrrr,xxxx, continuous, discrete | |
231 | REAL*8 rrin, xxin | |
232 | INTEGER nrlow, nrhigh, nxlow, nxhigh | |
233 | REAL*8 rrhigh, rrlow, rfraclow, rfrachigh | |
234 | REAL*8 xfraclow, xfrachigh | |
235 | REAL*8 clow, chigh | |
236 | * | |
237 | ||
238 | continuous=0.d0 | |
239 | discrete=0.d0 | |
240 | ||
241 | rrin = rrrr | |
242 | xxin = xxxx | |
243 | * | |
244 | do 666, nr=1,34 | |
245 | if (rrin.lt.rrr(nr)) then | |
246 | rrhigh = rrr(nr) | |
247 | else | |
248 | rrhigh = rrr(nr-1) | |
249 | rrlow = rrr(nr) | |
250 | nrlow = nr | |
251 | nrhigh = nr-1 | |
252 | goto 665 | |
253 | endif | |
254 | 666 enddo | |
255 | 665 continue | |
256 | * | |
257 | rfraclow = (rrhigh-rrin)/(rrhigh-rrlow) | |
258 | rfrachigh = (rrin-rrlow)/(rrhigh-rrlow) | |
259 | if (rrin.gt.10000d0) then | |
260 | rfraclow = dlog(rrhigh/rrin)/dlog(rrhigh/rrlow) | |
261 | rfrachigh = dlog(rrin/rrlow)/dlog(rrhigh/rrlow) | |
262 | endif | |
263 | * | |
264 | if (ipart.eq.1.and.rrin.ge.rrr(1)) then | |
265 | nrlow=1 | |
266 | nrhigh=1 | |
267 | rfraclow=1 | |
268 | rfrachigh=0 | |
269 | endif | |
270 | ||
271 | if (ipart.ne.1.and.rrin.ge.rrrg(1)) then | |
272 | nrlow=1 | |
273 | nrhigh=1 | |
274 | rfraclow=1 | |
275 | rfrachigh=0 | |
276 | endif | |
277 | ||
278 | if (xxxx.ge.xx(261)) go to 245 | |
279 | ||
280 | nxlow = int(xxin/0.01) + 1 | |
281 | nxhigh = nxlow + 1 | |
282 | xfraclow = (xx(nxhigh)-xxin)/0.01 | |
283 | xfrachigh = (xxin - xx(nxlow))/0.01 | |
284 | * | |
285 | if(ipart.eq.1) then | |
286 | clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh) | |
287 | chigh = xfraclow*caq(nrhigh,nxlow)+xfrachigh*caq(nrhigh,nxhigh) | |
288 | else | |
289 | clow = xfraclow*cag(nrlow,nxlow)+xfrachigh*cag(nrlow,nxhigh) | |
290 | chigh = xfraclow*cag(nrhigh,nxlow)+xfrachigh*cag(nrhigh,nxhigh) | |
291 | endif | |
292 | ||
293 | continuous = rfraclow*clow + rfrachigh*chigh | |
294 | ||
295 | 245 continue | |
296 | ||
297 | if(ipart.eq.1) then | |
298 | discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh) | |
299 | else | |
300 | discrete = rfraclow*dag(nrlow) + rfrachigh*dag(nrhigh) | |
301 | endif | |
302 | * | |
303 | END | |
304 | ||
305 | subroutine initmult | |
306 | REAL*8 xxq(400), daq(34), caq(34,261), rrr(34) | |
307 | COMMON /dataqua/ xxq, daq, caq, rrr | |
308 | * | |
309 | REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34) | |
310 | COMMON /dataglu/ xxg, dag, cag, rrrg | |
311 | * | |
312 | OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90) | |
313 | do 110 nn=1,261 | |
314 | read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn), | |
315 | + caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn), | |
316 | + caq(9,nn), caq(10,nn), caq(11,nn), caq(12,nn), | |
317 | + caq(13,nn), | |
318 | + caq(14,nn), caq(15,nn), caq(16,nn), caq(17,nn), | |
319 | + caq(18,nn), | |
320 | + caq(19,nn), caq(20,nn), caq(21,nn), caq(22,nn), | |
321 | + caq(23,nn), | |
322 | + caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn), | |
323 | + caq(28,nn), | |
324 | + caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn), | |
325 | + caq(33,nn), caq(34,nn) | |
326 | 110 continue | |
327 | do 111 nn=1,261 | |
328 | read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn), | |
329 | + cag(4,nn), cag(5,nn), cag(6,nn), cag(7,nn), cag(8,nn), | |
330 | + cag(9,nn), cag(10,nn), cag(11,nn), cag(12,nn), | |
331 | + cag(13,nn), | |
332 | + cag(14,nn), cag(15,nn), cag(16,nn), cag(17,nn), | |
333 | + cag(18,nn), | |
334 | + cag(19,nn), cag(20,nn), cag(21,nn), cag(22,nn), | |
335 | + cag(23,nn), | |
336 | + cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn), | |
337 | + cag(28,nn), | |
338 | + cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn), | |
339 | + cag(33,nn), cag(34,nn) | |
340 | 111 continue | |
341 | close(20) | |
342 | * | |
343 | OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91) | |
344 | do 112 nn=1,34 | |
345 | read (21,*) rrr(nn), daq(nn) | |
346 | 112 continue | |
347 | do 113 nn=1,34 | |
348 | read (21,*) rrrg(nn), dag(nn) | |
349 | 113 continue | |
350 | close(21) | |
351 | * | |
352 | goto 888 | |
353 | 90 PRINT*, 'input - output error' | |
354 | 91 PRINT*, 'input - output error #2' | |
355 | 888 continue | |
356 | ||
357 | end | |
358 | ||
359 | ||
360 | ======================================================================= | |
361 | ||
362 | Adapted to ROOT macro by A. Dainese - 13/07/2003 | |
363 | Ported to class by C. Loizides - 12/02/2004 | |
364 | New version for extended R values added - 06/03/2004 | |
365 | */ | |
366 | ||
367 | Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx, | |
368 | Double_t &continuous,Double_t &discrete) const | |
369 | { | |
370 | // Calculate Multiple Scattering approx. | |
371 | // weights for given parton type, | |
372 | // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2 | |
373 | ||
374 | //set result to zero | |
375 | continuous=0.; | |
376 | discrete=0.; | |
377 | ||
378 | // read-in data before first call | |
379 | if(!fTablesLoaded){ | |
380 | Error("CalcMult","Tables are not loaded."); | |
381 | return -1; | |
382 | } | |
383 | if(!fMultSoft){ | |
384 | Error("CalcMult","Tables are not loaded for Multiple Scattering."); | |
385 | return -1; | |
386 | } | |
387 | ||
388 | Double_t rrin = rrrr; | |
389 | Double_t xxin = xxxx; | |
390 | ||
391 | if(xxin>fxx[260]) return -1; | |
392 | Int_t nxlow = (Int_t)(xxin/0.01) + 1; | |
393 | Int_t nxhigh = nxlow + 1; | |
394 | Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.01; | |
395 | Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01; | |
396 | ||
397 | //why this? | |
398 | if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD | |
399 | if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD | |
400 | ||
401 | Int_t nrlow=0,nrhigh=0; | |
402 | Double_t rrhigh=0,rrlow=0; | |
403 | for(Int_t nr=1; nr<=34; nr++) { | |
404 | if(rrin<frrr[nr-1]) { | |
405 | rrhigh = frrr[nr-1]; | |
406 | } else { | |
407 | rrhigh = frrr[nr-1-1]; | |
408 | rrlow = frrr[nr-1]; | |
409 | nrlow = nr; | |
410 | nrhigh = nr-1; | |
411 | break; | |
412 | } | |
413 | } | |
414 | ||
415 | rrin = rrrr; // AD | |
416 | ||
417 | Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow); | |
418 | Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow); | |
419 | ||
420 | if(rrin>1.e4){ | |
421 | rfraclow = TMath::Log2(rrhigh/rrin)/TMath::Log2(rrhigh/rrlow); | |
422 | rfrachigh = TMath::Log2(rrin/rrlow)/TMath::Log2(rrhigh/rrlow); | |
423 | } | |
424 | if((ipart==1) && (rrin>=frrr[0])) | |
425 | { | |
426 | nrlow=1; | |
427 | nrhigh=1; | |
428 | rfraclow=1.; | |
429 | rfrachigh=0.; | |
430 | } | |
431 | if((ipart==2) && (rrin>=frrrg[0])) | |
432 | { | |
433 | nrlow=1; | |
434 | nrhigh=1; | |
435 | rfraclow=1.; | |
436 | rfrachigh=0.; | |
437 | } | |
438 | ||
439 | //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD | |
440 | ||
441 | Double_t clow=0,chigh=0; | |
442 | if(ipart==1) { | |
443 | clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1]; | |
444 | chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1]; | |
445 | } else { | |
446 | clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1]; | |
447 | chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1]; | |
448 | } | |
449 | ||
450 | continuous = rfraclow*clow + rfrachigh*chigh; | |
451 | //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n", | |
452 | // rfraclow,clow,rfrachigh,chigh,continuous); | |
453 | ||
454 | if(ipart==1) { | |
455 | discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1]; | |
456 | } else { | |
457 | discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1]; | |
458 | } | |
459 | ||
460 | return 0; | |
461 | } | |
462 | ||
463 | Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *discall) | |
464 | { | |
465 | // read in tables for Single Hard Approx. | |
466 | // path to continuum and to discrete part | |
467 | ||
468 | fTablesLoaded = kFALSE; | |
469 | fMultSoft=kFALSE; | |
470 | ||
471 | Char_t fname[1024]; | |
472 | sprintf(fname,"%s",gSystem->ExpandPathName(contall)); | |
473 | //PH ifstream fincont(fname); | |
474 | fstream fincont(fname,ios::in); | |
475 | #if defined(__HP_aCC) || defined(__DECCXX) | |
476 | if(!fincont.rdbuf()->is_open()) return -1; | |
477 | #else | |
478 | if(!fincont.is_open()) return -1; | |
479 | #endif | |
480 | ||
481 | Int_t nn=0; //quarks | |
482 | while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>> | |
483 | fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>> | |
484 | fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>> | |
485 | fcaq[13][nn]>> | |
486 | fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>> | |
487 | fcaq[18][nn]>> | |
488 | fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>> | |
489 | fcaq[23][nn]>> | |
490 | fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>> | |
491 | fcaq[28][nn]>> | |
492 | fcaq[29][nn]) | |
493 | { | |
494 | nn++; | |
495 | if(nn==261) break; | |
496 | } | |
497 | ||
498 | nn=0; //gluons | |
499 | while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>> | |
500 | fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>> | |
501 | fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>> | |
502 | fcag[13][nn]>> | |
503 | fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>> | |
504 | fcag[18][nn]>> | |
505 | fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>> | |
506 | fcag[23][nn]>> | |
507 | fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>> | |
508 | fcag[28][nn]>> | |
509 | fcag[29][nn]) { | |
510 | nn++; | |
511 | if(nn==261) break; | |
512 | } | |
513 | fincont.close(); | |
514 | ||
515 | sprintf(fname,"%s",gSystem->ExpandPathName(discall)); | |
516 | //PH ifstream findisc(fname); | |
517 | fstream findisc(fname,ios::in); | |
518 | #if defined(__HP_aCC) || defined(__DECCXX) | |
519 | if(!findisc.rdbuf()->is_open()) return -1; | |
520 | #else | |
521 | if(!findisc.is_open()) return -1; | |
522 | #endif | |
523 | ||
524 | nn=0; //quarks | |
525 | while(findisc>>frrr[nn]>>fdaq[nn]) { | |
526 | nn++; | |
527 | if(nn==30) break; | |
528 | } | |
529 | nn=0; //gluons | |
530 | while(findisc>>frrrg[nn]>>fdag[nn]) { | |
531 | nn++; | |
532 | if(nn==30) break; | |
533 | } | |
534 | findisc.close(); | |
535 | ||
536 | fTablesLoaded = kTRUE; | |
537 | return 0; | |
538 | } | |
539 | ||
540 | /* | |
541 | C*************************************************************************** | |
542 | C Quenching Weights for Single Hard Scattering | |
543 | C February 20, 2003 | |
544 | C | |
545 | C Refs: | |
546 | C | |
547 | C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184. | |
548 | C | |
549 | C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002. | |
550 | C | |
551 | C | |
552 | C This package contains quenching weights for gluon radiation in the | |
553 | C single hard scattering approximation. | |
554 | C | |
555 | C swqlin returns the quenching weight for a quark (ipart=1) or | |
556 | C a gluon (ipart=2) traversing a medium with Debye screening mass mu and | |
557 | C length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where | |
558 | C wc=0.5*mu^2*L and w is the energy radiated. The output values are | |
559 | C the continuous and discrete (prefactor of the delta function) parts | |
560 | C of the quenching weights. | |
561 | C | |
562 | C In order to use this routine, the files contlin.all and disclin.all | |
563 | C need to be in the working directory. | |
564 | C | |
565 | C An initialization of the tables is needed by doing call initlin before | |
566 | C using swqlin. | |
567 | C | |
568 | C Please, send us any comment: | |
569 | C | |
570 | C urs.wiedemann@cern.ch | |
571 | C carlos.salgado@cern.ch | |
572 | C | |
573 | C | |
574 | C------------------------------------------------------------------- | |
575 | ||
576 | ||
577 | SUBROUTINE swqlin(ipart,rrrr,xxxx,continuous,discrete) | |
578 | * | |
579 | REAL*8 xx(400), dalq(30), calq(30,261), rrr(30) | |
580 | COMMON /datalinqua/ xx, dalq, calq, rrr | |
581 | * | |
582 | REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30) | |
583 | COMMON /datalinglu/ xxlg, dalg, calg, rrrlg | |
584 | ||
585 | REAL*8 rrrr,xxxx, continuous, discrete | |
586 | REAL*8 rrin, xxin | |
587 | INTEGER nrlow, nrhigh, nxlow, nxhigh | |
588 | REAL*8 rrhigh, rrlow, rfraclow, rfrachigh | |
589 | REAL*8 xfraclow, xfrachigh | |
590 | REAL*8 clow, chigh | |
591 | * | |
592 | rrin = rrrr | |
593 | xxin = xxxx | |
594 | * | |
595 | nxlow = int(xxin/0.038) + 1 | |
596 | nxhigh = nxlow + 1 | |
597 | xfraclow = (xx(nxhigh)-xxin)/0.038 | |
598 | xfrachigh = (xxin - xx(nxlow))/0.038 | |
599 | * | |
600 | do 666, nr=1,30 | |
601 | if (rrin.lt.rrr(nr)) then | |
602 | rrhigh = rrr(nr) | |
603 | else | |
604 | rrhigh = rrr(nr-1) | |
605 | rrlow = rrr(nr) | |
606 | nrlow = nr | |
607 | nrhigh = nr-1 | |
608 | goto 665 | |
609 | endif | |
610 | 666 enddo | |
611 | 665 continue | |
612 | * | |
613 | rfraclow = (rrhigh-rrin)/(rrhigh-rrlow) | |
614 | rfrachigh = (rrin-rrlow)/(rrhigh-rrlow) | |
615 | * | |
616 | if(ipart.eq.1) then | |
617 | clow = xfraclow*calq(nrlow,nxlow)+xfrachigh*calq(nrlow,nxhigh) | |
618 | chigh = xfraclow*calq(nrhigh,nxlow)+xfrachigh*calq(nrhigh,nxhigh) | |
619 | else | |
620 | clow = xfraclow*calg(nrlow,nxlow)+xfrachigh*calg(nrlow,nxhigh) | |
621 | chigh = xfraclow*calg(nrhigh,nxlow)+xfrachigh*calg(nrhigh,nxhigh) | |
622 | endif | |
623 | ||
624 | continuous = rfraclow*clow + rfrachigh*chigh | |
625 | ||
626 | if(ipart.eq.1) then | |
627 | discrete = rfraclow*dalq(nrlow) + rfrachigh*dalq(nrhigh) | |
628 | else | |
629 | discrete = rfraclow*dalg(nrlow) + rfrachigh*dalg(nrhigh) | |
630 | endif | |
631 | * | |
632 | END | |
633 | ||
634 | subroutine initlin | |
635 | REAL*8 xxlq(400), dalq(30), calq(30,261), rrr(30) | |
636 | COMMON /datalinqua/ xxlq, dalq, calq, rrr | |
637 | * | |
638 | REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30) | |
639 | COMMON /datalinglu/ xxlg, dalg, calg, rrrlg | |
640 | * | |
641 | OPEN(UNIT=20,FILE='contlin.all',STATUS='OLD',ERR=90) | |
642 | do 110 nn=1,261 | |
643 | read (20,*) xxlq(nn), calq(1,nn), calq(2,nn), calq(3,nn), | |
644 | + calq(4,nn), calq(5,nn), calq(6,nn), calq(7,nn), calq(8,nn), | |
645 | + calq(9,nn), calq(10,nn), calq(11,nn), calq(12,nn), | |
646 | + calq(13,nn), | |
647 | + calq(14,nn), calq(15,nn), calq(16,nn), calq(17,nn), | |
648 | + calq(18,nn), | |
649 | + calq(19,nn), calq(20,nn), calq(21,nn), calq(22,nn), | |
650 | + calq(23,nn), | |
651 | + calq(24,nn), calq(25,nn), calq(26,nn), calq(27,nn), | |
652 | + calq(28,nn), | |
653 | + calq(29,nn), calq(30,nn) | |
654 | 110 continue | |
655 | do 111 nn=1,261 | |
656 | read (20,*) xxlg(nn), calg(1,nn), calg(2,nn), calg(3,nn), | |
657 | + calg(4,nn), calg(5,nn), calg(6,nn), calg(7,nn), calg(8,nn), | |
658 | + calg(9,nn), calg(10,nn), calg(11,nn), calg(12,nn), | |
659 | + calg(13,nn), | |
660 | + calg(14,nn), calg(15,nn), calg(16,nn), calg(17,nn), | |
661 | + calg(18,nn), | |
662 | + calg(19,nn), calg(20,nn), calg(21,nn), calg(22,nn), | |
663 | + calg(23,nn), | |
664 | + calg(24,nn), calg(25,nn), calg(26,nn), calg(27,nn), | |
665 | + calg(28,nn), | |
666 | + calg(29,nn), calg(30,nn) | |
667 | 111 continue | |
668 | close(20) | |
669 | * | |
670 | OPEN(UNIT=21,FILE='disclin.all',STATUS='OLD',ERR=91) | |
671 | do 112 nn=1,30 | |
672 | read (21,*) rrr(nn), dalq(nn) | |
673 | 112 continue | |
674 | do 113 nn=1,30 | |
675 | read (21,*) rrrlg(nn), dalg(nn) | |
676 | 113 continue | |
677 | close(21) | |
678 | * | |
679 | goto 888 | |
680 | 90 PRINT*, 'input - output error' | |
681 | 91 PRINT*, 'input - output error #2' | |
682 | 888 continue | |
683 | ||
684 | end | |
685 | ||
686 | ======================================================================= | |
687 | ||
688 | Ported to class by C. Loizides - 17/02/2004 | |
689 | ||
690 | */ | |
691 | ||
692 | Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx, | |
693 | Double_t &continuous,Double_t &discrete) const | |
694 | { | |
695 | // calculate Single Hard approx. | |
696 | // weights for given parton type, | |
697 | // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L | |
698 | ||
699 | // read-in data before first call | |
700 | if(!fTablesLoaded){ | |
701 | Error("CalcMult","Tables are not loaded."); | |
702 | return -1; | |
703 | } | |
704 | if(!fMultSoft){ | |
705 | Error("CalcMult","Tables are not loaded for Single Hard Scattering."); | |
706 | return -1; | |
707 | } | |
708 | ||
709 | Double_t rrin = rrrr; | |
710 | Double_t xxin = xxxx; | |
711 | ||
712 | Int_t nxlow = (Int_t)(xxin/0.038) + 1; | |
713 | Int_t nxhigh = nxlow + 1; | |
714 | Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.038; | |
715 | Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.038; | |
716 | ||
717 | //why this? | |
718 | if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD | |
719 | if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD | |
720 | ||
721 | Int_t nrlow=0,nrhigh=0; | |
722 | Double_t rrhigh=0,rrlow=0; | |
723 | for(Int_t nr=1; nr<=30; nr++) { | |
724 | if(rrin<frrr[nr-1]) { | |
725 | rrhigh = frrr[nr-1]; | |
726 | } else { | |
727 | rrhigh = frrr[nr-1-1]; | |
728 | rrlow = frrr[nr-1]; | |
729 | nrlow = nr; | |
730 | nrhigh = nr-1; | |
731 | break; | |
732 | } | |
733 | } | |
734 | ||
735 | rrin = rrrr; // AD | |
736 | ||
737 | Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow); | |
738 | Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow); | |
739 | ||
740 | //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD | |
741 | ||
742 | Double_t clow=0,chigh=0; | |
743 | if(ipart==1) { | |
744 | clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1]; | |
745 | chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1]; | |
746 | } else { | |
747 | clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1]; | |
748 | chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1]; | |
749 | } | |
750 | ||
751 | continuous = rfraclow*clow + rfrachigh*chigh; | |
752 | //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n", | |
753 | // rfraclow,clow,rfrachigh,chigh,continuous); | |
754 | ||
755 | if(ipart==1) { | |
756 | discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1]; | |
757 | } else { | |
758 | discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1]; | |
759 | } | |
760 | ||
761 | return 0; | |
762 | } | |
763 | ||
764 | Int_t AliQuenchingWeights::CalcMult(Int_t ipart, | |
765 | Double_t w,Double_t qtransp,Double_t length, | |
766 | Double_t &continuous,Double_t &discrete) const | |
767 | { | |
768 | // Calculate Multiple Scattering approx. | |
769 | // weights for given parton type, | |
770 | // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2 | |
771 | ||
772 | Double_t wc=CalcWC(qtransp,length); | |
773 | Double_t rrrr=CalcR(wc,length); | |
774 | Double_t xxxx=w/wc; | |
775 | return CalcMult(ipart,rrrr,xxxx,continuous,discrete); | |
776 | } | |
777 | ||
778 | Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, | |
779 | Double_t w,Double_t mu,Double_t length, | |
780 | Double_t &continuous,Double_t &discrete) const | |
781 | { | |
782 | // calculate Single Hard approx. | |
783 | // weights for given parton type, | |
784 | // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L | |
785 | ||
786 | Double_t wcbar=CalcWCbar(mu,length); | |
787 | Double_t rrrr=CalcR(wcbar,length); | |
788 | Double_t xxxx=w/wcbar; | |
789 | return CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete); | |
790 | } | |
791 | ||
792 | Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const | |
793 | { | |
794 | //calculate R value and | |
795 | //check if it is less then maximum | |
796 | ||
797 | Double_t R = wc*l*fgkConvFmToInvGeV; | |
798 | if(R>fgkRMax) { | |
799 | Warning("CalcR","Value of R = %.2f; should be less than %.2f",R,fgkRMax); | |
800 | return -R; | |
801 | } | |
802 | return R; | |
803 | } | |
804 | ||
805 | Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const | |
806 | { | |
807 | // return DeltaE for MS or SH scattering | |
808 | // for given parton type, length and energy | |
809 | // Dependant on ECM (energy constraint method) | |
810 | // e is used to determine where to set bins to zero. | |
811 | ||
812 | if(!fHistos){ | |
813 | Fatal("GetELossRandom","Call SampleEnergyLoss method before!"); | |
814 | return -1000.; | |
815 | } | |
816 | if((ipart<1) || (ipart>2)) { | |
817 | Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart); | |
818 | return -1000; | |
819 | } | |
820 | ||
821 | Int_t l=(Int_t)length; | |
822 | if((length-(Double_t)l)>0.5) l++; | |
823 | if(l<=0) return 0.; | |
824 | if(l>fLengthMax) l=fLengthMax; | |
825 | ||
826 | if(fECMethod==kReweight){ | |
827 | TH1F *dummy=new TH1F(*fHistos[ipart-1][l-1]); | |
828 | dummy->SetName("dummy"); | |
829 | for(Int_t bin=dummy->FindBin(e)+1;bin<=1100;bin++){ | |
830 | dummy->SetBinContent(bin,0.); | |
831 | } | |
832 | dummy->Scale(1./dummy->Integral()); | |
833 | Double_t ret=dummy->GetRandom(); | |
834 | delete dummy; | |
835 | return ret; | |
836 | //****** !! ALTERNATIVE WAY OF DOING IT !! *** | |
837 | //Double_t ret = 2.*e; | |
838 | //while(ret>e) ret=fHistos[ipart-1][l-1]->GetRandom(); | |
839 | //return ret; | |
840 | //******************************************** | |
841 | } else { //kDefault | |
842 | Double_t ret=fHistos[ipart-1][l-1]->GetRandom(); | |
843 | if(ret>e) return e; | |
844 | return ret; | |
845 | } | |
846 | } | |
847 | ||
848 | Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const | |
849 | { | |
850 | //return quenched parton energy | |
851 | //for given parton type, length and energy | |
852 | ||
853 | Double_t loss=GetELossRandom(ipart,length,e); | |
854 | return e-loss; | |
855 | } | |
856 | ||
857 | Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const | |
858 | { | |
859 | // return DeltaE for MS or SH scattering | |
860 | // for given parton type, length distribution and energy | |
861 | // Dependant on ECM (energy constraint method) | |
862 | // e is used to determine where to set bins to zero. | |
863 | ||
864 | if(!hell){ | |
865 | Warning("GetELossRandom","Pointer to length distribution is NULL."); | |
866 | return 0.; | |
867 | } | |
868 | Double_t ell=hell->GetRandom(); | |
869 | return GetELossRandom(ipart,ell,e); | |
870 | } | |
871 | ||
872 | Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e) const | |
873 | { | |
874 | //return quenched parton energy | |
875 | //for given parton type, length distribution and energy | |
876 | ||
877 | Double_t loss=GetELossRandom(ipart,hell,e); | |
878 | return e-loss; | |
879 | } | |
880 | ||
881 | Int_t AliQuenchingWeights::SampleEnergyLoss() | |
882 | { | |
883 | // Has to be called to fill the histograms | |
884 | // | |
885 | // For stored values fQTransport loop over | |
886 | // particle type and length = 1 to fMaxLength (fm) | |
887 | // to fill energy loss histos | |
888 | // | |
889 | // Take histogram of continuous weights | |
890 | // Take discrete_weight | |
891 | // If discrete_weight > 1, put all channels to 0, except channel 1 | |
892 | // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral | |
893 | ||
894 | // read-in data before first call | |
895 | if(!fTablesLoaded){ | |
896 | Error("CalcMult","Tables are not loaded."); | |
897 | return -1; | |
898 | } | |
899 | ||
900 | if(fMultSoft) { | |
901 | Int_t lmax=CalcLengthMax(fQTransport); | |
902 | if(fLengthMax>lmax){ | |
903 | Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax); | |
904 | fLengthMax=lmax; | |
905 | } | |
906 | } else { | |
907 | Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested."); | |
908 | } | |
909 | ||
910 | Reset(); | |
911 | fHistos=new TH1F**[2]; | |
912 | fHistos[0]=new TH1F*[fLengthMax]; | |
913 | fHistos[1]=new TH1F*[fLengthMax]; | |
914 | fLengthMaxOld=fLengthMax; //remember old value in case | |
915 | //user wants to reset | |
916 | ||
917 | Int_t medvalue=0; | |
918 | Char_t meddesc[100]; | |
919 | if(fMultSoft) { | |
920 | medvalue=(Int_t)(fQTransport*1000.); | |
921 | sprintf(meddesc,"MS"); | |
922 | } else { | |
923 | medvalue=(Int_t)(fMu*1000.); | |
924 | sprintf(meddesc,"SH"); | |
925 | } | |
926 | ||
927 | for(Int_t ipart=1;ipart<=2;ipart++){ | |
928 | for(Int_t l=1;l<=fLengthMax;l++){ | |
929 | ||
930 | Char_t hname[100]; | |
931 | sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l); | |
932 | Double_t wc = CalcWC(l); | |
933 | fHistos[ipart-1][l-1] = new TH1F(hname,hname,1100,0.,1.1*wc); | |
934 | fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]"); | |
935 | fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)"); | |
936 | fHistos[ipart-1][l-1]->SetLineColor(4); | |
937 | ||
938 | Double_t rrrr = CalcR(wc,l); | |
939 | Double_t discrete=0.; | |
940 | // loop on histogram channels | |
941 | for(Int_t bin=1; bin<=1100; bin++) { | |
942 | Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc; | |
943 | Double_t continuous; | |
944 | CalcMult(ipart,rrrr,xxxx,continuous,discrete); | |
945 | fHistos[ipart-1][l-1]->SetBinContent(bin,continuous); | |
946 | } | |
947 | // add discrete part to distribution | |
948 | if(discrete>=1.) | |
949 | for(Int_t bin=2;bin<=1100;bin++) | |
950 | fHistos[ipart-1][l-1]->SetBinContent(bin,0.); | |
951 | else { | |
952 | Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,1100); | |
953 | fHistos[ipart-1][l-1]->Fill(0.,val); | |
954 | } | |
955 | Double_t hint=fHistos[ipart-1][l-1]->Integral(1,1100); | |
956 | fHistos[ipart-1][l-1]->Scale(1./hint); | |
957 | } | |
958 | } | |
959 | return 0; | |
960 | } | |
961 | ||
962 | const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Int_t l) const | |
963 | { | |
964 | //return quenching histograms | |
965 | //for ipart and length | |
966 | ||
967 | if(!fHistos){ | |
968 | Fatal("GetELossRandom","Call SampleEnergyLoss method before!"); | |
969 | return 0; | |
970 | } | |
971 | if((ipart<1) || (ipart>2)) { | |
972 | Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart); | |
973 | return 0; | |
974 | } | |
975 | ||
976 | if(l<=0) return 0; | |
977 | if(l>fLengthMax) l=fLengthMax; | |
978 | ||
979 | return fHistos[ipart-1][l-1]; | |
980 | } | |
981 | ||
982 | TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const | |
983 | { | |
984 | // ipart = 1 for quark, 2 for gluon | |
985 | // medval a) qtransp = transport coefficient (GeV^2/fm) | |
986 | // b) mu = Debye mass (GeV) | |
987 | // length = path length in medium (fm) | |
988 | // Get from SW tables: | |
989 | // - discrete weight (probability to have NO energy loss) | |
990 | // - continuous weight, as a function of dE | |
991 | // compute up to max dE = 1.1*wc | |
992 | // (as an histogram named hContQW_<ipart>_<medval*1000>_<l> | |
993 | // and range [0,1.1*wc] (1100 channels) | |
994 | ||
995 | Double_t wc = 0; | |
996 | Char_t meddesc[100]; | |
997 | if(fMultSoft) { | |
998 | wc=CalcWC(medval,length); | |
999 | sprintf(meddesc,"MS"); | |
1000 | } else { | |
1001 | wc=CalcWCbar(medval,length); | |
1002 | sprintf(meddesc,"SH"); | |
1003 | } | |
1004 | ||
1005 | Char_t hname[100]; | |
1006 | sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart, | |
1007 | (Int_t)(medval*1000.),(Int_t)length); | |
1008 | ||
1009 | TH1F *hist = new TH1F("hist",hname,1100,0.,1.1*wc); | |
1010 | hist->SetXTitle("#Delta E [GeV]"); | |
1011 | hist->SetYTitle("p(#Delta E)"); | |
1012 | hist->SetLineColor(4); | |
1013 | ||
1014 | Double_t rrrr = CalcR(wc,length); | |
1015 | // loop on histogram channels | |
1016 | for(Int_t bin=1; bin<=1100; bin++) { | |
1017 | Double_t xxxx = hist->GetBinCenter(bin)/wc; | |
1018 | Double_t continuous,discrete; | |
1019 | Int_t ret=0; | |
1020 | if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete); | |
1021 | else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete); | |
1022 | if(ret!=0){ | |
1023 | delete hist; | |
1024 | return 0; | |
1025 | }; | |
1026 | hist->SetBinContent(bin,continuous); | |
1027 | } | |
1028 | return hist; | |
1029 | } | |
1030 | ||
1031 | TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const | |
1032 | { | |
1033 | // ipart = 1 for quark, 2 for gluon | |
1034 | // medval a) qtransp = transport coefficient (GeV^2/fm) | |
1035 | // b) mu = Debye mass (GeV) | |
1036 | // length = path length in medium (fm) | |
1037 | // Get from SW tables: | |
1038 | // - discrete weight (probability to have NO energy loss) | |
1039 | // - continuous weight, as a function of dE | |
1040 | // compute up to max dE = 1.1*wc | |
1041 | // (as an histogram named hContQW_<ipart>_<medval*1000>_<l> | |
1042 | // and range [0,1.1] (1100 channels) | |
1043 | ||
1044 | Double_t wc = 0; | |
1045 | Char_t meddesc[100]; | |
1046 | if(fMultSoft) { | |
1047 | wc=CalcWC(medval,length); | |
1048 | sprintf(meddesc,"MS"); | |
1049 | } else { | |
1050 | wc=CalcWCbar(medval,length); | |
1051 | sprintf(meddesc,"SH"); | |
1052 | } | |
1053 | ||
1054 | Char_t hname[100]; | |
1055 | sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart, | |
1056 | (Int_t)(medval*1000.),(Int_t)length); | |
1057 | ||
1058 | TH1F *histx = new TH1F("histx",hname,1100,0.,1.1); | |
1059 | histx->SetXTitle("x = #Delta E/#omega_{c}"); | |
1060 | if(fMultSoft) | |
1061 | histx->SetYTitle("p(#Delta E/#omega_{c})"); | |
1062 | else | |
1063 | histx->SetYTitle("p(#Delta E/#bar#omega_{c})"); | |
1064 | histx->SetLineColor(4); | |
1065 | ||
1066 | Double_t rrrr = CalcR(wc,length); | |
1067 | // loop on histogram channels | |
1068 | for(Int_t bin=1; bin<=1100; bin++) { | |
1069 | Double_t xxxx = histx->GetBinCenter(bin); | |
1070 | Double_t continuous,discrete; | |
1071 | Int_t ret=0; | |
1072 | if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete); | |
1073 | else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete); | |
1074 | if(ret!=0){ | |
1075 | delete histx; | |
1076 | return 0; | |
1077 | }; | |
1078 | histx->SetBinContent(bin,continuous); | |
1079 | } | |
1080 | return histx; | |
1081 | } | |
1082 | ||
1083 | TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const | |
1084 | { | |
1085 | //compute energy loss histogram for | |
1086 | //parton type, medium value, length and energy | |
1087 | ||
1088 | AliQuenchingWeights *dummy=new AliQuenchingWeights(*this); | |
1089 | if(fMultSoft){ | |
1090 | dummy->SetQTransport(medval); | |
1091 | dummy->InitMult(); | |
1092 | } else { | |
1093 | dummy->SetMu(medval); | |
1094 | dummy->InitSingleHard(); | |
1095 | } | |
1096 | dummy->SampleEnergyLoss(); | |
1097 | ||
1098 | Char_t name[100]; | |
1099 | Char_t hname[100]; | |
1100 | if(ipart==1){ | |
1101 | sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); | |
1102 | sprintf(hname,"hLossQuarks"); | |
1103 | } else { | |
1104 | sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); | |
1105 | sprintf(hname,"hLossGluons"); | |
1106 | } | |
1107 | ||
1108 | TH1F *h = new TH1F(hname,name,250,0,250); | |
1109 | for(Int_t i=0;i<100000;i++){ | |
1110 | //if(i % 1000 == 0) cout << "." << flush; | |
1111 | Double_t loss=dummy->GetELossRandom(ipart,l,e); | |
1112 | h->Fill(loss); | |
1113 | } | |
1114 | ||
1115 | h->SetStats(kTRUE); | |
1116 | ||
1117 | delete dummy; | |
1118 | return h; | |
1119 | } | |
1120 | ||
1121 | TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const | |
1122 | { | |
1123 | //compute energy loss histogram for | |
1124 | //parton type, medium value, | |
1125 | //length distribution and energy | |
1126 | ||
1127 | AliQuenchingWeights *dummy=new AliQuenchingWeights(*this); | |
1128 | if(fMultSoft){ | |
1129 | dummy->SetQTransport(medval); | |
1130 | dummy->InitMult(); | |
1131 | } else { | |
1132 | dummy->SetMu(medval); | |
1133 | dummy->InitSingleHard(); | |
1134 | } | |
1135 | dummy->SampleEnergyLoss(); | |
1136 | ||
1137 | Char_t name[100]; | |
1138 | Char_t hname[100]; | |
1139 | if(ipart==1){ | |
1140 | sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); | |
1141 | sprintf(hname,"hLossQuarks"); | |
1142 | } else { | |
1143 | sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); | |
1144 | sprintf(hname,"hLossGluons"); | |
1145 | } | |
1146 | ||
1147 | TH1F *h = new TH1F(hname,name,250,0,250); | |
1148 | for(Int_t i=0;i<100000;i++){ | |
1149 | //if(i % 1000 == 0) cout << "." << flush; | |
1150 | Double_t loss=dummy->GetELossRandom(ipart,hEll,e); | |
1151 | h->Fill(loss); | |
1152 | } | |
1153 | ||
1154 | h->SetStats(kTRUE); | |
1155 | ||
1156 | delete dummy; | |
1157 | return h; | |
1158 | } | |
1159 | ||
1160 | void AliQuenchingWeights::PlotDiscreteWeights(Int_t len) const | |
1161 | { | |
1162 | //plot discrete weights for given length | |
1163 | ||
1164 | TCanvas *c; | |
1165 | if(fMultSoft) | |
1166 | c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400); | |
1167 | else | |
1168 | c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400); | |
1169 | c->cd(); | |
1170 | ||
1171 | TH2F *hframe = new TH2F("hdisc","",2,0,5.1,2,0,1); | |
1172 | hframe->SetStats(0); | |
1173 | if(fMultSoft) | |
1174 | hframe->SetXTitle("#hat{q} [GeV^{2}/fm]"); | |
1175 | else | |
1176 | hframe->SetXTitle("#mu [GeV]"); | |
1177 | hframe->SetYTitle("Probability #Delta E = 0 , p_{0}"); | |
1178 | hframe->Draw(); | |
1179 | ||
1180 | TGraph *gq=new TGraph(20); | |
1181 | Int_t i=0; | |
1182 | if(fMultSoft) { | |
1183 | for(Double_t q=0.05;q<=5.05;q+=0.25){ | |
1184 | Double_t disc,cont; | |
1185 | CalcMult(1,1.0, q, len, cont, disc); | |
1186 | gq->SetPoint(i,q,disc);i++; | |
1187 | } | |
1188 | } else { | |
1189 | for(Double_t m=0.05;m<=5.05;m+=0.25){ | |
1190 | Double_t disc,cont; | |
1191 | CalcSingleHard(1,1.0, m, len, cont, disc); | |
1192 | gq->SetPoint(i,m,disc);i++; | |
1193 | } | |
1194 | } | |
1195 | gq->SetMarkerStyle(20); | |
1196 | gq->Draw("pl"); | |
1197 | ||
1198 | TGraph *gg=new TGraph(20); | |
1199 | i=0; | |
1200 | if(fMultSoft){ | |
1201 | for(Double_t q=0.05;q<=5.05;q+=0.25){ | |
1202 | Double_t disc,cont; | |
1203 | CalcMult(2,1.0, q, 5., cont, disc); | |
1204 | gg->SetPoint(i,q,disc);i++; | |
1205 | } | |
1206 | } else { | |
1207 | for(Double_t m=0.05;m<=5.05;m+=0.25){ | |
1208 | Double_t disc,cont; | |
1209 | CalcSingleHard(2,1.0, m, 5., cont, disc); | |
1210 | gg->SetPoint(i,m,disc);i++; | |
1211 | } | |
1212 | } | |
1213 | gg->SetMarkerStyle(24); | |
1214 | gg->Draw("pl"); | |
1215 | ||
1216 | TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); | |
1217 | l1a->SetFillStyle(0); | |
1218 | l1a->SetBorderSize(0); | |
1219 | Char_t label[100]; | |
1220 | sprintf(label,"L = %d fm",len); | |
1221 | l1a->AddEntry(gq,label,""); | |
1222 | l1a->AddEntry(gq,"quark","pl"); | |
1223 | l1a->AddEntry(gg,"gluon","pl"); | |
1224 | l1a->Draw(); | |
1225 | ||
1226 | c->Update(); | |
1227 | } | |
1228 | ||
1229 | void AliQuenchingWeights::PlotContWeights(Int_t itype,Int_t ell) const | |
1230 | { | |
1231 | //plot continous weights for | |
1232 | //given parton type and length | |
1233 | ||
1234 | Float_t medvals[3]; | |
1235 | Char_t title[1024]; | |
1236 | Char_t name[1024]; | |
1237 | if(fMultSoft) { | |
1238 | if(itype==1) | |
1239 | sprintf(title,"Cont. Weight for Multiple Scattering - Quarks"); | |
1240 | else if(itype==2) | |
1241 | sprintf(title,"Cont. Weight for Multiple Scattering - Gluons"); | |
1242 | else return; | |
1243 | medvals[0]=4;medvals[1]=1;medvals[2]=0.5; | |
1244 | sprintf(name,"ccont-ms-%d",itype); | |
1245 | } else { | |
1246 | if(itype==1) | |
1247 | sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks"); | |
1248 | else if(itype==2) | |
1249 | sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons"); | |
1250 | else return; | |
1251 | medvals[0]=2;medvals[1]=1;medvals[2]=0.5; | |
1252 | sprintf(name,"ccont-ms-%d",itype); | |
1253 | } | |
1254 | ||
1255 | TCanvas *c = new TCanvas(name,title,0,0,500,400); | |
1256 | c->cd(); | |
1257 | TH1F *h1=ComputeQWHisto(itype,medvals[0],ell); | |
1258 | h1->SetName("h1"); | |
1259 | h1->SetTitle(title); | |
1260 | h1->SetStats(0); | |
1261 | h1->SetLineColor(1); | |
1262 | h1->DrawCopy(); | |
1263 | TH1F *h2=ComputeQWHisto(itype,medvals[1],ell); | |
1264 | h2->SetName("h2"); | |
1265 | h2->SetLineColor(2); | |
1266 | h2->DrawCopy("SAME"); | |
1267 | TH1F *h3=ComputeQWHisto(itype,medvals[2],ell); | |
1268 | h3->SetName("h3"); | |
1269 | h3->SetLineColor(3); | |
1270 | h3->DrawCopy("SAME"); | |
1271 | ||
1272 | TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); | |
1273 | l1a->SetFillStyle(0); | |
1274 | l1a->SetBorderSize(0); | |
1275 | Char_t label[100]; | |
1276 | sprintf(label,"L = %d fm",ell); | |
1277 | l1a->AddEntry(h1,label,""); | |
1278 | if(fMultSoft) { | |
1279 | sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]); | |
1280 | l1a->AddEntry(h1,label,"pl"); | |
1281 | sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]); | |
1282 | l1a->AddEntry(h2,label,"pl"); | |
1283 | sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]); | |
1284 | l1a->AddEntry(h3,label,"pl"); | |
1285 | } else { | |
1286 | sprintf(label,"#mu = %.1f GeV",medvals[0]); | |
1287 | l1a->AddEntry(h1,label,"pl"); | |
1288 | sprintf(label,"#mu = %.1f GeV",medvals[1]); | |
1289 | l1a->AddEntry(h2,label,"pl"); | |
1290 | sprintf(label,"#mu = %.1f GeV",medvals[2]); | |
1291 | l1a->AddEntry(h3,label,"pl"); | |
1292 | } | |
1293 | l1a->Draw(); | |
1294 | ||
1295 | c->Update(); | |
1296 | } | |
1297 | ||
1298 | void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const | |
1299 | { | |
1300 | //plot continous weights for | |
1301 | //given parton type and medium value | |
1302 | ||
1303 | Char_t title[1024]; | |
1304 | Char_t name[1024]; | |
1305 | if(fMultSoft) { | |
1306 | if(itype==1) | |
1307 | sprintf(title,"Cont. Weight for Multiple Scattering - Quarks"); | |
1308 | else if(itype==2) | |
1309 | sprintf(title,"Cont. Weight for Multiple Scattering - Gluons"); | |
1310 | else return; | |
1311 | sprintf(name,"ccont2-ms-%d",itype); | |
1312 | } else { | |
1313 | if(itype==1) | |
1314 | sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks"); | |
1315 | else if(itype==2) | |
1316 | sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons"); | |
1317 | else return; | |
1318 | sprintf(name,"ccont2-sh-%d",itype); | |
1319 | } | |
1320 | TCanvas *c = new TCanvas(name,title,0,0,500,400); | |
1321 | c->cd(); | |
1322 | TH1F *h1=ComputeQWHisto(itype,medval,8); | |
1323 | h1->SetName("h1"); | |
1324 | h1->SetTitle(title); | |
1325 | h1->SetStats(0); | |
1326 | h1->SetLineColor(1); | |
1327 | h1->DrawCopy(); | |
1328 | TH1F *h2=ComputeQWHisto(itype,medval,5); | |
1329 | h2->SetName("h2"); | |
1330 | h2->SetLineColor(2); | |
1331 | h2->DrawCopy("SAME"); | |
1332 | TH1F *h3=ComputeQWHisto(itype,medval,2); | |
1333 | h3->SetName("h3"); | |
1334 | h3->SetLineColor(3); | |
1335 | h3->DrawCopy("SAME"); | |
1336 | ||
1337 | TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); | |
1338 | l1a->SetFillStyle(0); | |
1339 | l1a->SetBorderSize(0); | |
1340 | Char_t label[100]; | |
1341 | if(fMultSoft) | |
1342 | sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval); | |
1343 | else | |
1344 | sprintf(label,"#mu = %.1f GeV",medval); | |
1345 | ||
1346 | l1a->AddEntry(h1,label,""); | |
1347 | l1a->AddEntry(h1,"L = 8 fm","pl"); | |
1348 | l1a->AddEntry(h2,"L = 5 fm","pl"); | |
1349 | l1a->AddEntry(h3,"L = 2 fm","pl"); | |
1350 | l1a->Draw(); | |
1351 | ||
1352 | c->Update(); | |
1353 | } | |
1354 | ||
1355 | void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e) const | |
1356 | { | |
1357 | //plot average energy loss for given length | |
1358 | //and parton energy | |
1359 | ||
1360 | if(!fTablesLoaded){ | |
1361 | Error("CalcMult","Tables are not loaded."); | |
1362 | return; | |
1363 | } | |
1364 | ||
1365 | Char_t title[1024]; | |
1366 | Char_t name[1024]; | |
1367 | if(fMultSoft){ | |
1368 | sprintf(title,"Average Loss for Multiple Scattering"); | |
1369 | sprintf(name,"cavgelossms"); | |
1370 | } else { | |
1371 | sprintf(title,"Average Loss for Single Hard Scattering"); | |
1372 | sprintf(name,"cavgelosssh"); | |
1373 | } | |
1374 | ||
1375 | TCanvas *c = new TCanvas(name,title,0,0,500,400); | |
1376 | c->cd(); | |
1377 | TH2F *hframe = new TH2F("avgloss",title,2,0,5.1,2,0,100); | |
1378 | hframe->SetStats(0); | |
1379 | if(fMultSoft) | |
1380 | hframe->SetXTitle("#hat{q} [GeV^{2}/fm]"); | |
1381 | else | |
1382 | hframe->SetXTitle("#mu [GeV]"); | |
1383 | hframe->SetYTitle("<E_{loss}> [GeV]"); | |
1384 | hframe->Draw(); | |
1385 | ||
1386 | TGraph *gq=new TGraph(20); | |
1387 | Int_t i=0; | |
1388 | for(Double_t v=0.05;v<=5.05;v+=0.25){ | |
1389 | TH1F *dummy=ComputeELossHisto(1,v,len,e); | |
1390 | Double_t avgloss=dummy->GetMean(); | |
1391 | gq->SetPoint(i,v,avgloss);i++; | |
1392 | delete dummy; | |
1393 | } | |
1394 | gq->SetMarkerStyle(20); | |
1395 | gq->Draw("pl"); | |
1396 | ||
1397 | TGraph *gg=new TGraph(20); | |
1398 | i=0; | |
1399 | for(Double_t v=0.05;v<=5.05;v+=0.25){ | |
1400 | TH1F *dummy=ComputeELossHisto(2,v,len,e); | |
1401 | Double_t avgloss=dummy->GetMean(); | |
1402 | gg->SetPoint(i,v,avgloss);i++; | |
1403 | delete dummy; | |
1404 | } | |
1405 | gg->SetMarkerStyle(24); | |
1406 | gg->Draw("pl"); | |
1407 | ||
1408 | TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); | |
1409 | l1a->SetFillStyle(0); | |
1410 | l1a->SetBorderSize(0); | |
1411 | Char_t label[100]; | |
1412 | sprintf(label,"L = %d fm",len); | |
1413 | l1a->AddEntry(gq,label,""); | |
1414 | l1a->AddEntry(gq,"quark","pl"); | |
1415 | l1a->AddEntry(gg,"gluon","pl"); | |
1416 | l1a->Draw(); | |
1417 | ||
1418 | c->Update(); | |
1419 | } | |
1420 | ||
1421 | void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const | |
1422 | { | |
1423 | //plot average energy loss for given | |
1424 | //length distribution and parton energy | |
1425 | ||
1426 | if(!fTablesLoaded){ | |
1427 | Error("CalcMult","Tables are not loaded."); | |
1428 | return; | |
1429 | } | |
1430 | ||
1431 | Char_t title[1024]; | |
1432 | Char_t name[1024]; | |
1433 | if(fMultSoft){ | |
1434 | sprintf(title,"Average Loss for Multiple Scattering"); | |
1435 | sprintf(name,"cavgelossms2"); | |
1436 | } else { | |
1437 | sprintf(title,"Average Loss for Single Hard Scattering"); | |
1438 | sprintf(name,"cavgelosssh2"); | |
1439 | } | |
1440 | ||
1441 | TCanvas *c = new TCanvas(name,title,0,0,500,400); | |
1442 | c->cd(); | |
1443 | TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100); | |
1444 | hframe->SetStats(0); | |
1445 | if(fMultSoft) | |
1446 | hframe->SetXTitle("#hat{q} [GeV^{2}/fm]"); | |
1447 | else | |
1448 | hframe->SetXTitle("#mu [GeV]"); | |
1449 | hframe->SetYTitle("<E_{loss}> [GeV]"); | |
1450 | hframe->Draw(); | |
1451 | ||
1452 | TGraph *gq=new TGraph(20); | |
1453 | Int_t i=0; | |
1454 | for(Double_t v=0.05;v<=5.05;v+=0.25){ | |
1455 | TH1F *dummy=ComputeELossHisto(1,v,hEll,e); | |
1456 | Double_t avgloss=dummy->GetMean(); | |
1457 | gq->SetPoint(i,v,avgloss);i++; | |
1458 | delete dummy; | |
1459 | } | |
1460 | gq->SetMarkerStyle(20); | |
1461 | gq->Draw("pl"); | |
1462 | ||
1463 | TGraph *gg=new TGraph(20); | |
1464 | i=0; | |
1465 | for(Double_t v=0.05;v<=5.05;v+=0.25){ | |
1466 | TH1F *dummy=ComputeELossHisto(2,v,hEll,e); | |
1467 | Double_t avgloss=dummy->GetMean(); | |
1468 | gg->SetPoint(i,v,avgloss);i++; | |
1469 | delete dummy; | |
1470 | } | |
1471 | gg->SetMarkerStyle(24); | |
1472 | gg->Draw("pl"); | |
1473 | ||
1474 | TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); | |
1475 | l1a->SetFillStyle(0); | |
1476 | l1a->SetBorderSize(0); | |
1477 | Char_t label[100]; | |
1478 | sprintf(label,"<L> = %.2f fm",hEll->GetMean()); | |
1479 | l1a->AddEntry(gq,label,""); | |
1480 | l1a->AddEntry(gq,"quark","pl"); | |
1481 | l1a->AddEntry(gg,"gluon","pl"); | |
1482 | l1a->Draw(); | |
1483 | ||
1484 | c->Update(); | |
1485 | } | |
1486 | ||
1487 | void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Int_t len) const | |
1488 | { | |
1489 | //plot relative energy loss for given | |
1490 | //length and parton energy versus pt | |
1491 | ||
1492 | if(!fTablesLoaded){ | |
1493 | Error("CalcMult","Tables are not loaded."); | |
1494 | return; | |
1495 | } | |
1496 | ||
1497 | Char_t title[1024]; | |
1498 | Char_t name[1024]; | |
1499 | if(fMultSoft){ | |
1500 | sprintf(title,"Relative Loss for Multiple Scattering"); | |
1501 | sprintf(name,"cavgelossvsptms"); | |
1502 | } else { | |
1503 | sprintf(title,"Relative Loss for Single Hard Scattering"); | |
1504 | sprintf(name,"cavgelossvsptsh"); | |
1505 | } | |
1506 | ||
1507 | TCanvas *c = new TCanvas(name,title,0,0,500,400); | |
1508 | c->cd(); | |
1509 | TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1); | |
1510 | hframe->SetStats(0); | |
1511 | hframe->SetXTitle("p_{T} [GeV]"); | |
1512 | hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]"); | |
1513 | hframe->Draw(); | |
1514 | ||
1515 | TGraph *gq=new TGraph(40); | |
1516 | Int_t i=0; | |
1517 | for(Double_t pt=2.5;pt<=100.05;pt+=2.5){ | |
1518 | TH1F *dummy=ComputeELossHisto(1,medval,len,pt); | |
1519 | Double_t avgloss=dummy->GetMean(); | |
1520 | gq->SetPoint(i,pt,avgloss/pt);i++; | |
1521 | delete dummy; | |
1522 | } | |
1523 | gq->SetMarkerStyle(20); | |
1524 | gq->Draw("pl"); | |
1525 | ||
1526 | TGraph *gg=new TGraph(40); | |
1527 | i=0; | |
1528 | for(Double_t pt=2.5;pt<=100.05;pt+=2.5){ | |
1529 | TH1F *dummy=ComputeELossHisto(2,medval,len,pt); | |
1530 | Double_t avgloss=dummy->GetMean(); | |
1531 | gg->SetPoint(i,pt,avgloss/pt);i++; | |
1532 | delete dummy; | |
1533 | } | |
1534 | gg->SetMarkerStyle(24); | |
1535 | gg->Draw("pl"); | |
1536 | ||
1537 | TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); | |
1538 | l1a->SetFillStyle(0); | |
1539 | l1a->SetBorderSize(0); | |
1540 | Char_t label[100]; | |
1541 | sprintf(label,"L = %d fm",len); | |
1542 | l1a->AddEntry(gq,label,""); | |
1543 | l1a->AddEntry(gq,"quark","pl"); | |
1544 | l1a->AddEntry(gg,"gluon","pl"); | |
1545 | l1a->Draw(); | |
1546 | ||
1547 | c->Update(); | |
1548 | } | |
1549 | ||
1550 | void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const | |
1551 | { | |
1552 | //plot relative energy loss for given | |
1553 | //length distribution and parton energy versus pt | |
1554 | ||
1555 | if(!fTablesLoaded){ | |
1556 | Error("CalcMult","Tables are not loaded."); | |
1557 | return; | |
1558 | } | |
1559 | ||
1560 | Char_t title[1024]; | |
1561 | Char_t name[1024]; | |
1562 | if(fMultSoft){ | |
1563 | sprintf(title,"Relative Loss for Multiple Scattering"); | |
1564 | sprintf(name,"cavgelossvsptms2"); | |
1565 | } else { | |
1566 | sprintf(title,"Relative Loss for Single Hard Scattering"); | |
1567 | sprintf(name,"cavgelossvsptsh2"); | |
1568 | } | |
1569 | TCanvas *c = new TCanvas(name,title,0,0,500,400); | |
1570 | c->cd(); | |
1571 | TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1); | |
1572 | hframe->SetStats(0); | |
1573 | hframe->SetXTitle("p_{T} [GeV]"); | |
1574 | hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]"); | |
1575 | hframe->Draw(); | |
1576 | ||
1577 | TGraph *gq=new TGraph(40); | |
1578 | Int_t i=0; | |
1579 | for(Double_t pt=2.5;pt<=100.05;pt+=2.5){ | |
1580 | TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt); | |
1581 | Double_t avgloss=dummy->GetMean(); | |
1582 | gq->SetPoint(i,pt,avgloss/pt);i++; | |
1583 | delete dummy; | |
1584 | } | |
1585 | gq->SetMarkerStyle(20); | |
1586 | gq->Draw("pl"); | |
1587 | ||
1588 | TGraph *gg=new TGraph(40); | |
1589 | i=0; | |
1590 | for(Double_t pt=2.5;pt<=100.05;pt+=2.5){ | |
1591 | TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt); | |
1592 | Double_t avgloss=dummy->GetMean(); | |
1593 | gg->SetPoint(i,pt,avgloss/pt);i++; | |
1594 | delete dummy; | |
1595 | } | |
1596 | gg->SetMarkerStyle(24); | |
1597 | gg->Draw("pl"); | |
1598 | ||
1599 | TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); | |
1600 | l1a->SetFillStyle(0); | |
1601 | l1a->SetBorderSize(0); | |
1602 | Char_t label[100]; | |
1603 | sprintf(label,"<L> = %.2f fm",hEll->GetMean()); | |
1604 | l1a->AddEntry(gq,label,""); | |
1605 | l1a->AddEntry(gq,"quark","pl"); | |
1606 | l1a->AddEntry(gg,"gluon","pl"); | |
1607 | l1a->Draw(); | |
1608 | ||
1609 | c->Update(); | |
1610 | } | |
1611 |