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