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