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