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