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