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