]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - FASTSIM/AliQuenchingWeights.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Avoid warning messages.
[u/mrichter/AliRoot.git] / FASTSIM / AliQuenchingWeights.cxx
... / ...
CommitLineData
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//
26//
27// Origin: C. Loizides constantinos.loizides@cern.ch
28// A. Dainese andrea.dainese@pd.infn.it
29//
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//-----------------------------------------------------------------------------
35
36#include <Riostream.h>
37#include <TF1.h>
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
50const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197;
51
52// maximum value of R
53const Double_t AliQuenchingWeights::fgkRMax = 1.e6;
54
55// hist binning
56const Int_t AliQuenchingWeights::fgkBins = 1300;
57const Double_t AliQuenchingWeights::fgkMaxBin = 1.3;
58
59// counter for histogram labels
60Int_t AliQuenchingWeights::fgCounter = 0;
61
62
63AliQuenchingWeights::AliQuenchingWeights()
64 : TObject()
65{
66 //default constructor
67
68 fTablesLoaded=kFALSE;
69 fMultSoft=kTRUE;
70 fHistos=0;
71 SetMu();
72 SetQTransport();
73 SetK();
74 fECMethod=kDefault;
75 SetLengthMax();
76 fLengthMaxOld=0;
77 fInstanceNumber=fgCounter++;
78 Char_t name[100];
79 sprintf(name,"hhistoqw_%d",fInstanceNumber);
80 fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin);
81 for(Int_t bin=1;bin<=fgkBins;bin++)
82 fHisto->SetBinContent(bin,0.);
83}
84
85AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a)
86 : TObject()
87{
88 // copy constructor
89
90 fTablesLoaded=kFALSE;
91 fHistos=0;
92 fLengthMaxOld=0;
93 fMultSoft=a.GetMultSoft();;
94 fMu=a.GetMu();
95 fK=a.GetK();
96 fQTransport=a.GetQTransport();
97 fECMethod=(kECMethod)a.GetECMethod();
98 fLengthMax=a.GetLengthMax();
99 fInstanceNumber=fgCounter++;
100 Char_t name[100];
101 sprintf(name,"hhistoqw_%d",fInstanceNumber);
102 fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin);
103 for(Int_t bin=1;bin<=fgkBins;bin++)
104 fHisto->SetBinContent(bin,0.);
105
106 //Missing in the class is the pathname
107 //to the tables, can be added if needed
108}
109
110AliQuenchingWeights::~AliQuenchingWeights()
111{
112 Reset();
113 delete fHisto;
114}
115
116void AliQuenchingWeights::Reset()
117{
118 //reset tables if there were used
119
120 if(!fHistos) return;
121 for(Int_t l=0;l<4*fLengthMaxOld;l++){
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{
132 //set energy constraint method
133
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.");
137 else if(fECMethod==kReweight)
138 Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy.");
139 else Info("SetECMethod","Energy Constraint Method set to REWEIGHTCONT:\nRequire sampled energy loss <= parton energy (only implemented for FAST method.");
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));
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
157 if(!fincont.is_open()) return -1;
158#endif
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]>>
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])
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]>>
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 {
182 nn++;
183 if(nn==261) break;
184 }
185 fincont.close();
186
187 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
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
193 if(!findisc.is_open()) return -1;
194#endif
195
196 nn=0; //quarks
197 while(findisc>>frrr[nn]>>fdaq[nn]) {
198 nn++;
199 if(nn==34) break;
200 }
201 nn=0; //gluons
202 while(findisc>>frrrg[nn]>>fdag[nn]) {
203 nn++;
204 if(nn==34) break;
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
241C urs.wiedemann@cern.ch
242C carlos.salgado@cern.ch
243C
244C
245C-------------------------------------------------------------------
246
247 SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete)
248*
249 REAL*8 xx(400), daq(34), caq(34,261), rrr(34)
250 COMMON /dataqua/ xx, daq, caq, rrr
251*
252 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
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*
262
263 continuous=0.d0
264 discrete=0.d0
265
266 rrin = rrrr
267 xxin = xxxx
268*
269 do 666, nr=1,34
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)
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
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
320245 continue
321
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
331 REAL*8 xxq(400), daq(34), caq(34,261), rrr(34)
332 COMMON /dataqua/ xxq, daq, caq, rrr
333*
334 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
335 COMMON /dataglu/ xxg, dag, cag, rrrg
336*
337 OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90)
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),
349 + caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn),
350 + caq(33,nn), caq(34,nn)
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),
363 + cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn),
364 + cag(33,nn), cag(34,nn)
365 111 continue
366 close(20)
367*
368 OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91)
369 do 112 nn=1,34
370 read (21,*) rrr(nn), daq(nn)
371 112 continue
372 do 113 nn=1,34
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
384
385=======================================================================
386
387 Adapted to ROOT macro by A. Dainese - 13/07/2003
388 Ported to class by C. Loizides - 12/02/2004
389 New version for extended R values added - 06/03/2004
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
399 //set result to zero
400 continuous=0.;
401 discrete=0.;
402
403 //read-in data before first call
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
416 if(xxin>fxx[260]) return -1;
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?
423 if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD
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;
428 for(Int_t nr=1; nr<=34; nr++) {
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
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
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",
477 //rfraclow,clow,rfrachigh,chigh,continuous);
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));
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
503 if(!fincont.is_open()) return -1;
504#endif
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));
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
546 if(!findisc.is_open()) return -1;
547#endif
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){
726 Error("CalcSingleHard","Tables are not loaded.");
727 return -1;
728 }
729 if(!fMultSoft){
730 Error("CalcSingleHard","Tables are not loaded for Single Hard Scattering.");
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{
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
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{
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
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{
819 //calculate R value and
820 //check if it is less then maximum
821
822 Double_t R = wc*l*fgkConvFmToInvGeV;
823 if(R>=fgkRMax) {
824 Warning("CalcR","Value of R = %.2f; should be less than %.2f",R,fgkRMax);
825 return fgkRMax-1;
826 }
827 return R;
828}
829
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
835 Double_t R = fgkRMax-1;
836 if(I0>0)
837 R = 2*I1*I1/I0*k;
838 if(R>=fgkRMax) {
839 Warning("CalcRk","Value of R = %.2f; should be less than %.2f",R,fgkRMax);
840 return fgkRMax-1;
841 }
842 return R;
843}
844
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);
858 return -1000.;
859 }
860
861 Int_t l=GetIndex(length);
862 if(l<=0) return 0.;
863
864 if(fECMethod==kReweight){
865 Double_t ret = 2.*e;
866 Int_t ws=0;
867 while(ret>e){
868 ret=fHistos[ipart-1][l-1]->GetRandom();
869 if(++ws==1e6){
870 Warning("GetELossRandom",
871 "Aborted reweighting; maximum loss assigned after 1e6 trials.");
872 return e;
873 }
874 }
875 return ret;
876 }
877 //kDefault
878 Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
879 if(ret>e) return e;
880 return ret;
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
892Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const
893{
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
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
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);
934 if(r<0.){
935 Fatal("GetELossRandomK","R should not be negative");
936 return -1000.;
937 }
938 Double_t wc=CalcWCk(I1);
939 if(wc<=0.){
940 Fatal("GetELossRandomK","wc should be greater than zero");
941 return -1000.;
942 }
943 if(SampleEnergyLoss(ipart,r)!=0){
944 Fatal("GetELossRandomK","Could not sample energy loss");
945 return -1000.;
946 }
947
948 if(fECMethod==kReweight){
949 Double_t ret = 2.*e;
950 Int_t ws=0;
951 while(ret>e){
952 ret=fHisto->GetRandom();
953 if(++ws==1e6){
954 Warning("GetELossRandomK",
955 "Aborted reweighting; maximum loss assigned after 1e6 trials.");
956 return e;
957 }
958 }
959 return ret;
960 }
961
962 //kDefault
963 Double_t ret=fHisto->GetRandom()*wc;
964 if(ret>e) return e;
965 return ret;
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
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);
988 if(r<=0.){
989 return 0.;
990 }
991
992 Double_t wc=CalcWCk(I1);
993 if(wc<=0.){
994 return 0.;
995 }
996
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 }
1013
1014 Double_t discrete=0.;
1015 Double_t continuous=0.;
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
1023 if(discrete>=1.0) {
1024 return 0.; //no energy loss
1025 }
1026
1027 fHisto->SetBinContent(bin,continuous);
1028 Int_t kbinmax=fHisto->FindBin(e/wc);
1029 if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
1030 if(kbinmax==1) return e; //maximum energy loss
1031
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 }
1045
1046 if(fECMethod==kReweight){
1047 fHisto->SetBinContent(kbinmax+1,0);
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);
1058 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1059 }
1060 for(Int_t bin=kbinmax+2; bin<=fgkBins; bin++) {
1061 fHisto->SetBinContent(bin,0);
1062 }
1063 //cout << kbinmax << " " << discrete << " " << fHisto->Integral() << endl;
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
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 &prwcont,
1109 Int_t ipart,Double_t I0,Double_t I1,Double_t e)
1110{
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.){
1116 return;
1117 }
1118 Double_t wc=CalcWCk(I1);
1119 if(wc<=0.){
1120 return;
1121 }
1122 GetZeroLossProbR(p,prw,prwcont,ipart,r,wc,e);
1123}
1124
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)
1127{
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;
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);
1142 CalcMult(ipart,r,xxxx,continuous,discrete);
1143 fHisto->SetBinContent(bin,continuous);
1144 }
1145 } else {
1146 for(Int_t bin=1; bin<=kbinmax; bin++) {
1147 Double_t xxxx = fHisto->GetBinCenter(bin);
1148 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
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);
1166 CalcMult(ipart,r,xxxx,continuous,discrete);
1167 fHisto->SetBinContent(1,continuous);
1168 hint=fHisto->Integral(1,kbinmax);
1169 fHisto->Scale(1./hint*(1-discrete));
1170 fHisto->Fill(0.,discrete);
1171 prwcont=fHisto->GetBinContent(1);
1172}
1173
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){
1189 Error("SampleEnergyLoss","Tables are not loaded.");
1190 return -1;
1191 }
1192
1193 if(fMultSoft) {
1194 Int_t lmax=CalcLengthMax(fQTransport);
1195 if(fLengthMax>lmax){
1196 Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
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];
1205 fHistos[0]=new TH1F*[4*fLengthMax];
1206 fHistos[1]=new TH1F*[4*fLengthMax];
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++){
1221 for(Int_t l=1;l<=4*fLengthMax;l++){
1222 Char_t hname[100];
1223 sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
1224 Double_t len=l/4.;
1225 Double_t wc = CalcWC(len);
1226 fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgkBins,0.,fgkMaxBin*wc);
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
1231 Double_t rrrr = CalcR(wc,len);
1232 Double_t discrete=0.;
1233 // loop on histogram channels
1234 for(Int_t bin=1; bin<=fgkBins; bin++) {
1235 Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
1236 Double_t continuous;
1237 if(fMultSoft)
1238 CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1239 else
1240 CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1241 fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
1242 }
1243 // add discrete part to distribution
1244 if(discrete>=1.)
1245 for(Int_t bin=2;bin<=fgkBins;bin++)
1246 fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
1247 else {
1248 Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1249 fHistos[ipart-1][l-1]->Fill(0.,val);
1250 }
1251 Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1252 fHistos[ipart-1][l-1]->Scale(1./hint);
1253 }
1254 }
1255 return 0;
1256}
1257
1258Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t r)
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)
1274 CalcMult(ipart,r,xxxx,continuous,discrete);
1275 else
1276 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1277
1278 if(discrete>=1.) {
1279 fHisto->SetBinContent(1,1.);
1280 for(Int_t bin=2;bin<=fgkBins;bin++)
1281 fHisto->SetBinContent(bin,0.);
1282 return 0;
1283 }
1284
1285 fHisto->SetBinContent(bin,continuous);
1286 for(Int_t bin=2; bin<=fgkBins; bin++) {
1287 xxxx = fHisto->GetBinCenter(bin);
1288 if(fMultSoft)
1289 CalcMult(ipart,r,xxxx,continuous,discrete);
1290 else
1291 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1292 fHisto->SetBinContent(bin,continuous);
1293 }
1294
1295 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
1296 fHisto->Fill(0.,val);
1297 Double_t hint=fHisto->Integral(1,fgkBins);
1298 if(hint!=0)
1299 fHisto->Scale(1./hint);
1300 else {
1301 //cout << discrete << " " << hint << " " << continuous << endl;
1302 return -1;
1303 }
1304 return 0;
1305}
1306
1307const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) const
1308{
1309 //return quenching histograms
1310 //for ipart and length
1311
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
1321 Int_t l=GetIndex(length);
1322 if(l<=0) return 0;
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:
1333 // - continuous weight, as a function of dE/wc
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
1349 TH1F *hist = new TH1F("hist",hname,fgkBins,0.,fgkMaxBin*wc);
1350 hist->SetXTitle("#Delta E [GeV]");
1351 hist->SetYTitle("p(#Delta E)");
1352 hist->SetLineColor(4);
1353
1354 Double_t rrrr = CalcR(wc,length);
1355 //loop on histogram channels
1356 for(Int_t bin=1; bin<=fgkBins; bin++) {
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:
1378 // - continuous weight, as a function of dE/wc
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
1394 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
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);
1403 //loop on histogram channels
1404 for(Int_t bin=1; bin<=fgkBins; bin++) {
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
1419TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t r) const
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];
1433 sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,r);
1434 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
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
1442 Double_t rrrr = r;
1443 Double_t continuous=0.,discrete=0.;
1444 //loop on histogram channels
1445 for(Int_t bin=1; bin<=fgkBins; bin++) {
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.)
1459 for(Int_t bin=2;bin<=fgkBins;bin++)
1460 histx->SetBinContent(bin,0.);
1461 else {
1462 Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins);
1463 histx->Fill(0.,val);
1464 }
1465 Double_t hint=histx->Integral(1,fgkBins);
1466 if(hint!=0) histx->Scale(1./hint);
1467
1468 return histx;
1469}
1470
1471TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
1472{
1473 // compute energy loss histogram for
1474 // parton type, medium value, length and energy
1475
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 }
1502 h->SetStats(kTRUE);
1503 delete dummy;
1504 return h;
1505}
1506
1507TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
1508{
1509 // compute energy loss histogram for
1510 // parton type, medium value,
1511 // length distribution and energy
1512
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 }
1539 h->SetStats(kTRUE);
1540 delete dummy;
1541 return h;
1542}
1543
1544TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t r) const
1545{
1546 // compute energy loss histogram for
1547 // parton type and given R
1548
1549 TH1F *dummy = ComputeQWHistoX(ipart,r);
1550 if(!dummy) return 0;
1551
1552 Char_t hname[100];
1553 sprintf(hname,"hELossHistox_%d_%.2f",ipart,r);
1554 TH1F *histx = new TH1F("histxr",hname,fgkBins,0.,fgkMaxBin);
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
1589Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t r) const
1590{
1591 // compute average energy loss over wc
1592 // for parton type and given R
1593
1594 TH1F *dummy = ComputeELossHisto(ipart,r);
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,Double_t qm) const
1602{
1603 // plot discrete weights for given length
1604
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,qm+.1,2,0,1.25);
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->SetYTitle("p_{0} (discrete weight)");
1620 hframe->Draw();
1621
1622 Int_t points=(Int_t)qm*4;
1623 TGraph *gq=new TGraph(points);
1624 Int_t i=0;
1625 if(fMultSoft) {
1626 for(Double_t q=0.05;q<=qm+.05;q+=0.25){
1627 Double_t disc,cont;
1628 CalcMult(1,1.0,q,len,cont,disc);
1629 gq->SetPoint(i,q,disc);i++;
1630 }
1631 } else {
1632 for(Double_t m=0.05;m<=qm+.05;m+=0.25){
1633 Double_t disc,cont;
1634 CalcSingleHard(1,1.0,m,len,cont, disc);
1635 gq->SetPoint(i,m,disc);i++;
1636 }
1637 }
1638 gq->SetMarkerStyle(20);
1639 gq->SetMarkerColor(1);
1640 gq->SetLineStyle(1);
1641 gq->SetLineColor(1);
1642 gq->Draw("l");
1643
1644 TGraph *gg=new TGraph(points);
1645 i=0;
1646 if(fMultSoft){
1647 for(Double_t q=0.05;q<=qm+.05;q+=0.25){
1648 Double_t disc,cont;
1649 CalcMult(2,1.0,q,len,cont,disc);
1650 gg->SetPoint(i,q,disc);i++;
1651 }
1652 } else {
1653 for(Double_t m=0.05;m<=qm+.05;m+=0.25){
1654 Double_t disc,cont;
1655 CalcSingleHard(2,1.0,m,len,cont,disc);
1656 gg->SetPoint(i,m,disc);i++;
1657 }
1658 }
1659 gg->SetMarkerStyle(24);
1660 gg->SetMarkerColor(2);
1661 gg->SetLineStyle(2);
1662 gg->SetLineColor(2);
1663 gg->Draw("l");
1664
1665 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1666 l1a->SetFillStyle(0);
1667 l1a->SetBorderSize(0);
1668 Char_t label[100];
1669 sprintf(label,"L = %.1f fm",len);
1670 l1a->AddEntry(gq,label,"");
1671 l1a->AddEntry(gq,"quark projectile","l");
1672 l1a->AddEntry(gg,"gluon projectile","l");
1673 l1a->Draw();
1674
1675 c->Update();
1676}
1677
1678void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
1679{
1680 // plot continous weights for
1681 // given parton type and length
1682
1683 Float_t medvals[3];
1684 Char_t title[1024];
1685 Char_t name[1024];
1686 if(fMultSoft) {
1687 if(itype==1)
1688 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1689 else if(itype==2)
1690 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1691 else return;
1692 medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1693 sprintf(name,"ccont-ms-%d",itype);
1694 } else {
1695 if(itype==1)
1696 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1697 else if(itype==2)
1698 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1699 else return;
1700 medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1701 sprintf(name,"ccont-ms-%d",itype);
1702 }
1703
1704 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1705 c->cd();
1706 TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1707 h1->SetName("h1");
1708 h1->SetTitle(title);
1709 h1->SetStats(0);
1710 h1->SetLineColor(1);
1711 h1->DrawCopy();
1712 TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1713 h2->SetName("h2");
1714 h2->SetLineColor(2);
1715 h2->DrawCopy("SAME");
1716 TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1717 h3->SetName("h3");
1718 h3->SetLineColor(3);
1719 h3->DrawCopy("SAME");
1720
1721 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1722 l1a->SetFillStyle(0);
1723 l1a->SetBorderSize(0);
1724 Char_t label[100];
1725 sprintf(label,"L = %.1f fm",ell);
1726 l1a->AddEntry(h1,label,"");
1727 if(fMultSoft) {
1728 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1729 l1a->AddEntry(h1,label,"pl");
1730 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1731 l1a->AddEntry(h2,label,"pl");
1732 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1733 l1a->AddEntry(h3,label,"pl");
1734 } else {
1735 sprintf(label,"#mu = %.1f GeV",medvals[0]);
1736 l1a->AddEntry(h1,label,"pl");
1737 sprintf(label,"#mu = %.1f GeV",medvals[1]);
1738 l1a->AddEntry(h2,label,"pl");
1739 sprintf(label,"#mu = %.1f GeV",medvals[2]);
1740 l1a->AddEntry(h3,label,"pl");
1741 }
1742 l1a->Draw();
1743
1744 c->Update();
1745}
1746
1747void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
1748{
1749 // plot continous weights for
1750 // given parton type and medium value
1751
1752 Char_t title[1024];
1753 Char_t name[1024];
1754 if(fMultSoft) {
1755 if(itype==1)
1756 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1757 else if(itype==2)
1758 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1759 else return;
1760 sprintf(name,"ccont2-ms-%d",itype);
1761 } else {
1762 if(itype==1)
1763 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1764 else if(itype==2)
1765 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1766 else return;
1767 sprintf(name,"ccont2-sh-%d",itype);
1768 }
1769 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1770 c->cd();
1771 TH1F *h1=ComputeQWHisto(itype,medval,8);
1772 h1->SetName("h1");
1773 h1->SetTitle(title);
1774 h1->SetStats(0);
1775 h1->SetLineColor(1);
1776 h1->DrawCopy();
1777 TH1F *h2=ComputeQWHisto(itype,medval,5);
1778 h2->SetName("h2");
1779 h2->SetLineColor(2);
1780 h2->DrawCopy("SAME");
1781 TH1F *h3=ComputeQWHisto(itype,medval,2);
1782 h3->SetName("h3");
1783 h3->SetLineColor(3);
1784 h3->DrawCopy("SAME");
1785
1786 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1787 l1a->SetFillStyle(0);
1788 l1a->SetBorderSize(0);
1789 Char_t label[100];
1790 if(fMultSoft)
1791 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval);
1792 else
1793 sprintf(label,"#mu = %.1f GeV",medval);
1794
1795 l1a->AddEntry(h1,label,"");
1796 l1a->AddEntry(h1,"L = 8 fm","pl");
1797 l1a->AddEntry(h2,"L = 5 fm","pl");
1798 l1a->AddEntry(h3,"L = 2 fm","pl");
1799 l1a->Draw();
1800
1801 c->Update();
1802}
1803
1804void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t qm,Double_t e) const
1805{
1806 // plot average energy loss for given length
1807 // and parton energy
1808
1809 if(!fTablesLoaded){
1810 Error("PlotAvgELoss","Tables are not loaded.");
1811 return;
1812 }
1813
1814 Char_t title[1024];
1815 Char_t name[1024];
1816 if(fMultSoft){
1817 sprintf(title,"Average Energy Loss for Multiple Scattering");
1818 sprintf(name,"cavgelossms");
1819 } else {
1820 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1821 sprintf(name,"cavgelosssh");
1822 }
1823
1824 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1825 c->cd();
1826 TH2F *hframe = new TH2F("avgloss","",2,0,qm+.1,2,0,100);
1827 hframe->SetStats(0);
1828 if(fMultSoft)
1829 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1830 else
1831 hframe->SetXTitle("#mu [GeV]");
1832 hframe->SetYTitle("<E_{loss}> [GeV]");
1833 hframe->Draw();
1834
1835 TGraph *gq=new TGraph(20);
1836 Int_t i=0;
1837 for(Double_t v=0.05;v<=qm+.05;v+=0.25){
1838 TH1F *dummy=ComputeELossHisto(1,v,len,e);
1839 Double_t avgloss=dummy->GetMean();
1840 gq->SetPoint(i,v,avgloss);i++;
1841 delete dummy;
1842 }
1843 gq->SetMarkerStyle(21);
1844 gq->Draw("pl");
1845
1846 Int_t points=(Int_t)qm*4;
1847 TGraph *gg=new TGraph(points);
1848 i=0;
1849 for(Double_t v=0.05;v<=qm+.05;v+=0.25){
1850 TH1F *dummy=ComputeELossHisto(2,v,len,e);
1851 Double_t avgloss=dummy->GetMean();
1852 gg->SetPoint(i,v,avgloss);i++;
1853 delete dummy;
1854 }
1855 gg->SetMarkerStyle(20);
1856 gg->SetMarkerColor(2);
1857 gg->Draw("pl");
1858
1859 TGraph *gratio=new TGraph(points);
1860 for(Int_t i=0;i<points;i++){
1861 Double_t x,y,x2,y2;
1862 gg->GetPoint(i,x,y);
1863 gq->GetPoint(i,x2,y2);
1864 if(y2>0)
1865 gratio->SetPoint(i,x,y/y2*10/2.25);
1866 else gratio->SetPoint(i,x,0);
1867 }
1868 gratio->SetLineStyle(4);
1869 gratio->Draw();
1870 TLegend *l1a = new TLegend(0.15,0.60,0.50,0.90);
1871 l1a->SetFillStyle(0);
1872 l1a->SetBorderSize(0);
1873 Char_t label[100];
1874 sprintf(label,"L = %.1f fm",len);
1875 l1a->AddEntry(gq,label,"");
1876 l1a->AddEntry(gq,"quark projectile","pl");
1877 l1a->AddEntry(gg,"gluon projectile","pl");
1878 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1879 l1a->Draw();
1880
1881 c->Update();
1882}
1883
1884void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
1885{
1886 // plot average energy loss for given
1887 // length distribution and parton energy
1888
1889 if(!fTablesLoaded){
1890 Error("PlotAvgELossVs","Tables are not loaded.");
1891 return;
1892 }
1893
1894 Char_t title[1024];
1895 Char_t name[1024];
1896 if(fMultSoft){
1897 sprintf(title,"Average Energy Loss for Multiple Scattering");
1898 sprintf(name,"cavgelossms2");
1899 } else {
1900 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1901 sprintf(name,"cavgelosssh2");
1902 }
1903
1904 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1905 c->cd();
1906 TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1907 hframe->SetStats(0);
1908 if(fMultSoft)
1909 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1910 else
1911 hframe->SetXTitle("#mu [GeV]");
1912 hframe->SetYTitle("<E_{loss}> [GeV]");
1913 hframe->Draw();
1914
1915 TGraph *gq=new TGraph(20);
1916 Int_t i=0;
1917 for(Double_t v=0.05;v<=5.05;v+=0.25){
1918 TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1919 Double_t avgloss=dummy->GetMean();
1920 gq->SetPoint(i,v,avgloss);i++;
1921 delete dummy;
1922 }
1923 gq->SetMarkerStyle(20);
1924 gq->Draw("pl");
1925
1926 TGraph *gg=new TGraph(20);
1927 i=0;
1928 for(Double_t v=0.05;v<=5.05;v+=0.25){
1929 TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1930 Double_t avgloss=dummy->GetMean();
1931 gg->SetPoint(i,v,avgloss);i++;
1932 delete dummy;
1933 }
1934 gg->SetMarkerStyle(24);
1935 gg->Draw("pl");
1936
1937 TGraph *gratio=new TGraph(20);
1938 for(Int_t i=0;i<20;i++){
1939 Double_t x,y,x2,y2;
1940 gg->GetPoint(i,x,y);
1941 gq->GetPoint(i,x2,y2);
1942 if(y2>0)
1943 gratio->SetPoint(i,x,y/y2*10/2.25);
1944 else gratio->SetPoint(i,x,0);
1945 }
1946 gratio->SetLineStyle(4);
1947 //gratio->Draw();
1948
1949 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1950 l1a->SetFillStyle(0);
1951 l1a->SetBorderSize(0);
1952 Char_t label[100];
1953 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1954 l1a->AddEntry(gq,label,"");
1955 l1a->AddEntry(gq,"quark","pl");
1956 l1a->AddEntry(gg,"gluon","pl");
1957 //l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1958 l1a->Draw();
1959
1960 c->Update();
1961}
1962
1963void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const
1964{
1965 // plot average energy loss versus ell
1966
1967 if(!fTablesLoaded){
1968 Error("PlotAvgELossVsEll","Tables are not loaded.");
1969 return;
1970 }
1971
1972 Char_t title[1024];
1973 Char_t name[1024];
1974 Float_t medval;
1975 if(fMultSoft){
1976 sprintf(title,"Average Energy Loss for Multiple Scattering");
1977 sprintf(name,"cavgelosslms");
1978 medval=fQTransport;
1979 } else {
1980 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1981 sprintf(name,"cavgelosslsh");
1982 medval=fMu;
1983 }
1984
1985 TCanvas *c = new TCanvas(name,title,0,0,600,400);
1986 c->cd();
1987 TH2F *hframe = new TH2F("avglossell",title,2,0,fLengthMax,2,0,250);
1988 hframe->SetStats(0);
1989 hframe->SetXTitle("length [fm]");
1990 hframe->SetYTitle("<E_{loss}> [GeV]");
1991 hframe->Draw();
1992
1993 TGraph *gq=new TGraph((Int_t)fLengthMax*4);
1994 Int_t i=0;
1995 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
1996 TH1F *dummy=ComputeELossHisto(1,medval,v,e);
1997 Double_t avgloss=dummy->GetMean();
1998 gq->SetPoint(i,v,avgloss);i++;
1999 delete dummy;
2000 }
2001 gq->SetMarkerStyle(20);
2002 gq->Draw("pl");
2003
2004 TGraph *gg=new TGraph((Int_t)fLengthMax*4);
2005 i=0;
2006 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
2007 TH1F *dummy=ComputeELossHisto(2,medval,v,e);
2008 Double_t avgloss=dummy->GetMean();
2009 gg->SetPoint(i,v,avgloss);i++;
2010 delete dummy;
2011 }
2012 gg->SetMarkerStyle(24);
2013 gg->Draw("pl");
2014
2015 TGraph *gratio=new TGraph((Int_t)fLengthMax*4);
2016 for(Int_t i=0;i<=(Int_t)fLengthMax*4;i++){
2017 Double_t x,y,x2,y2;
2018 gg->GetPoint(i,x,y);
2019 gq->GetPoint(i,x2,y2);
2020 if(y2>0)
2021 gratio->SetPoint(i,x,y/y2*100/2.25);
2022 else gratio->SetPoint(i,x,0);
2023 }
2024 gratio->SetLineStyle(1);
2025 gratio->SetLineWidth(2);
2026 gratio->Draw();
2027 TLegend *l1a = new TLegend(0.15,0.65,.60,0.85);
2028 l1a->SetFillStyle(0);
2029 l1a->SetBorderSize(0);
2030 Char_t label[100];
2031 if(fMultSoft)
2032 sprintf(label,"#hat{q} = %.2f GeV^{2}/fm",medval);
2033 else
2034 sprintf(label,"#mu = %.2f GeV",medval);
2035 l1a->AddEntry(gq,label,"");
2036 l1a->AddEntry(gq,"quark","pl");
2037 l1a->AddEntry(gg,"gluon","pl");
2038 l1a->AddEntry(gratio,"gluon/quark/2.25*100","pl");
2039 l1a->Draw();
2040
2041 TF1 *f=new TF1("f","100",0,fLengthMax);
2042 f->SetLineStyle(4);
2043 f->SetLineWidth(1);
2044 f->Draw("same");
2045 c->Update();
2046}
2047
2048void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
2049{
2050 // plot relative energy loss for given
2051 // length and parton energy versus pt
2052
2053 if(!fTablesLoaded){
2054 Error("PlotAvgELossVsPt","Tables are not loaded.");
2055 return;
2056 }
2057
2058 Char_t title[1024];
2059 Char_t name[1024];
2060 if(fMultSoft){
2061 sprintf(title,"Relative Energy Loss for Multiple Scattering");
2062 sprintf(name,"cavgelossvsptms");
2063 } else {
2064 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
2065 sprintf(name,"cavgelossvsptsh");
2066 }
2067
2068 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2069 c->cd();
2070 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2071 hframe->SetStats(0);
2072 hframe->SetXTitle("p_{T} [GeV]");
2073 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2074 hframe->Draw();
2075
2076 TGraph *gq=new TGraph(40);
2077 Int_t i=0;
2078 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2079 TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
2080 Double_t avgloss=dummy->GetMean();
2081 gq->SetPoint(i,pt,avgloss/pt);i++;
2082 delete dummy;
2083 }
2084 gq->SetMarkerStyle(20);
2085 gq->Draw("pl");
2086
2087 TGraph *gg=new TGraph(40);
2088 i=0;
2089 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2090 TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
2091 Double_t avgloss=dummy->GetMean();
2092 gg->SetPoint(i,pt,avgloss/pt);i++;
2093 delete dummy;
2094 }
2095 gg->SetMarkerStyle(24);
2096 gg->Draw("pl");
2097
2098 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2099 l1a->SetFillStyle(0);
2100 l1a->SetBorderSize(0);
2101 Char_t label[100];
2102 sprintf(label,"L = %.1f fm",len);
2103 l1a->AddEntry(gq,label,"");
2104 l1a->AddEntry(gq,"quark","pl");
2105 l1a->AddEntry(gg,"gluon","pl");
2106 l1a->Draw();
2107
2108 c->Update();
2109}
2110
2111void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
2112{
2113 // plot relative energy loss for given
2114 // length distribution and parton energy versus pt
2115
2116 if(!fTablesLoaded){
2117 Error("PlotAvgELossVsPt","Tables are not loaded.");
2118 return;
2119 }
2120
2121 Char_t title[1024];
2122 Char_t name[1024];
2123 if(fMultSoft){
2124 sprintf(title,"Relative Energy Loss for Multiple Scattering");
2125 sprintf(name,"cavgelossvsptms2");
2126 } else {
2127 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
2128 sprintf(name,"cavgelossvsptsh2");
2129 }
2130 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2131 c->cd();
2132 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2133 hframe->SetStats(0);
2134 hframe->SetXTitle("p_{T} [GeV]");
2135 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2136 hframe->Draw();
2137
2138 TGraph *gq=new TGraph(40);
2139 Int_t i=0;
2140 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2141 TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
2142 Double_t avgloss=dummy->GetMean();
2143 gq->SetPoint(i,pt,avgloss/pt);i++;
2144 delete dummy;
2145 }
2146 gq->SetMarkerStyle(20);
2147 gq->Draw("pl");
2148
2149 TGraph *gg=new TGraph(40);
2150 i=0;
2151 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2152 TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
2153 Double_t avgloss=dummy->GetMean();
2154 gg->SetPoint(i,pt,avgloss/pt);i++;
2155 delete dummy;
2156 }
2157 gg->SetMarkerStyle(24);
2158 gg->Draw("pl");
2159
2160 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2161 l1a->SetFillStyle(0);
2162 l1a->SetBorderSize(0);
2163 Char_t label[100];
2164 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
2165 l1a->AddEntry(gq,label,"");
2166 l1a->AddEntry(gq,"quark","pl");
2167 l1a->AddEntry(gg,"gluon","pl");
2168 l1a->Draw();
2169
2170 c->Update();
2171}
2172
2173Int_t AliQuenchingWeights::GetIndex(Double_t len) const
2174{
2175 //get the index according to length
2176 if(len>fLengthMax) len=fLengthMax;
2177
2178 Int_t l=Int_t(len/0.25);
2179 if((len-l*0.25)>0.125) l++;
2180 return l;
2181}
2182