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