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