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