Cosmetics and new GetLengthDef() method
[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
2552c51a 1010 fHisto->SetBinContent(bin,continuous);
8ab8044e 1011#if 1 /* fast new version*/
1012 const Int_t kbinmax=fHisto->FindBin(e/wc);
1013 if(fMultSoft) {
1014 for(Int_t bin=2; bin<=kbinmax; 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<=kbinmax; bin++) {
1021 xxxx = fHisto->GetBinCenter(bin);
1022 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1023 fHisto->SetBinContent(bin,continuous);
1024 }
1025 }
2552c51a 1026
8ab8044e 1027 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1028
1029 Double_t val=discrete*fgkBins/fgkMaxBin;
1030 fHisto->Fill(0.,val);
1031
1032 if(fECMethod==kReweight)
1033 fHisto->SetBinContent(kbinmax+1,0);
1034 else
1035 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1036
1037 for(Int_t bin=kbinmax+2; bin<=fgkBins; bin++) {
1038 fHisto->SetBinContent(bin,0);
1039 }
1040
1041#else
2552c51a 1042 if(fMultSoft) {
1043 for(Int_t bin=2; bin<=fgkBins; bin++) {
1044 xxxx = fHisto->GetBinCenter(bin);
1045 CalcMult(ipart,R,xxxx,continuous,discrete);
1046 fHisto->SetBinContent(bin,continuous);
1047 }
1048 } else {
1049 for(Int_t bin=2; bin<=fgkBins; bin++) {
1050 xxxx = fHisto->GetBinCenter(bin);
1051 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1052 fHisto->SetBinContent(bin,continuous);
1053 }
1054 }
1055
1056 // add discrete part to distribution
1057 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
1058 fHisto->Fill(0.,val);
1059
1060 if(fECMethod==kReweight) {
1061 for(Int_t bin=fHisto->FindBin(e/wc)+1; bin<=fgkBins; bin++) {
1062 fHisto->SetBinContent(bin,0);
1063 }
1064 }
8ab8044e 1065#endif
2552c51a 1066
1067 Double_t ret=fHisto->GetRandom()*wc;
1068 if(ret>e) return e;
1069 return ret;
1070}
1071
1072Double_t AliQuenchingWeights::CalcQuenchedEnergyKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
1073{
1074 //return quenched parton energy (fast method)
1075 //for given parton type, I0 and I1 value and energy
1076
1077 Double_t loss=GetELossRandomKFast(ipart,I0,I1,e);
1078 return e-loss;
1079}
1080
b6d061b7 1081Int_t AliQuenchingWeights::SampleEnergyLoss()
1082{
1083 // Has to be called to fill the histograms
1084 //
1085 // For stored values fQTransport loop over
1086 // particle type and length = 1 to fMaxLength (fm)
1087 // to fill energy loss histos
1088 //
1089 // Take histogram of continuous weights
1090 // Take discrete_weight
1091 // If discrete_weight > 1, put all channels to 0, except channel 1
1092 // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral
1093
1094 // read-in data before first call
1095 if(!fTablesLoaded){
7258586f 1096 Error("SampleEnergyLoss","Tables are not loaded.");
b6d061b7 1097 return -1;
1098 }
1099
1100 if(fMultSoft) {
1101 Int_t lmax=CalcLengthMax(fQTransport);
1102 if(fLengthMax>lmax){
7a76a12e 1103 Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
b6d061b7 1104 fLengthMax=lmax;
1105 }
1106 } else {
1107 Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested.");
1108 }
1109
1110 Reset();
1111 fHistos=new TH1F**[2];
9d851d20 1112 fHistos[0]=new TH1F*[4*fLengthMax];
1113 fHistos[1]=new TH1F*[4*fLengthMax];
b6d061b7 1114 fLengthMaxOld=fLengthMax; //remember old value in case
1115 //user wants to reset
1116
1117 Int_t medvalue=0;
1118 Char_t meddesc[100];
1119 if(fMultSoft) {
1120 medvalue=(Int_t)(fQTransport*1000.);
1121 sprintf(meddesc,"MS");
1122 } else {
1123 medvalue=(Int_t)(fMu*1000.);
1124 sprintf(meddesc,"SH");
1125 }
1126
1127 for(Int_t ipart=1;ipart<=2;ipart++){
9d851d20 1128 for(Int_t l=1;l<=4*fLengthMax;l++){
b6d061b7 1129 Char_t hname[100];
1130 sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
9d851d20 1131 Double_t len=l/4.;
7258586f 1132 Double_t wc = CalcWC(len);
2552c51a 1133 fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgkBins,0.,fgkMaxBin*wc);
b6d061b7 1134 fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
1135 fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
1136 fHistos[ipart-1][l-1]->SetLineColor(4);
1137
7258586f 1138 Double_t rrrr = CalcR(wc,len);
b6d061b7 1139 Double_t discrete=0.;
1140 // loop on histogram channels
2552c51a 1141 for(Int_t bin=1; bin<=fgkBins; bin++) {
b6d061b7 1142 Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
1143 Double_t continuous;
7258586f 1144 if(fMultSoft)
1145 CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1146 else
1147 CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
b6d061b7 1148 fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
1149 }
1150 // add discrete part to distribution
1151 if(discrete>=1.)
2552c51a 1152 for(Int_t bin=2;bin<=fgkBins;bin++)
b6d061b7 1153 fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
1154 else {
2552c51a 1155 Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgkBins);
b6d061b7 1156 fHistos[ipart-1][l-1]->Fill(0.,val);
1157 }
2552c51a 1158 Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins);
b6d061b7 1159 fHistos[ipart-1][l-1]->Scale(1./hint);
1160 }
1161 }
1162 return 0;
1163}
1164
7258586f 1165Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t R)
1166{
1167 // Sample energy loss directly for one particle type
1168 // choses R (safe it and keep it until next call of function)
1169
1170 // read-in data before first call
1171 if(!fTablesLoaded){
1172 Error("SampleEnergyLoss","Tables are not loaded.");
1173 return -1;
1174 }
1175
1176 Double_t discrete=0.;
1177 Double_t continuous=0;;
1178 Int_t bin=1;
1179 Double_t xxxx = fHisto->GetBinCenter(bin);
1180 if(fMultSoft)
1181 CalcMult(ipart,R,xxxx,continuous,discrete);
1182 else
1183 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1184
1185 if(discrete>=1.) {
1186 fHisto->SetBinContent(1,1.);
2552c51a 1187 for(Int_t bin=2;bin<=fgkBins;bin++)
7258586f 1188 fHisto->SetBinContent(bin,0.);
1189 return 0;
1190 }
1191
1192 fHisto->SetBinContent(bin,continuous);
2552c51a 1193 for(Int_t bin=2; bin<=fgkBins; bin++) {
7258586f 1194 xxxx = fHisto->GetBinCenter(bin);
1195 if(fMultSoft)
1196 CalcMult(ipart,R,xxxx,continuous,discrete);
1197 else
1198 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1199 fHisto->SetBinContent(bin,continuous);
1200 }
8ab8044e 1201
2552c51a 1202 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
7258586f 1203 fHisto->Fill(0.,val);
2552c51a 1204 Double_t hint=fHisto->Integral(1,fgkBins);
7258586f 1205 if(hint!=0)
1206 fHisto->Scale(1./hint);
1207 else {
2552c51a 1208 //cout << discrete << " " << hint << " " << continuous << endl;
7258586f 1209 return -1;
1210 }
1211 return 0;
1212}
1213
1214const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) const
b6d061b7 1215{
7a76a12e 1216 //return quenching histograms
1217 //for ipart and length
1218
b6d061b7 1219 if(!fHistos){
1220 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
1221 return 0;
1222 }
1223 if((ipart<1) || (ipart>2)) {
1224 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
1225 return 0;
1226 }
1227
7258586f 1228 Int_t l=GetIndex(length);
b6d061b7 1229 if(l<=0) return 0;
b6d061b7 1230 return fHistos[ipart-1][l-1];
1231}
1232
1233TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const
1234{
1235 // ipart = 1 for quark, 2 for gluon
1236 // medval a) qtransp = transport coefficient (GeV^2/fm)
1237 // b) mu = Debye mass (GeV)
1238 // length = path length in medium (fm)
1239 // Get from SW tables:
7258586f 1240 // - continuous weight, as a function of dE/wc
b6d061b7 1241
1242 Double_t wc = 0;
1243 Char_t meddesc[100];
1244 if(fMultSoft) {
1245 wc=CalcWC(medval,length);
1246 sprintf(meddesc,"MS");
1247 } else {
1248 wc=CalcWCbar(medval,length);
1249 sprintf(meddesc,"SH");
1250 }
1251
1252 Char_t hname[100];
1253 sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart,
1254 (Int_t)(medval*1000.),(Int_t)length);
1255
2552c51a 1256 TH1F *hist = new TH1F("hist",hname,fgkBins,0.,fgkMaxBin*wc);
b6d061b7 1257 hist->SetXTitle("#Delta E [GeV]");
1258 hist->SetYTitle("p(#Delta E)");
1259 hist->SetLineColor(4);
1260
1261 Double_t rrrr = CalcR(wc,length);
7258586f 1262 //loop on histogram channels
2552c51a 1263 for(Int_t bin=1; bin<=fgkBins; bin++) {
b6d061b7 1264 Double_t xxxx = hist->GetBinCenter(bin)/wc;
1265 Double_t continuous,discrete;
1266 Int_t ret=0;
1267 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1268 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1269 if(ret!=0){
1270 delete hist;
1271 return 0;
1272 };
1273 hist->SetBinContent(bin,continuous);
1274 }
1275 return hist;
1276}
1277
1278TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const
1279{
1280 // ipart = 1 for quark, 2 for gluon
1281 // medval a) qtransp = transport coefficient (GeV^2/fm)
1282 // b) mu = Debye mass (GeV)
1283 // length = path length in medium (fm)
1284 // Get from SW tables:
7258586f 1285 // - continuous weight, as a function of dE/wc
b6d061b7 1286
1287 Double_t wc = 0;
1288 Char_t meddesc[100];
1289 if(fMultSoft) {
1290 wc=CalcWC(medval,length);
1291 sprintf(meddesc,"MS");
1292 } else {
1293 wc=CalcWCbar(medval,length);
1294 sprintf(meddesc,"SH");
1295 }
1296
1297 Char_t hname[100];
1298 sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart,
1299 (Int_t)(medval*1000.),(Int_t)length);
1300
2552c51a 1301 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
b6d061b7 1302 histx->SetXTitle("x = #Delta E/#omega_{c}");
1303 if(fMultSoft)
1304 histx->SetYTitle("p(#Delta E/#omega_{c})");
1305 else
1306 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1307 histx->SetLineColor(4);
1308
1309 Double_t rrrr = CalcR(wc,length);
7258586f 1310 //loop on histogram channels
2552c51a 1311 for(Int_t bin=1; bin<=fgkBins; bin++) {
b6d061b7 1312 Double_t xxxx = histx->GetBinCenter(bin);
1313 Double_t continuous,discrete;
1314 Int_t ret=0;
1315 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1316 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1317 if(ret!=0){
1318 delete histx;
1319 return 0;
1320 };
1321 histx->SetBinContent(bin,continuous);
1322 }
1323 return histx;
1324}
1325
7258586f 1326TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t R) const
1327{
1328 // compute P(E) distribution for
1329 // given ipart = 1 for quark, 2 for gluon
1330 // and R
1331
1332 Char_t meddesc[100];
1333 if(fMultSoft) {
1334 sprintf(meddesc,"MS");
1335 } else {
1336 sprintf(meddesc,"SH");
1337 }
1338
1339 Char_t hname[100];
1340 sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,R);
2552c51a 1341 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
7258586f 1342 histx->SetXTitle("x = #Delta E/#omega_{c}");
1343 if(fMultSoft)
1344 histx->SetYTitle("p(#Delta E/#omega_{c})");
1345 else
1346 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1347 histx->SetLineColor(4);
1348
1349 Double_t rrrr = R;
1350 Double_t continuous=0.,discrete=0.;
1351 //loop on histogram channels
2552c51a 1352 for(Int_t bin=1; bin<=fgkBins; bin++) {
7258586f 1353 Double_t xxxx = histx->GetBinCenter(bin);
1354 Int_t ret=0;
1355 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1356 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1357 if(ret!=0){
1358 delete histx;
1359 return 0;
1360 };
1361 histx->SetBinContent(bin,continuous);
1362 }
1363
1364 //add discrete part to distribution
1365 if(discrete>=1.)
2552c51a 1366 for(Int_t bin=2;bin<=fgkBins;bin++)
7258586f 1367 histx->SetBinContent(bin,0.);
1368 else {
2552c51a 1369 Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins);
7258586f 1370 histx->Fill(0.,val);
1371 }
2552c51a 1372 Double_t hint=histx->Integral(1,fgkBins);
7258586f 1373 if(hint!=0) histx->Scale(1./hint);
1374
1375 return histx;
1376}
1377
e99e3ed5 1378TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
b6d061b7 1379{
7258586f 1380 // compute energy loss histogram for
1381 // parton type, medium value, length and energy
7a76a12e 1382
b6d061b7 1383 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1384 if(fMultSoft){
1385 dummy->SetQTransport(medval);
1386 dummy->InitMult();
1387 } else {
1388 dummy->SetMu(medval);
1389 dummy->InitSingleHard();
1390 }
1391 dummy->SampleEnergyLoss();
1392
1393 Char_t name[100];
1394 Char_t hname[100];
1395 if(ipart==1){
1396 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1397 sprintf(hname,"hLossQuarks");
1398 } else {
1399 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1400 sprintf(hname,"hLossGluons");
1401 }
1402
1403 TH1F *h = new TH1F(hname,name,250,0,250);
1404 for(Int_t i=0;i<100000;i++){
1405 //if(i % 1000 == 0) cout << "." << flush;
1406 Double_t loss=dummy->GetELossRandom(ipart,l,e);
1407 h->Fill(loss);
1408 }
b6d061b7 1409 h->SetStats(kTRUE);
b6d061b7 1410 delete dummy;
1411 return h;
1412}
1413
e99e3ed5 1414TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
b6d061b7 1415{
7258586f 1416 // compute energy loss histogram for
1417 // parton type, medium value,
1418 // length distribution and energy
7a76a12e 1419
b6d061b7 1420 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1421 if(fMultSoft){
1422 dummy->SetQTransport(medval);
1423 dummy->InitMult();
1424 } else {
1425 dummy->SetMu(medval);
1426 dummy->InitSingleHard();
1427 }
1428 dummy->SampleEnergyLoss();
1429
1430 Char_t name[100];
1431 Char_t hname[100];
1432 if(ipart==1){
1433 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1434 sprintf(hname,"hLossQuarks");
1435 } else {
1436 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1437 sprintf(hname,"hLossGluons");
1438 }
1439
1440 TH1F *h = new TH1F(hname,name,250,0,250);
1441 for(Int_t i=0;i<100000;i++){
1442 //if(i % 1000 == 0) cout << "." << flush;
1443 Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
1444 h->Fill(loss);
1445 }
b6d061b7 1446 h->SetStats(kTRUE);
b6d061b7 1447 delete dummy;
1448 return h;
1449}
1450
7258586f 1451TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t R) const
1452{
1453 // compute energy loss histogram for
1454 // parton type and given R
1455
1456 TH1F *dummy = ComputeQWHistoX(ipart,R);
1457 if(!dummy) return 0;
1458
1459 Char_t hname[100];
1460 sprintf(hname,"hELossHistox_%d_%.2f",ipart,R);
2552c51a 1461 TH1F *histx = new TH1F("histxr",hname,fgkBins,0.,fgkMaxBin);
7258586f 1462 for(Int_t i=0;i<100000;i++){
1463 //if(i % 1000 == 0) cout << "." << flush;
1464 Double_t loss=dummy->GetRandom();
1465 histx->Fill(loss);
1466 }
1467 delete dummy;
1468 return histx;
1469}
1470
1471Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,Double_t l) const
1472{
1473 // compute average energy loss for
1474 // parton type, medium value, length and energy
1475
1476 TH1F *dummy = ComputeELossHisto(ipart,medval,l);
1477 if(!dummy) return 0;
1478 Double_t ret=dummy->GetMean();
1479 delete dummy;
1480 return ret;
1481}
1482
1483Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEll) const
1484{
1485 // compute average energy loss for
1486 // parton type, medium value,
1487 // length distribution and energy
1488
1489 TH1F *dummy = ComputeELossHisto(ipart,medval,hEll);
1490 if(!dummy) return 0;
1491 Double_t ret=dummy->GetMean();
1492 delete dummy;
1493 return ret;
1494}
1495
1496Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t R) const
1497{
1498 // compute average energy loss over wc
1499 // for parton type and given R
1500
1501 TH1F *dummy = ComputeELossHisto(ipart,R);
1502 if(!dummy) return 0;
1503 Double_t ret=dummy->GetMean();
1504 delete dummy;
1505 return ret;
1506}
1507
1508void AliQuenchingWeights::PlotDiscreteWeights(Double_t len) const
b6d061b7 1509{
7258586f 1510 // plot discrete weights for given length
7a76a12e 1511
b6d061b7 1512 TCanvas *c;
1513 if(fMultSoft)
1514 c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
1515 else
1516 c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
1517 c->cd();
1518
1519 TH2F *hframe = new TH2F("hdisc","",2,0,5.1,2,0,1);
1520 hframe->SetStats(0);
1521 if(fMultSoft)
1522 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1523 else
1524 hframe->SetXTitle("#mu [GeV]");
1525 hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
1526 hframe->Draw();
1527
1528 TGraph *gq=new TGraph(20);
1529 Int_t i=0;
1530 if(fMultSoft) {
1531 for(Double_t q=0.05;q<=5.05;q+=0.25){
1532 Double_t disc,cont;
7258586f 1533 CalcMult(1,1.0,q,len,cont,disc);
b6d061b7 1534 gq->SetPoint(i,q,disc);i++;
1535 }
1536 } else {
1537 for(Double_t m=0.05;m<=5.05;m+=0.25){
1538 Double_t disc,cont;
7258586f 1539 CalcSingleHard(1,1.0,m,len,cont, disc);
b6d061b7 1540 gq->SetPoint(i,m,disc);i++;
1541 }
1542 }
1543 gq->SetMarkerStyle(20);
1544 gq->Draw("pl");
1545
1546 TGraph *gg=new TGraph(20);
1547 i=0;
1548 if(fMultSoft){
1549 for(Double_t q=0.05;q<=5.05;q+=0.25){
1550 Double_t disc,cont;
7258586f 1551 CalcMult(2,1.0,q,len,cont,disc);
b6d061b7 1552 gg->SetPoint(i,q,disc);i++;
1553 }
1554 } else {
1555 for(Double_t m=0.05;m<=5.05;m+=0.25){
1556 Double_t disc,cont;
7258586f 1557 CalcSingleHard(2,1.0,m,len,cont,disc);
b6d061b7 1558 gg->SetPoint(i,m,disc);i++;
1559 }
1560 }
1561 gg->SetMarkerStyle(24);
1562 gg->Draw("pl");
1563
1564 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1565 l1a->SetFillStyle(0);
1566 l1a->SetBorderSize(0);
1567 Char_t label[100];
7258586f 1568 sprintf(label,"L = %.1f fm",len);
b6d061b7 1569 l1a->AddEntry(gq,label,"");
1570 l1a->AddEntry(gq,"quark","pl");
1571 l1a->AddEntry(gg,"gluon","pl");
1572 l1a->Draw();
1573
1574 c->Update();
1575}
1576
7258586f 1577void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
b6d061b7 1578{
7258586f 1579 // plot continous weights for
1580 // given parton type and length
7a76a12e 1581
b6d061b7 1582 Float_t medvals[3];
1583 Char_t title[1024];
1584 Char_t name[1024];
1585 if(fMultSoft) {
1586 if(itype==1)
1587 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1588 else if(itype==2)
1589 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1590 else return;
1591 medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1592 sprintf(name,"ccont-ms-%d",itype);
1593 } else {
1594 if(itype==1)
1595 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1596 else if(itype==2)
1597 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1598 else return;
1599 medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1600 sprintf(name,"ccont-ms-%d",itype);
1601 }
1602
1603 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1604 c->cd();
1605 TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1606 h1->SetName("h1");
1607 h1->SetTitle(title);
1608 h1->SetStats(0);
1609 h1->SetLineColor(1);
1610 h1->DrawCopy();
1611 TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1612 h2->SetName("h2");
1613 h2->SetLineColor(2);
1614 h2->DrawCopy("SAME");
1615 TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1616 h3->SetName("h3");
1617 h3->SetLineColor(3);
1618 h3->DrawCopy("SAME");
1619
1620 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1621 l1a->SetFillStyle(0);
1622 l1a->SetBorderSize(0);
1623 Char_t label[100];
7258586f 1624 sprintf(label,"L = %.1f fm",ell);
b6d061b7 1625 l1a->AddEntry(h1,label,"");
1626 if(fMultSoft) {
1627 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1628 l1a->AddEntry(h1,label,"pl");
1629 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1630 l1a->AddEntry(h2,label,"pl");
1631 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1632 l1a->AddEntry(h3,label,"pl");
1633 } else {
1634 sprintf(label,"#mu = %.1f GeV",medvals[0]);
1635 l1a->AddEntry(h1,label,"pl");
1636 sprintf(label,"#mu = %.1f GeV",medvals[1]);
1637 l1a->AddEntry(h2,label,"pl");
1638 sprintf(label,"#mu = %.1f GeV",medvals[2]);
1639 l1a->AddEntry(h3,label,"pl");
1640 }
1641 l1a->Draw();
1642
1643 c->Update();
1644}
1645
7258586f 1646void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
b6d061b7 1647{
7258586f 1648 // plot continous weights for
1649 // given parton type and medium value
7a76a12e 1650
b6d061b7 1651 Char_t title[1024];
1652 Char_t name[1024];
1653 if(fMultSoft) {
1654 if(itype==1)
1655 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1656 else if(itype==2)
1657 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1658 else return;
1659 sprintf(name,"ccont2-ms-%d",itype);
1660 } else {
1661 if(itype==1)
1662 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1663 else if(itype==2)
1664 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1665 else return;
1666 sprintf(name,"ccont2-sh-%d",itype);
1667 }
1668 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1669 c->cd();
1670 TH1F *h1=ComputeQWHisto(itype,medval,8);
1671 h1->SetName("h1");
1672 h1->SetTitle(title);
1673 h1->SetStats(0);
1674 h1->SetLineColor(1);
1675 h1->DrawCopy();
1676 TH1F *h2=ComputeQWHisto(itype,medval,5);
1677 h2->SetName("h2");
1678 h2->SetLineColor(2);
1679 h2->DrawCopy("SAME");
1680 TH1F *h3=ComputeQWHisto(itype,medval,2);
1681 h3->SetName("h3");
1682 h3->SetLineColor(3);
1683 h3->DrawCopy("SAME");
1684
1685 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1686 l1a->SetFillStyle(0);
1687 l1a->SetBorderSize(0);
1688 Char_t label[100];
1689 if(fMultSoft)
1690 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval);
1691 else
1692 sprintf(label,"#mu = %.1f GeV",medval);
1693
1694 l1a->AddEntry(h1,label,"");
1695 l1a->AddEntry(h1,"L = 8 fm","pl");
1696 l1a->AddEntry(h2,"L = 5 fm","pl");
1697 l1a->AddEntry(h3,"L = 2 fm","pl");
1698 l1a->Draw();
1699
1700 c->Update();
1701}
1702
7258586f 1703void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t e) const
b6d061b7 1704{
7258586f 1705 // plot average energy loss for given length
1706 // and parton energy
7a76a12e 1707
b6d061b7 1708 if(!fTablesLoaded){
7258586f 1709 Error("PlotAvgELoss","Tables are not loaded.");
b6d061b7 1710 return;
1711 }
1712
1713 Char_t title[1024];
1714 Char_t name[1024];
1715 if(fMultSoft){
7258586f 1716 sprintf(title,"Average Energy Loss for Multiple Scattering");
b6d061b7 1717 sprintf(name,"cavgelossms");
1718 } else {
7258586f 1719 sprintf(title,"Average Energy Loss for Single Hard Scattering");
b6d061b7 1720 sprintf(name,"cavgelosssh");
1721 }
1722
1723 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1724 c->cd();
1725 TH2F *hframe = new TH2F("avgloss",title,2,0,5.1,2,0,100);
1726 hframe->SetStats(0);
1727 if(fMultSoft)
1728 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1729 else
1730 hframe->SetXTitle("#mu [GeV]");
1731 hframe->SetYTitle("<E_{loss}> [GeV]");
1732 hframe->Draw();
1733
1734 TGraph *gq=new TGraph(20);
1735 Int_t i=0;
1736 for(Double_t v=0.05;v<=5.05;v+=0.25){
1737 TH1F *dummy=ComputeELossHisto(1,v,len,e);
1738 Double_t avgloss=dummy->GetMean();
1739 gq->SetPoint(i,v,avgloss);i++;
1740 delete dummy;
1741 }
1742 gq->SetMarkerStyle(20);
1743 gq->Draw("pl");
1744
1745 TGraph *gg=new TGraph(20);
1746 i=0;
1747 for(Double_t v=0.05;v<=5.05;v+=0.25){
1748 TH1F *dummy=ComputeELossHisto(2,v,len,e);
1749 Double_t avgloss=dummy->GetMean();
1750 gg->SetPoint(i,v,avgloss);i++;
1751 delete dummy;
1752 }
1753 gg->SetMarkerStyle(24);
1754 gg->Draw("pl");
1755
7258586f 1756 TGraph *gratio=new TGraph(20);
1757 for(Int_t i=0;i<20;i++){
1758 Double_t x,y,x2,y2;
1759 gg->GetPoint(i,x,y);
1760 gq->GetPoint(i,x2,y2);
1761 if(y2>0)
1762 gratio->SetPoint(i,x,y/y2*10/2.25);
1763 else gratio->SetPoint(i,x,0);
1764 }
1765 gratio->SetLineStyle(4);
1766 gratio->Draw();
b6d061b7 1767 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1768 l1a->SetFillStyle(0);
1769 l1a->SetBorderSize(0);
1770 Char_t label[100];
7258586f 1771 sprintf(label,"L = %.1f fm",len);
b6d061b7 1772 l1a->AddEntry(gq,label,"");
1773 l1a->AddEntry(gq,"quark","pl");
1774 l1a->AddEntry(gg,"gluon","pl");
7258586f 1775 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
b6d061b7 1776 l1a->Draw();
1777
1778 c->Update();
1779}
1780
e99e3ed5 1781void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
b6d061b7 1782{
7258586f 1783 // plot average energy loss for given
1784 // length distribution and parton energy
7a76a12e 1785
b6d061b7 1786 if(!fTablesLoaded){
7258586f 1787 Error("PlotAvgELossVs","Tables are not loaded.");
b6d061b7 1788 return;
1789 }
1790
1791 Char_t title[1024];
1792 Char_t name[1024];
1793 if(fMultSoft){
7258586f 1794 sprintf(title,"Average Energy Loss for Multiple Scattering");
b6d061b7 1795 sprintf(name,"cavgelossms2");
1796 } else {
7258586f 1797 sprintf(title,"Average Energy Loss for Single Hard Scattering");
b6d061b7 1798 sprintf(name,"cavgelosssh2");
1799 }
1800
1801 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1802 c->cd();
1803 TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1804 hframe->SetStats(0);
1805 if(fMultSoft)
1806 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1807 else
1808 hframe->SetXTitle("#mu [GeV]");
1809 hframe->SetYTitle("<E_{loss}> [GeV]");
1810 hframe->Draw();
1811
1812 TGraph *gq=new TGraph(20);
1813 Int_t i=0;
1814 for(Double_t v=0.05;v<=5.05;v+=0.25){
1815 TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1816 Double_t avgloss=dummy->GetMean();
1817 gq->SetPoint(i,v,avgloss);i++;
1818 delete dummy;
1819 }
1820 gq->SetMarkerStyle(20);
1821 gq->Draw("pl");
1822
1823 TGraph *gg=new TGraph(20);
1824 i=0;
1825 for(Double_t v=0.05;v<=5.05;v+=0.25){
1826 TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1827 Double_t avgloss=dummy->GetMean();
1828 gg->SetPoint(i,v,avgloss);i++;
1829 delete dummy;
1830 }
1831 gg->SetMarkerStyle(24);
1832 gg->Draw("pl");
1833
7258586f 1834 TGraph *gratio=new TGraph(20);
1835 for(Int_t i=0;i<20;i++){
1836 Double_t x,y,x2,y2;
1837 gg->GetPoint(i,x,y);
1838 gq->GetPoint(i,x2,y2);
1839 if(y2>0)
1840 gratio->SetPoint(i,x,y/y2*10/2.25);
1841 else gratio->SetPoint(i,x,0);
1842 }
1843 gratio->SetLineStyle(4);
1844 gratio->Draw();
1845
b6d061b7 1846 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1847 l1a->SetFillStyle(0);
1848 l1a->SetBorderSize(0);
1849 Char_t label[100];
1850 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1851 l1a->AddEntry(gq,label,"");
1852 l1a->AddEntry(gq,"quark","pl");
1853 l1a->AddEntry(gg,"gluon","pl");
7258586f 1854 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
b6d061b7 1855 l1a->Draw();
1856
1857 c->Update();
1858}
1859
9d851d20 1860void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const
b6d061b7 1861{
7258586f 1862 // plot average energy loss versus ell
7a76a12e 1863
b6d061b7 1864 if(!fTablesLoaded){
7258586f 1865 Error("PlotAvgELossVsEll","Tables are not loaded.");
1866 return;
1867 }
1868
1869 Char_t title[1024];
1870 Char_t name[1024];
1871 Float_t medval;
1872 if(fMultSoft){
1873 sprintf(title,"Average Energy Loss for Multiple Scattering");
1874 sprintf(name,"cavgelosslms");
1875 medval=fQTransport;
1876 } else {
1877 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1878 sprintf(name,"cavgelosslsh");
1879 medval=fMu;
1880 }
1881
1882 TCanvas *c = new TCanvas(name,title,0,0,600,400);
1883 c->cd();
1884 TH2F *hframe = new TH2F("avglossell",title,2,0,fLengthMax,2,0,250);
1885 hframe->SetStats(0);
1886 hframe->SetXTitle("length [fm]");
1887 hframe->SetYTitle("<E_{loss}> [GeV]");
1888 hframe->Draw();
1889
1890 TGraph *gq=new TGraph((Int_t)fLengthMax*4);
1891 Int_t i=0;
1892 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
1893 TH1F *dummy=ComputeELossHisto(1,medval,v,e);
1894 Double_t avgloss=dummy->GetMean();
1895 gq->SetPoint(i,v,avgloss);i++;
1896 delete dummy;
1897 }
1898 gq->SetMarkerStyle(20);
1899 gq->Draw("pl");
1900
1901 TGraph *gg=new TGraph((Int_t)fLengthMax*4);
1902 i=0;
1903 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
1904 TH1F *dummy=ComputeELossHisto(2,medval,v,e);
1905 Double_t avgloss=dummy->GetMean();
1906 gg->SetPoint(i,v,avgloss);i++;
1907 delete dummy;
1908 }
1909 gg->SetMarkerStyle(24);
1910 gg->Draw("pl");
1911
1912 TGraph *gratio=new TGraph((Int_t)fLengthMax*4);
1913 for(Int_t i=0;i<=(Int_t)fLengthMax*4;i++){
1914 Double_t x,y,x2,y2;
1915 gg->GetPoint(i,x,y);
1916 gq->GetPoint(i,x2,y2);
1917 if(y2>0)
1918 gratio->SetPoint(i,x,y/y2*100/2.25);
1919 else gratio->SetPoint(i,x,0);
1920 }
1921 gratio->SetLineStyle(1);
1922 gratio->SetLineWidth(2);
1923 gratio->Draw();
1924 TLegend *l1a = new TLegend(0.15,0.65,.60,0.85);
1925 l1a->SetFillStyle(0);
1926 l1a->SetBorderSize(0);
1927 Char_t label[100];
1928 if(fMultSoft)
1929 sprintf(label,"#hat{q} = %.2f GeV^{2}/fm",medval);
1930 else
1931 sprintf(label,"#mu = %.2f GeV",medval);
1932 l1a->AddEntry(gq,label,"");
1933 l1a->AddEntry(gq,"quark","pl");
1934 l1a->AddEntry(gg,"gluon","pl");
1935 l1a->AddEntry(gratio,"gluon/quark/2.25*100","pl");
1936 l1a->Draw();
1937
1938 TF1 *f=new TF1("f","100",0,fLengthMax);
1939 f->SetLineStyle(4);
1940 f->SetLineWidth(1);
1941 f->Draw("same");
1942 c->Update();
1943}
1944
1945void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
1946{
1947 // plot relative energy loss for given
1948 // length and parton energy versus pt
1949
1950 if(!fTablesLoaded){
1951 Error("PlotAvgELossVsPt","Tables are not loaded.");
b6d061b7 1952 return;
1953 }
1954
1955 Char_t title[1024];
1956 Char_t name[1024];
1957 if(fMultSoft){
7258586f 1958 sprintf(title,"Relative Energy Loss for Multiple Scattering");
b6d061b7 1959 sprintf(name,"cavgelossvsptms");
1960 } else {
7258586f 1961 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
b6d061b7 1962 sprintf(name,"cavgelossvsptsh");
1963 }
1964
1965 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1966 c->cd();
1967 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
1968 hframe->SetStats(0);
1969 hframe->SetXTitle("p_{T} [GeV]");
1970 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
1971 hframe->Draw();
1972
1973 TGraph *gq=new TGraph(40);
1974 Int_t i=0;
1975 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1976 TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
1977 Double_t avgloss=dummy->GetMean();
1978 gq->SetPoint(i,pt,avgloss/pt);i++;
1979 delete dummy;
1980 }
1981 gq->SetMarkerStyle(20);
1982 gq->Draw("pl");
1983
1984 TGraph *gg=new TGraph(40);
1985 i=0;
1986 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1987 TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
1988 Double_t avgloss=dummy->GetMean();
1989 gg->SetPoint(i,pt,avgloss/pt);i++;
1990 delete dummy;
1991 }
1992 gg->SetMarkerStyle(24);
1993 gg->Draw("pl");
1994
1995 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1996 l1a->SetFillStyle(0);
1997 l1a->SetBorderSize(0);
1998 Char_t label[100];
7258586f 1999 sprintf(label,"L = %.1f fm",len);
b6d061b7 2000 l1a->AddEntry(gq,label,"");
2001 l1a->AddEntry(gq,"quark","pl");
2002 l1a->AddEntry(gg,"gluon","pl");
2003 l1a->Draw();
2004
2005 c->Update();
2006}
2007
2008void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
2009{
7258586f 2010 // plot relative energy loss for given
2011 // length distribution and parton energy versus pt
7a76a12e 2012
b6d061b7 2013 if(!fTablesLoaded){
7258586f 2014 Error("PlotAvgELossVsPt","Tables are not loaded.");
b6d061b7 2015 return;
2016 }
2017
2018 Char_t title[1024];
2019 Char_t name[1024];
2020 if(fMultSoft){
7258586f 2021 sprintf(title,"Relative Energy Loss for Multiple Scattering");
b6d061b7 2022 sprintf(name,"cavgelossvsptms2");
2023 } else {
7258586f 2024 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
b6d061b7 2025 sprintf(name,"cavgelossvsptsh2");
2026 }
2027 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2028 c->cd();
2029 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2030 hframe->SetStats(0);
2031 hframe->SetXTitle("p_{T} [GeV]");
2032 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2033 hframe->Draw();
2034
2035 TGraph *gq=new TGraph(40);
2036 Int_t i=0;
2037 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2038 TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
2039 Double_t avgloss=dummy->GetMean();
2040 gq->SetPoint(i,pt,avgloss/pt);i++;
2041 delete dummy;
2042 }
2043 gq->SetMarkerStyle(20);
2044 gq->Draw("pl");
2045
2046 TGraph *gg=new TGraph(40);
2047 i=0;
2048 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2049 TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
2050 Double_t avgloss=dummy->GetMean();
2051 gg->SetPoint(i,pt,avgloss/pt);i++;
2052 delete dummy;
2053 }
2054 gg->SetMarkerStyle(24);
2055 gg->Draw("pl");
2056
2057 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2058 l1a->SetFillStyle(0);
2059 l1a->SetBorderSize(0);
2060 Char_t label[100];
2061 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
2062 l1a->AddEntry(gq,label,"");
2063 l1a->AddEntry(gq,"quark","pl");
2064 l1a->AddEntry(gg,"gluon","pl");
2065 l1a->Draw();
2066
2067 c->Update();
2068}
2069
7258586f 2070Int_t AliQuenchingWeights::GetIndex(Double_t len) const
2071{
2072 if(len>fLengthMax) len=fLengthMax;
2073
9d851d20 2074 Int_t l=Int_t(len/0.25);
2075 if((len-l*0.5)>0.125) l++;
7258586f 2076 return l;
2077}
ea16e52f 2078