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