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