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