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