]> 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
Coding violations.
[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
b90de01a 933 Double_t r=CalcRk(I0,I1);
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 }
b90de01a 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
b90de01a 987 Double_t r=CalcRk(I0,I1);
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
b90de01a 997 return GetELossRandomKFastR(ipart,r,wc,e);
cc885e36 998}
999
b90de01a 1000Double_t AliQuenchingWeights::GetELossRandomKFastR(Int_t ipart, Double_t r, Double_t wc, Double_t e)
cc885e36 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
b90de01a 1010 if(r>=fgkRMax) {
1011 r=fgkRMax-1;
cc885e36 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)
b90de01a 1019 CalcMult(ipart,r,xxxx,continuous,discrete);
2552c51a 1020 else
b90de01a 1021 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
2552c51a 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);
b90de01a 1035 CalcMult(ipart,r,xxxx,continuous,discrete);
8ab8044e 1036 fHisto->SetBinContent(bin,continuous);
1037 }
1038 } else {
1039 for(Int_t bin=2; bin<=kbinmax; bin++) {
1040 xxxx = fHisto->GetBinCenter(bin);
b90de01a 1041 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
8ab8044e 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
b90de01a 1082 Double_t r=CalcRk(I0,I1);
1083 if(r<=0.){
b677de56 1084 return 1.;
1085 }
b90de01a 1086 return GetDiscreteWeightR(ipart,r);
b677de56 1087}
1088
b90de01a 1089Double_t AliQuenchingWeights::GetDiscreteWeightR(Int_t ipart, Double_t r)
b677de56 1090{
1091 // return discrete weight
1092
b90de01a 1093 if(r>=fgkRMax) {
1094 r=fgkRMax-1;
b677de56 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)
b90de01a 1102 CalcMult(ipart,r,xxxx,continuous,discrete);
b677de56 1103 else
b90de01a 1104 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
b677de56 1105 return discrete;
1106}
1107
b90de01a 1108void AliQuenchingWeights::GetZeroLossProb(Double_t &p,Double_t &prw,Double_t &prwcont,
b677de56 1109 Int_t ipart,Double_t I0,Double_t I1,Double_t e)
1110{
b90de01a 1111 //calculate the probabilty that there is no energy
1112 //loss for different ways of energy constraint
1113 p=1.;prw=1.;prwcont=1.;
1114 Double_t r=CalcRk(I0,I1);
1115 if(r<=0.){
b677de56 1116 return;
1117 }
1118 Double_t wc=CalcWCk(I1);
1119 if(wc<=0.){
1120 return;
1121 }
b90de01a 1122 GetZeroLossProbR(p,prw,prwcont,ipart,r,wc,e);
b677de56 1123}
1124
b90de01a 1125void AliQuenchingWeights::GetZeroLossProbR(Double_t &p,Double_t &prw,Double_t &prwcont,
1126 Int_t ipart, Double_t r,Double_t wc,Double_t e)
b677de56 1127{
b90de01a 1128 //calculate the probabilty that there is no energy
1129 //loss for different ways of energy constraint
1130 if(r>=fgkRMax) {
1131 r=fgkRMax-1;
b677de56 1132 }
1133
1134 Double_t discrete=0.;
1135 Double_t continuous=0.;
1136
1137 Int_t kbinmax=fHisto->FindBin(e/wc);
1138 if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
1139 if(fMultSoft) {
1140 for(Int_t bin=1; bin<=kbinmax; bin++) {
1141 Double_t xxxx = fHisto->GetBinCenter(bin);
b90de01a 1142 CalcMult(ipart,r,xxxx,continuous,discrete);
b677de56 1143 fHisto->SetBinContent(bin,continuous);
1144 }
1145 } else {
1146 for(Int_t bin=1; bin<=kbinmax; bin++) {
1147 Double_t xxxx = fHisto->GetBinCenter(bin);
b90de01a 1148 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
b677de56 1149 fHisto->SetBinContent(bin,continuous);
1150 }
1151 }
1152
1153 //non-reweighted P(Delta E = 0)
1154 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1155 Double_t val=discrete*fgkBins/fgkMaxBin;
1156 fHisto->Fill(0.,val);
1157 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1158 Double_t hint=fHisto->Integral(1,kbinmax+1);
1159 p=fHisto->GetBinContent(1)/hint;
1160
1161 // reweighted
1162 hint=fHisto->Integral(1,kbinmax);
1163 prw=fHisto->GetBinContent(1)/hint;
1164
1165 Double_t xxxx = fHisto->GetBinCenter(1);
b90de01a 1166 CalcMult(ipart,r,xxxx,continuous,discrete);
b677de56 1167 fHisto->SetBinContent(1,continuous);
1168 hint=fHisto->Integral(1,kbinmax);
1169 fHisto->Scale(1./hint*(1-discrete));
1170 fHisto->Fill(0.,discrete);
b90de01a 1171 prwcont=fHisto->GetBinContent(1);
b677de56 1172}
1173
b6d061b7 1174Int_t AliQuenchingWeights::SampleEnergyLoss()
1175{
1176 // Has to be called to fill the histograms
1177 //
1178 // For stored values fQTransport loop over
1179 // particle type and length = 1 to fMaxLength (fm)
1180 // to fill energy loss histos
1181 //
1182 // Take histogram of continuous weights
1183 // Take discrete_weight
1184 // If discrete_weight > 1, put all channels to 0, except channel 1
1185 // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral
1186
1187 // read-in data before first call
1188 if(!fTablesLoaded){
7258586f 1189 Error("SampleEnergyLoss","Tables are not loaded.");
b6d061b7 1190 return -1;
1191 }
1192
1193 if(fMultSoft) {
1194 Int_t lmax=CalcLengthMax(fQTransport);
1195 if(fLengthMax>lmax){
7a76a12e 1196 Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
b6d061b7 1197 fLengthMax=lmax;
1198 }
1199 } else {
1200 Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested.");
1201 }
1202
1203 Reset();
1204 fHistos=new TH1F**[2];
9d851d20 1205 fHistos[0]=new TH1F*[4*fLengthMax];
1206 fHistos[1]=new TH1F*[4*fLengthMax];
b6d061b7 1207 fLengthMaxOld=fLengthMax; //remember old value in case
1208 //user wants to reset
1209
1210 Int_t medvalue=0;
1211 Char_t meddesc[100];
1212 if(fMultSoft) {
1213 medvalue=(Int_t)(fQTransport*1000.);
1214 sprintf(meddesc,"MS");
1215 } else {
1216 medvalue=(Int_t)(fMu*1000.);
1217 sprintf(meddesc,"SH");
1218 }
1219
1220 for(Int_t ipart=1;ipart<=2;ipart++){
9d851d20 1221 for(Int_t l=1;l<=4*fLengthMax;l++){
b6d061b7 1222 Char_t hname[100];
1223 sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
9d851d20 1224 Double_t len=l/4.;
7258586f 1225 Double_t wc = CalcWC(len);
2552c51a 1226 fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgkBins,0.,fgkMaxBin*wc);
b6d061b7 1227 fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
1228 fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
1229 fHistos[ipart-1][l-1]->SetLineColor(4);
1230
7258586f 1231 Double_t rrrr = CalcR(wc,len);
b6d061b7 1232 Double_t discrete=0.;
1233 // loop on histogram channels
2552c51a 1234 for(Int_t bin=1; bin<=fgkBins; bin++) {
b6d061b7 1235 Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
1236 Double_t continuous;
7258586f 1237 if(fMultSoft)
1238 CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1239 else
1240 CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
b6d061b7 1241 fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
1242 }
1243 // add discrete part to distribution
1244 if(discrete>=1.)
2552c51a 1245 for(Int_t bin=2;bin<=fgkBins;bin++)
b6d061b7 1246 fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
1247 else {
2552c51a 1248 Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgkBins);
b6d061b7 1249 fHistos[ipart-1][l-1]->Fill(0.,val);
1250 }
2552c51a 1251 Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins);
b6d061b7 1252 fHistos[ipart-1][l-1]->Scale(1./hint);
1253 }
1254 }
1255 return 0;
1256}
1257
b90de01a 1258Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t r)
7258586f 1259{
1260 // Sample energy loss directly for one particle type
1261 // choses R (safe it and keep it until next call of function)
1262
1263 // read-in data before first call
1264 if(!fTablesLoaded){
1265 Error("SampleEnergyLoss","Tables are not loaded.");
1266 return -1;
1267 }
1268
1269 Double_t discrete=0.;
1270 Double_t continuous=0;;
1271 Int_t bin=1;
1272 Double_t xxxx = fHisto->GetBinCenter(bin);
1273 if(fMultSoft)
b90de01a 1274 CalcMult(ipart,r,xxxx,continuous,discrete);
7258586f 1275 else
b90de01a 1276 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
7258586f 1277
1278 if(discrete>=1.) {
1279 fHisto->SetBinContent(1,1.);
2552c51a 1280 for(Int_t bin=2;bin<=fgkBins;bin++)
7258586f 1281 fHisto->SetBinContent(bin,0.);
1282 return 0;
1283 }
1284
1285 fHisto->SetBinContent(bin,continuous);
2552c51a 1286 for(Int_t bin=2; bin<=fgkBins; bin++) {
7258586f 1287 xxxx = fHisto->GetBinCenter(bin);
1288 if(fMultSoft)
b90de01a 1289 CalcMult(ipart,r,xxxx,continuous,discrete);
7258586f 1290 else
b90de01a 1291 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
7258586f 1292 fHisto->SetBinContent(bin,continuous);
1293 }
8ab8044e 1294
2552c51a 1295 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
7258586f 1296 fHisto->Fill(0.,val);
2552c51a 1297 Double_t hint=fHisto->Integral(1,fgkBins);
7258586f 1298 if(hint!=0)
1299 fHisto->Scale(1./hint);
1300 else {
2552c51a 1301 //cout << discrete << " " << hint << " " << continuous << endl;
7258586f 1302 return -1;
1303 }
1304 return 0;
1305}
1306
1307const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) const
b6d061b7 1308{
7a76a12e 1309 //return quenching histograms
1310 //for ipart and length
1311
b6d061b7 1312 if(!fHistos){
1313 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
1314 return 0;
1315 }
1316 if((ipart<1) || (ipart>2)) {
1317 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
1318 return 0;
1319 }
1320
7258586f 1321 Int_t l=GetIndex(length);
b6d061b7 1322 if(l<=0) return 0;
b6d061b7 1323 return fHistos[ipart-1][l-1];
1324}
1325
1326TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const
1327{
1328 // ipart = 1 for quark, 2 for gluon
1329 // medval a) qtransp = transport coefficient (GeV^2/fm)
1330 // b) mu = Debye mass (GeV)
1331 // length = path length in medium (fm)
1332 // Get from SW tables:
7258586f 1333 // - continuous weight, as a function of dE/wc
b6d061b7 1334
1335 Double_t wc = 0;
1336 Char_t meddesc[100];
1337 if(fMultSoft) {
1338 wc=CalcWC(medval,length);
1339 sprintf(meddesc,"MS");
1340 } else {
1341 wc=CalcWCbar(medval,length);
1342 sprintf(meddesc,"SH");
1343 }
1344
1345 Char_t hname[100];
1346 sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart,
1347 (Int_t)(medval*1000.),(Int_t)length);
1348
2552c51a 1349 TH1F *hist = new TH1F("hist",hname,fgkBins,0.,fgkMaxBin*wc);
b6d061b7 1350 hist->SetXTitle("#Delta E [GeV]");
1351 hist->SetYTitle("p(#Delta E)");
1352 hist->SetLineColor(4);
1353
1354 Double_t rrrr = CalcR(wc,length);
7258586f 1355 //loop on histogram channels
2552c51a 1356 for(Int_t bin=1; bin<=fgkBins; bin++) {
b6d061b7 1357 Double_t xxxx = hist->GetBinCenter(bin)/wc;
1358 Double_t continuous,discrete;
1359 Int_t ret=0;
1360 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1361 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1362 if(ret!=0){
1363 delete hist;
1364 return 0;
1365 };
1366 hist->SetBinContent(bin,continuous);
1367 }
1368 return hist;
1369}
1370
1371TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const
1372{
1373 // ipart = 1 for quark, 2 for gluon
1374 // medval a) qtransp = transport coefficient (GeV^2/fm)
1375 // b) mu = Debye mass (GeV)
1376 // length = path length in medium (fm)
1377 // Get from SW tables:
7258586f 1378 // - continuous weight, as a function of dE/wc
b6d061b7 1379
1380 Double_t wc = 0;
1381 Char_t meddesc[100];
1382 if(fMultSoft) {
1383 wc=CalcWC(medval,length);
1384 sprintf(meddesc,"MS");
1385 } else {
1386 wc=CalcWCbar(medval,length);
1387 sprintf(meddesc,"SH");
1388 }
1389
1390 Char_t hname[100];
1391 sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart,
1392 (Int_t)(medval*1000.),(Int_t)length);
1393
2552c51a 1394 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
b6d061b7 1395 histx->SetXTitle("x = #Delta E/#omega_{c}");
1396 if(fMultSoft)
1397 histx->SetYTitle("p(#Delta E/#omega_{c})");
1398 else
1399 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1400 histx->SetLineColor(4);
1401
1402 Double_t rrrr = CalcR(wc,length);
7258586f 1403 //loop on histogram channels
2552c51a 1404 for(Int_t bin=1; bin<=fgkBins; bin++) {
b6d061b7 1405 Double_t xxxx = histx->GetBinCenter(bin);
1406 Double_t continuous,discrete;
1407 Int_t ret=0;
1408 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1409 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1410 if(ret!=0){
1411 delete histx;
1412 return 0;
1413 };
1414 histx->SetBinContent(bin,continuous);
1415 }
1416 return histx;
1417}
1418
b90de01a 1419TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t r) const
7258586f 1420{
1421 // compute P(E) distribution for
1422 // given ipart = 1 for quark, 2 for gluon
1423 // and R
1424
1425 Char_t meddesc[100];
1426 if(fMultSoft) {
1427 sprintf(meddesc,"MS");
1428 } else {
1429 sprintf(meddesc,"SH");
1430 }
1431
1432 Char_t hname[100];
b90de01a 1433 sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,r);
2552c51a 1434 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
7258586f 1435 histx->SetXTitle("x = #Delta E/#omega_{c}");
1436 if(fMultSoft)
1437 histx->SetYTitle("p(#Delta E/#omega_{c})");
1438 else
1439 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1440 histx->SetLineColor(4);
1441
b90de01a 1442 Double_t rrrr = r;
7258586f 1443 Double_t continuous=0.,discrete=0.;
1444 //loop on histogram channels
2552c51a 1445 for(Int_t bin=1; bin<=fgkBins; bin++) {
7258586f 1446 Double_t xxxx = histx->GetBinCenter(bin);
1447 Int_t ret=0;
1448 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1449 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1450 if(ret!=0){
1451 delete histx;
1452 return 0;
1453 };
1454 histx->SetBinContent(bin,continuous);
1455 }
1456
1457 //add discrete part to distribution
1458 if(discrete>=1.)
2552c51a 1459 for(Int_t bin=2;bin<=fgkBins;bin++)
7258586f 1460 histx->SetBinContent(bin,0.);
1461 else {
2552c51a 1462 Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins);
7258586f 1463 histx->Fill(0.,val);
1464 }
2552c51a 1465 Double_t hint=histx->Integral(1,fgkBins);
7258586f 1466 if(hint!=0) histx->Scale(1./hint);
1467
1468 return histx;
1469}
1470
e99e3ed5 1471TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
b6d061b7 1472{
7258586f 1473 // compute energy loss histogram for
1474 // parton type, medium value, length and energy
7a76a12e 1475
b6d061b7 1476 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1477 if(fMultSoft){
1478 dummy->SetQTransport(medval);
1479 dummy->InitMult();
1480 } else {
1481 dummy->SetMu(medval);
1482 dummy->InitSingleHard();
1483 }
1484 dummy->SampleEnergyLoss();
1485
1486 Char_t name[100];
1487 Char_t hname[100];
1488 if(ipart==1){
1489 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1490 sprintf(hname,"hLossQuarks");
1491 } else {
1492 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1493 sprintf(hname,"hLossGluons");
1494 }
1495
1496 TH1F *h = new TH1F(hname,name,250,0,250);
1497 for(Int_t i=0;i<100000;i++){
1498 //if(i % 1000 == 0) cout << "." << flush;
1499 Double_t loss=dummy->GetELossRandom(ipart,l,e);
1500 h->Fill(loss);
1501 }
b6d061b7 1502 h->SetStats(kTRUE);
b6d061b7 1503 delete dummy;
1504 return h;
1505}
1506
e99e3ed5 1507TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
b6d061b7 1508{
7258586f 1509 // compute energy loss histogram for
1510 // parton type, medium value,
1511 // length distribution and energy
7a76a12e 1512
b6d061b7 1513 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1514 if(fMultSoft){
1515 dummy->SetQTransport(medval);
1516 dummy->InitMult();
1517 } else {
1518 dummy->SetMu(medval);
1519 dummy->InitSingleHard();
1520 }
1521 dummy->SampleEnergyLoss();
1522
1523 Char_t name[100];
1524 Char_t hname[100];
1525 if(ipart==1){
1526 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1527 sprintf(hname,"hLossQuarks");
1528 } else {
1529 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1530 sprintf(hname,"hLossGluons");
1531 }
1532
1533 TH1F *h = new TH1F(hname,name,250,0,250);
1534 for(Int_t i=0;i<100000;i++){
1535 //if(i % 1000 == 0) cout << "." << flush;
1536 Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
1537 h->Fill(loss);
1538 }
b6d061b7 1539 h->SetStats(kTRUE);
b6d061b7 1540 delete dummy;
1541 return h;
1542}
1543
b90de01a 1544TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t r) const
7258586f 1545{
1546 // compute energy loss histogram for
1547 // parton type and given R
1548
b90de01a 1549 TH1F *dummy = ComputeQWHistoX(ipart,r);
7258586f 1550 if(!dummy) return 0;
1551
1552 Char_t hname[100];
b90de01a 1553 sprintf(hname,"hELossHistox_%d_%.2f",ipart,r);
2552c51a 1554 TH1F *histx = new TH1F("histxr",hname,fgkBins,0.,fgkMaxBin);
7258586f 1555 for(Int_t i=0;i<100000;i++){
1556 //if(i % 1000 == 0) cout << "." << flush;
1557 Double_t loss=dummy->GetRandom();
1558 histx->Fill(loss);
1559 }
1560 delete dummy;
1561 return histx;
1562}
1563
1564Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,Double_t l) const
1565{
1566 // compute average energy loss for
1567 // parton type, medium value, length and energy
1568
1569 TH1F *dummy = ComputeELossHisto(ipart,medval,l);
1570 if(!dummy) return 0;
1571 Double_t ret=dummy->GetMean();
1572 delete dummy;
1573 return ret;
1574}
1575
1576Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEll) const
1577{
1578 // compute average energy loss for
1579 // parton type, medium value,
1580 // length distribution and energy
1581
1582 TH1F *dummy = ComputeELossHisto(ipart,medval,hEll);
1583 if(!dummy) return 0;
1584 Double_t ret=dummy->GetMean();
1585 delete dummy;
1586 return ret;
1587}
1588
b90de01a 1589Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t r) const
7258586f 1590{
1591 // compute average energy loss over wc
1592 // for parton type and given R
1593
b90de01a 1594 TH1F *dummy = ComputeELossHisto(ipart,r);
7258586f 1595 if(!dummy) return 0;
1596 Double_t ret=dummy->GetMean();
1597 delete dummy;
1598 return ret;
1599}
1600
1601void AliQuenchingWeights::PlotDiscreteWeights(Double_t len) const
b6d061b7 1602{
7258586f 1603 // plot discrete weights for given length
7a76a12e 1604
b6d061b7 1605 TCanvas *c;
1606 if(fMultSoft)
1607 c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
1608 else
1609 c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
1610 c->cd();
1611
1612 TH2F *hframe = new TH2F("hdisc","",2,0,5.1,2,0,1);
1613 hframe->SetStats(0);
1614 if(fMultSoft)
1615 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1616 else
1617 hframe->SetXTitle("#mu [GeV]");
1618 hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
1619 hframe->Draw();
1620
1621 TGraph *gq=new TGraph(20);
1622 Int_t i=0;
1623 if(fMultSoft) {
1624 for(Double_t q=0.05;q<=5.05;q+=0.25){
1625 Double_t disc,cont;
7258586f 1626 CalcMult(1,1.0,q,len,cont,disc);
b6d061b7 1627 gq->SetPoint(i,q,disc);i++;
1628 }
1629 } else {
1630 for(Double_t m=0.05;m<=5.05;m+=0.25){
1631 Double_t disc,cont;
7258586f 1632 CalcSingleHard(1,1.0,m,len,cont, disc);
b6d061b7 1633 gq->SetPoint(i,m,disc);i++;
1634 }
1635 }
1636 gq->SetMarkerStyle(20);
1637 gq->Draw("pl");
1638
1639 TGraph *gg=new TGraph(20);
1640 i=0;
1641 if(fMultSoft){
1642 for(Double_t q=0.05;q<=5.05;q+=0.25){
1643 Double_t disc,cont;
7258586f 1644 CalcMult(2,1.0,q,len,cont,disc);
b6d061b7 1645 gg->SetPoint(i,q,disc);i++;
1646 }
1647 } else {
1648 for(Double_t m=0.05;m<=5.05;m+=0.25){
1649 Double_t disc,cont;
7258586f 1650 CalcSingleHard(2,1.0,m,len,cont,disc);
b6d061b7 1651 gg->SetPoint(i,m,disc);i++;
1652 }
1653 }
1654 gg->SetMarkerStyle(24);
1655 gg->Draw("pl");
1656
1657 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1658 l1a->SetFillStyle(0);
1659 l1a->SetBorderSize(0);
1660 Char_t label[100];
7258586f 1661 sprintf(label,"L = %.1f fm",len);
b6d061b7 1662 l1a->AddEntry(gq,label,"");
1663 l1a->AddEntry(gq,"quark","pl");
1664 l1a->AddEntry(gg,"gluon","pl");
1665 l1a->Draw();
1666
1667 c->Update();
1668}
1669
7258586f 1670void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
b6d061b7 1671{
7258586f 1672 // plot continous weights for
1673 // given parton type and length
7a76a12e 1674
b6d061b7 1675 Float_t medvals[3];
1676 Char_t title[1024];
1677 Char_t name[1024];
1678 if(fMultSoft) {
1679 if(itype==1)
1680 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1681 else if(itype==2)
1682 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1683 else return;
1684 medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1685 sprintf(name,"ccont-ms-%d",itype);
1686 } else {
1687 if(itype==1)
1688 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1689 else if(itype==2)
1690 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1691 else return;
1692 medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1693 sprintf(name,"ccont-ms-%d",itype);
1694 }
1695
1696 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1697 c->cd();
1698 TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1699 h1->SetName("h1");
1700 h1->SetTitle(title);
1701 h1->SetStats(0);
1702 h1->SetLineColor(1);
1703 h1->DrawCopy();
1704 TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1705 h2->SetName("h2");
1706 h2->SetLineColor(2);
1707 h2->DrawCopy("SAME");
1708 TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1709 h3->SetName("h3");
1710 h3->SetLineColor(3);
1711 h3->DrawCopy("SAME");
1712
1713 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1714 l1a->SetFillStyle(0);
1715 l1a->SetBorderSize(0);
1716 Char_t label[100];
7258586f 1717 sprintf(label,"L = %.1f fm",ell);
b6d061b7 1718 l1a->AddEntry(h1,label,"");
1719 if(fMultSoft) {
1720 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1721 l1a->AddEntry(h1,label,"pl");
1722 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1723 l1a->AddEntry(h2,label,"pl");
1724 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1725 l1a->AddEntry(h3,label,"pl");
1726 } else {
1727 sprintf(label,"#mu = %.1f GeV",medvals[0]);
1728 l1a->AddEntry(h1,label,"pl");
1729 sprintf(label,"#mu = %.1f GeV",medvals[1]);
1730 l1a->AddEntry(h2,label,"pl");
1731 sprintf(label,"#mu = %.1f GeV",medvals[2]);
1732 l1a->AddEntry(h3,label,"pl");
1733 }
1734 l1a->Draw();
1735
1736 c->Update();
1737}
1738
7258586f 1739void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
b6d061b7 1740{
7258586f 1741 // plot continous weights for
1742 // given parton type and medium value
7a76a12e 1743
b6d061b7 1744 Char_t title[1024];
1745 Char_t name[1024];
1746 if(fMultSoft) {
1747 if(itype==1)
1748 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1749 else if(itype==2)
1750 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1751 else return;
1752 sprintf(name,"ccont2-ms-%d",itype);
1753 } else {
1754 if(itype==1)
1755 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1756 else if(itype==2)
1757 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1758 else return;
1759 sprintf(name,"ccont2-sh-%d",itype);
1760 }
1761 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1762 c->cd();
1763 TH1F *h1=ComputeQWHisto(itype,medval,8);
1764 h1->SetName("h1");
1765 h1->SetTitle(title);
1766 h1->SetStats(0);
1767 h1->SetLineColor(1);
1768 h1->DrawCopy();
1769 TH1F *h2=ComputeQWHisto(itype,medval,5);
1770 h2->SetName("h2");
1771 h2->SetLineColor(2);
1772 h2->DrawCopy("SAME");
1773 TH1F *h3=ComputeQWHisto(itype,medval,2);
1774 h3->SetName("h3");
1775 h3->SetLineColor(3);
1776 h3->DrawCopy("SAME");
1777
1778 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1779 l1a->SetFillStyle(0);
1780 l1a->SetBorderSize(0);
1781 Char_t label[100];
1782 if(fMultSoft)
1783 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval);
1784 else
1785 sprintf(label,"#mu = %.1f GeV",medval);
1786
1787 l1a->AddEntry(h1,label,"");
1788 l1a->AddEntry(h1,"L = 8 fm","pl");
1789 l1a->AddEntry(h2,"L = 5 fm","pl");
1790 l1a->AddEntry(h3,"L = 2 fm","pl");
1791 l1a->Draw();
1792
1793 c->Update();
1794}
1795
7258586f 1796void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t e) const
b6d061b7 1797{
7258586f 1798 // plot average energy loss for given length
1799 // and parton energy
7a76a12e 1800
b6d061b7 1801 if(!fTablesLoaded){
7258586f 1802 Error("PlotAvgELoss","Tables are not loaded.");
b6d061b7 1803 return;
1804 }
1805
1806 Char_t title[1024];
1807 Char_t name[1024];
1808 if(fMultSoft){
7258586f 1809 sprintf(title,"Average Energy Loss for Multiple Scattering");
b6d061b7 1810 sprintf(name,"cavgelossms");
1811 } else {
7258586f 1812 sprintf(title,"Average Energy Loss for Single Hard Scattering");
b6d061b7 1813 sprintf(name,"cavgelosssh");
1814 }
1815
1816 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1817 c->cd();
1818 TH2F *hframe = new TH2F("avgloss",title,2,0,5.1,2,0,100);
1819 hframe->SetStats(0);
1820 if(fMultSoft)
1821 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1822 else
1823 hframe->SetXTitle("#mu [GeV]");
1824 hframe->SetYTitle("<E_{loss}> [GeV]");
1825 hframe->Draw();
1826
1827 TGraph *gq=new TGraph(20);
1828 Int_t i=0;
1829 for(Double_t v=0.05;v<=5.05;v+=0.25){
1830 TH1F *dummy=ComputeELossHisto(1,v,len,e);
1831 Double_t avgloss=dummy->GetMean();
1832 gq->SetPoint(i,v,avgloss);i++;
1833 delete dummy;
1834 }
1835 gq->SetMarkerStyle(20);
1836 gq->Draw("pl");
1837
1838 TGraph *gg=new TGraph(20);
1839 i=0;
1840 for(Double_t v=0.05;v<=5.05;v+=0.25){
1841 TH1F *dummy=ComputeELossHisto(2,v,len,e);
1842 Double_t avgloss=dummy->GetMean();
1843 gg->SetPoint(i,v,avgloss);i++;
1844 delete dummy;
1845 }
1846 gg->SetMarkerStyle(24);
1847 gg->Draw("pl");
1848
7258586f 1849 TGraph *gratio=new TGraph(20);
1850 for(Int_t i=0;i<20;i++){
1851 Double_t x,y,x2,y2;
1852 gg->GetPoint(i,x,y);
1853 gq->GetPoint(i,x2,y2);
1854 if(y2>0)
1855 gratio->SetPoint(i,x,y/y2*10/2.25);
1856 else gratio->SetPoint(i,x,0);
1857 }
1858 gratio->SetLineStyle(4);
1859 gratio->Draw();
b6d061b7 1860 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1861 l1a->SetFillStyle(0);
1862 l1a->SetBorderSize(0);
1863 Char_t label[100];
7258586f 1864 sprintf(label,"L = %.1f fm",len);
b6d061b7 1865 l1a->AddEntry(gq,label,"");
1866 l1a->AddEntry(gq,"quark","pl");
1867 l1a->AddEntry(gg,"gluon","pl");
7258586f 1868 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
b6d061b7 1869 l1a->Draw();
1870
1871 c->Update();
1872}
1873
e99e3ed5 1874void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
b6d061b7 1875{
7258586f 1876 // plot average energy loss for given
1877 // length distribution and parton energy
7a76a12e 1878
b6d061b7 1879 if(!fTablesLoaded){
7258586f 1880 Error("PlotAvgELossVs","Tables are not loaded.");
b6d061b7 1881 return;
1882 }
1883
1884 Char_t title[1024];
1885 Char_t name[1024];
1886 if(fMultSoft){
7258586f 1887 sprintf(title,"Average Energy Loss for Multiple Scattering");
b6d061b7 1888 sprintf(name,"cavgelossms2");
1889 } else {
7258586f 1890 sprintf(title,"Average Energy Loss for Single Hard Scattering");
b6d061b7 1891 sprintf(name,"cavgelosssh2");
1892 }
1893
1894 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1895 c->cd();
1896 TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1897 hframe->SetStats(0);
1898 if(fMultSoft)
1899 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1900 else
1901 hframe->SetXTitle("#mu [GeV]");
1902 hframe->SetYTitle("<E_{loss}> [GeV]");
1903 hframe->Draw();
1904
1905 TGraph *gq=new TGraph(20);
1906 Int_t i=0;
1907 for(Double_t v=0.05;v<=5.05;v+=0.25){
1908 TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1909 Double_t avgloss=dummy->GetMean();
1910 gq->SetPoint(i,v,avgloss);i++;
1911 delete dummy;
1912 }
1913 gq->SetMarkerStyle(20);
1914 gq->Draw("pl");
1915
1916 TGraph *gg=new TGraph(20);
1917 i=0;
1918 for(Double_t v=0.05;v<=5.05;v+=0.25){
1919 TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1920 Double_t avgloss=dummy->GetMean();
1921 gg->SetPoint(i,v,avgloss);i++;
1922 delete dummy;
1923 }
1924 gg->SetMarkerStyle(24);
1925 gg->Draw("pl");
1926
7258586f 1927 TGraph *gratio=new TGraph(20);
1928 for(Int_t i=0;i<20;i++){
1929 Double_t x,y,x2,y2;
1930 gg->GetPoint(i,x,y);
1931 gq->GetPoint(i,x2,y2);
1932 if(y2>0)
1933 gratio->SetPoint(i,x,y/y2*10/2.25);
1934 else gratio->SetPoint(i,x,0);
1935 }
1936 gratio->SetLineStyle(4);
1937 gratio->Draw();
1938
b6d061b7 1939 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1940 l1a->SetFillStyle(0);
1941 l1a->SetBorderSize(0);
1942 Char_t label[100];
1943 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1944 l1a->AddEntry(gq,label,"");
1945 l1a->AddEntry(gq,"quark","pl");
1946 l1a->AddEntry(gg,"gluon","pl");
7258586f 1947 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
b6d061b7 1948 l1a->Draw();
1949
1950 c->Update();
1951}
1952
9d851d20 1953void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const
b6d061b7 1954{
7258586f 1955 // plot average energy loss versus ell
7a76a12e 1956
b6d061b7 1957 if(!fTablesLoaded){
7258586f 1958 Error("PlotAvgELossVsEll","Tables are not loaded.");
1959 return;
1960 }
1961
1962 Char_t title[1024];
1963 Char_t name[1024];
1964 Float_t medval;
1965 if(fMultSoft){
1966 sprintf(title,"Average Energy Loss for Multiple Scattering");
1967 sprintf(name,"cavgelosslms");
1968 medval=fQTransport;
1969 } else {
1970 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1971 sprintf(name,"cavgelosslsh");
1972 medval=fMu;
1973 }
1974
1975 TCanvas *c = new TCanvas(name,title,0,0,600,400);
1976 c->cd();
1977 TH2F *hframe = new TH2F("avglossell",title,2,0,fLengthMax,2,0,250);
1978 hframe->SetStats(0);
1979 hframe->SetXTitle("length [fm]");
1980 hframe->SetYTitle("<E_{loss}> [GeV]");
1981 hframe->Draw();
1982
1983 TGraph *gq=new TGraph((Int_t)fLengthMax*4);
1984 Int_t i=0;
1985 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
1986 TH1F *dummy=ComputeELossHisto(1,medval,v,e);
1987 Double_t avgloss=dummy->GetMean();
1988 gq->SetPoint(i,v,avgloss);i++;
1989 delete dummy;
1990 }
1991 gq->SetMarkerStyle(20);
1992 gq->Draw("pl");
1993
1994 TGraph *gg=new TGraph((Int_t)fLengthMax*4);
1995 i=0;
1996 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
1997 TH1F *dummy=ComputeELossHisto(2,medval,v,e);
1998 Double_t avgloss=dummy->GetMean();
1999 gg->SetPoint(i,v,avgloss);i++;
2000 delete dummy;
2001 }
2002 gg->SetMarkerStyle(24);
2003 gg->Draw("pl");
2004
2005 TGraph *gratio=new TGraph((Int_t)fLengthMax*4);
2006 for(Int_t i=0;i<=(Int_t)fLengthMax*4;i++){
2007 Double_t x,y,x2,y2;
2008 gg->GetPoint(i,x,y);
2009 gq->GetPoint(i,x2,y2);
2010 if(y2>0)
2011 gratio->SetPoint(i,x,y/y2*100/2.25);
2012 else gratio->SetPoint(i,x,0);
2013 }
2014 gratio->SetLineStyle(1);
2015 gratio->SetLineWidth(2);
2016 gratio->Draw();
2017 TLegend *l1a = new TLegend(0.15,0.65,.60,0.85);
2018 l1a->SetFillStyle(0);
2019 l1a->SetBorderSize(0);
2020 Char_t label[100];
2021 if(fMultSoft)
2022 sprintf(label,"#hat{q} = %.2f GeV^{2}/fm",medval);
2023 else
2024 sprintf(label,"#mu = %.2f GeV",medval);
2025 l1a->AddEntry(gq,label,"");
2026 l1a->AddEntry(gq,"quark","pl");
2027 l1a->AddEntry(gg,"gluon","pl");
2028 l1a->AddEntry(gratio,"gluon/quark/2.25*100","pl");
2029 l1a->Draw();
2030
2031 TF1 *f=new TF1("f","100",0,fLengthMax);
2032 f->SetLineStyle(4);
2033 f->SetLineWidth(1);
2034 f->Draw("same");
2035 c->Update();
2036}
2037
2038void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
2039{
2040 // plot relative energy loss for given
2041 // length and parton energy versus pt
2042
2043 if(!fTablesLoaded){
2044 Error("PlotAvgELossVsPt","Tables are not loaded.");
b6d061b7 2045 return;
2046 }
2047
2048 Char_t title[1024];
2049 Char_t name[1024];
2050 if(fMultSoft){
7258586f 2051 sprintf(title,"Relative Energy Loss for Multiple Scattering");
b6d061b7 2052 sprintf(name,"cavgelossvsptms");
2053 } else {
7258586f 2054 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
b6d061b7 2055 sprintf(name,"cavgelossvsptsh");
2056 }
2057
2058 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2059 c->cd();
2060 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2061 hframe->SetStats(0);
2062 hframe->SetXTitle("p_{T} [GeV]");
2063 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2064 hframe->Draw();
2065
2066 TGraph *gq=new TGraph(40);
2067 Int_t i=0;
2068 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2069 TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
2070 Double_t avgloss=dummy->GetMean();
2071 gq->SetPoint(i,pt,avgloss/pt);i++;
2072 delete dummy;
2073 }
2074 gq->SetMarkerStyle(20);
2075 gq->Draw("pl");
2076
2077 TGraph *gg=new TGraph(40);
2078 i=0;
2079 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2080 TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
2081 Double_t avgloss=dummy->GetMean();
2082 gg->SetPoint(i,pt,avgloss/pt);i++;
2083 delete dummy;
2084 }
2085 gg->SetMarkerStyle(24);
2086 gg->Draw("pl");
2087
2088 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2089 l1a->SetFillStyle(0);
2090 l1a->SetBorderSize(0);
2091 Char_t label[100];
7258586f 2092 sprintf(label,"L = %.1f fm",len);
b6d061b7 2093 l1a->AddEntry(gq,label,"");
2094 l1a->AddEntry(gq,"quark","pl");
2095 l1a->AddEntry(gg,"gluon","pl");
2096 l1a->Draw();
2097
2098 c->Update();
2099}
2100
2101void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
2102{
7258586f 2103 // plot relative energy loss for given
2104 // length distribution and parton energy versus pt
7a76a12e 2105
b6d061b7 2106 if(!fTablesLoaded){
7258586f 2107 Error("PlotAvgELossVsPt","Tables are not loaded.");
b6d061b7 2108 return;
2109 }
2110
2111 Char_t title[1024];
2112 Char_t name[1024];
2113 if(fMultSoft){
7258586f 2114 sprintf(title,"Relative Energy Loss for Multiple Scattering");
b6d061b7 2115 sprintf(name,"cavgelossvsptms2");
2116 } else {
7258586f 2117 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
b6d061b7 2118 sprintf(name,"cavgelossvsptsh2");
2119 }
2120 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2121 c->cd();
2122 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2123 hframe->SetStats(0);
2124 hframe->SetXTitle("p_{T} [GeV]");
2125 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2126 hframe->Draw();
2127
2128 TGraph *gq=new TGraph(40);
2129 Int_t i=0;
2130 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2131 TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
2132 Double_t avgloss=dummy->GetMean();
2133 gq->SetPoint(i,pt,avgloss/pt);i++;
2134 delete dummy;
2135 }
2136 gq->SetMarkerStyle(20);
2137 gq->Draw("pl");
2138
2139 TGraph *gg=new TGraph(40);
2140 i=0;
2141 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2142 TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
2143 Double_t avgloss=dummy->GetMean();
2144 gg->SetPoint(i,pt,avgloss/pt);i++;
2145 delete dummy;
2146 }
2147 gg->SetMarkerStyle(24);
2148 gg->Draw("pl");
2149
2150 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2151 l1a->SetFillStyle(0);
2152 l1a->SetBorderSize(0);
2153 Char_t label[100];
2154 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
2155 l1a->AddEntry(gq,label,"");
2156 l1a->AddEntry(gq,"quark","pl");
2157 l1a->AddEntry(gg,"gluon","pl");
2158 l1a->Draw();
2159
2160 c->Update();
2161}
2162
7258586f 2163Int_t AliQuenchingWeights::GetIndex(Double_t len) const
2164{
b90de01a 2165 //get the index according to length
7258586f 2166 if(len>fLengthMax) len=fLengthMax;
2167
9d851d20 2168 Int_t l=Int_t(len/0.25);
cc885e36 2169 if((len-l*0.25)>0.125) l++;
7258586f 2170 return l;
2171}
ea16e52f 2172