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