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