Made more robust
[u/mrichter/AliRoot.git] / FASTSIM / AliQuenchingWeights.cxx
CommitLineData
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
40ClassImp(AliQuenchingWeights)
41
42// conversion from fm to GeV^-1: 1 fm = fmGeV GeV^-1
43const Double_t AliQuenchingWeights::gkConvFmToInvGeV = 1./0.197;
44
45// maximum value of R
46const Double_t AliQuenchingWeights::gkRMax = 40000.;
47
48// counter for histogram labels
49Int_t AliQuenchingWeights::gCounter = 0;
50
51AliQuenchingWeights::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
65AliQuenchingWeights::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
81AliQuenchingWeights::~AliQuenchingWeights()
82{
83 Reset();
84}
85
86void 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
98void 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
107Int_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/*
174C***************************************************************************
175C Quenching Weights for Multiple Soft Scattering
176C February 10, 2003
177C
178C Refs:
179C
180C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
181C
182C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
183C
184C
185C This package contains quenching weights for gluon radiation in the
186C multiple soft scattering approximation.
187C
188C swqmult returns the quenching weight for a quark (ipart=1) or
189C a gluon (ipart=2) traversing a medium with transport coeficient q and
190C length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where
191C wc=0.5*q*L^2 and w is the energy radiated. The output values are
192C the continuous and discrete (prefactor of the delta function) parts
193C of the quenching weights.
194C
195C In order to use this routine, the files cont.all and disc.all need to be
196C in the working directory.
197C
198C An initialization of the tables is needed by doing call initmult before
199C using swqmult.
200C
201C Please, send us any comment:
202C
203C urs.wiedemann@cern.ch
204C carlos.salgado@cern.ch
205C
206C
207C-------------------------------------------------------------------
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
324Int_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
396Int_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/*
464C***************************************************************************
465C Quenching Weights for Single Hard Scattering
466C February 20, 2003
467C
468C Refs:
469C
470C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
471C
472C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
473C
474C
475C This package contains quenching weights for gluon radiation in the
476C single hard scattering approximation.
477C
478C swqlin returns the quenching weight for a quark (ipart=1) or
479C a gluon (ipart=2) traversing a medium with Debye screening mass mu and
480C length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where
481C wc=0.5*mu^2*L and w is the energy radiated. The output values are
482C the continuous and discrete (prefactor of the delta function) parts
483C of the quenching weights.
484C
485C In order to use this routine, the files contlin.all and disclin.all
486C need to be in the working directory.
487C
488C An initialization of the tables is needed by doing call initlin before
489C using swqlin.
490C
491C Please, send us any comment:
492C
493C urs.wiedemann@cern.ch
494C carlos.salgado@cern.ch
495C
496C
497C-------------------------------------------------------------------
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
615Int_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
687Int_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
697Int_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
707Double_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
717Double_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
760Double_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
e99e3ed5 769Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const
b6d061b7 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
779Double_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
788Int_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
869const 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
886TH1F* 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
935TH1F* 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
e99e3ed5 987TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
b6d061b7 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
e99e3ed5 1022TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
b6d061b7 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
1057void 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
1124void 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
1190void 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
e99e3ed5 1244void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e) const
b6d061b7 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
e99e3ed5 1307void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
b6d061b7 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
1370void 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
1430void 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