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