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