1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
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]
26 // Origin: C. Loizides constantinos.loizides@cern.ch
27 // A. Dainese andrea.dainese@pd.infn.it
29 //=================== Added by C. Loizides 27/03/04 ===========================
31 // Added support for k-Quenching, where wc=I1*k and R=2I1^2/I0*k
32 // (see the AliFastGlauber class for definition of I0/I1)
33 //-----------------------------------------------------------------------------
35 #include <Riostream.h>
45 #include "AliQuenchingWeights.h"
49 ClassImp(AliQuenchingWeights)
51 // conversion from fm to GeV^-1: 1 fm = fmGeV GeV^-1
52 const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197;
55 const Double_t AliQuenchingWeights::fgkRMax = 1.e6;
58 const Int_t AliQuenchingWeights::fgkBins = 1300;
59 const Double_t AliQuenchingWeights::fgkMaxBin = 1.3;
61 // counter for histogram labels
62 Int_t AliQuenchingWeights::fgCounter = 0;
65 AliQuenchingWeights::AliQuenchingWeights()
67 fInstanceNumber(fgCounter++),
83 AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a)
85 fInstanceNumber(fgCounter++),
102 fMultSoft=a.GetMultSoft();;
105 fQTransport=a.GetQTransport();
106 fECMethod=(kECMethod)a.GetECMethod();
107 fLengthMax=a.GetLengthMax();
108 fInstanceNumber=fgCounter++;
110 snprintf(name,100, "hhistoqw_%d",fInstanceNumber);
111 fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin);
112 for(Int_t bin=1;bin<=fgkBins;bin++)
113 fHisto->SetBinContent(bin,0.);
115 //Missing in the class is the pathname
116 //to the tables, can be added if needed
119 AliQuenchingWeights::~AliQuenchingWeights()
125 void AliQuenchingWeights::Init()
129 fHisto = new TH1F(Form("hhistoqw_%d",fInstanceNumber),"",fgkBins,0.,fgkMaxBin);
130 for(Int_t bin=1;bin<=fgkBins;bin++)
131 fHisto->SetBinContent(bin,0.);
134 void AliQuenchingWeights::Reset()
136 //reset tables if there were used
139 for(Int_t l=0;l<4*fLengthMaxOld;l++){
140 delete fHistos[0][l];
141 delete fHistos[1][l];
148 void AliQuenchingWeights::SetECMethod(kECMethod type)
150 //set energy constraint method
153 if(fECMethod==kDefault)
154 Info("SetECMethod","Energy Constraint Method set to DEFAULT:\nIf (sampled energy loss > parton energy) then sampled energy loss = parton energy.");
155 else if(fECMethod==kReweight)
156 Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy.");
157 else Info("SetECMethod","Energy Constraint Method set to REWEIGHTCONT:\nRequire sampled energy loss <= parton energy (only implemented for FAST method.");
160 Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
162 // read in tables for multiple scattering approximation
163 // path to continuum and to discrete part
165 fTablesLoaded = kFALSE;
169 snprintf(fname,1024, "%s",gSystem->ExpandPathName(contall));
170 //PH ifstream fincont(fname);
171 fstream fincont(fname,ios::in);
172 #if defined(__HP_aCC) || defined(__DECCXX)
173 if(!fincont.rdbuf()->is_open()) return -1;
175 if(!fincont.is_open()) return -1;
179 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
180 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
181 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>fcaq[13][nn]>>
182 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>fcaq[18][nn]>>
183 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>fcaq[23][nn]>>
184 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>fcaq[28][nn]>>
185 fcaq[29][nn]>>fcaq[30][nn]>>fcaq[31][nn]>>fcaq[32][nn]>>fcaq[33][nn])
192 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
193 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
194 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>fcag[13][nn]>>
195 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>fcag[18][nn]>>
196 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>fcag[23][nn]>>
197 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>fcag[28][nn]>>
198 fcag[29][nn]>>fcag[30][nn]>>fcag[31][nn]>>fcag[32][nn]>>fcag[33][nn])
205 snprintf(fname,1024, "%s",gSystem->ExpandPathName(discall));
206 //PH ifstream findisc(fname);
207 fstream findisc(fname,ios::in);
208 #if defined(__HP_aCC) || defined(__DECCXX)
209 if(!findisc.rdbuf()->is_open()) return -1;
211 if(!findisc.is_open()) return -1;
215 while(findisc>>frrr[nn]>>fdaq[nn]) {
220 while(findisc>>frrrg[nn]>>fdag[nn]) {
225 fTablesLoaded = kTRUE;
230 C***************************************************************************
231 C Quenching Weights for Multiple Soft Scattering
236 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
238 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
241 C This package contains quenching weights for gluon radiation in the
242 C multiple soft scattering approximation.
244 C swqmult returns the quenching weight for a quark (ipart=1) or
245 C a gluon (ipart=2) traversing a medium with transport coeficient q and
246 C length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where
247 C wc=0.5*q*L^2 and w is the energy radiated. The output values are
248 C the continuous and discrete (prefactor of the delta function) parts
249 C of the quenching weights.
251 C In order to use this routine, the files cont.all and disc.all need to be
252 C in the working directory.
254 C An initialization of the tables is needed by doing call initmult before
257 C Please, send us any comment:
259 C urs.wiedemann@cern.ch
260 C carlos.salgado@cern.ch
263 C-------------------------------------------------------------------
265 SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete)
267 REAL*8 xx(400), daq(34), caq(34,261), rrr(34)
268 COMMON /dataqua/ xx, daq, caq, rrr
270 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
271 COMMON /dataglu/ xxg, dag, cag, rrrg
273 REAL*8 rrrr,xxxx, continuous, discrete
275 INTEGER nrlow, nrhigh, nxlow, nxhigh
276 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
277 REAL*8 xfraclow, xfrachigh
288 if (rrin.lt.rrr(nr)) then
300 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
301 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
302 if (rrin.gt.10000d0) then
303 rfraclow = dlog(rrhigh/rrin)/dlog(rrhigh/rrlow)
304 rfrachigh = dlog(rrin/rrlow)/dlog(rrhigh/rrlow)
307 if (ipart.eq.1.and.rrin.ge.rrr(1)) then
314 if (ipart.ne.1.and.rrin.ge.rrrg(1)) then
321 if (xxxx.ge.xx(261)) go to 245
323 nxlow = int(xxin/0.01) + 1
325 xfraclow = (xx(nxhigh)-xxin)/0.01
326 xfrachigh = (xxin - xx(nxlow))/0.01
329 clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh)
330 chigh = xfraclow*caq(nrhigh,nxlow)+xfrachigh*caq(nrhigh,nxhigh)
332 clow = xfraclow*cag(nrlow,nxlow)+xfrachigh*cag(nrlow,nxhigh)
333 chigh = xfraclow*cag(nrhigh,nxlow)+xfrachigh*cag(nrhigh,nxhigh)
336 continuous = rfraclow*clow + rfrachigh*chigh
341 discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh)
343 discrete = rfraclow*dag(nrlow) + rfrachigh*dag(nrhigh)
349 REAL*8 xxq(400), daq(34), caq(34,261), rrr(34)
350 COMMON /dataqua/ xxq, daq, caq, rrr
352 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
353 COMMON /dataglu/ xxg, dag, cag, rrrg
355 OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90)
357 read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn),
358 + caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn),
359 + caq(9,nn), caq(10,nn), caq(11,nn), caq(12,nn),
361 + caq(14,nn), caq(15,nn), caq(16,nn), caq(17,nn),
363 + caq(19,nn), caq(20,nn), caq(21,nn), caq(22,nn),
365 + caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn),
367 + caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn),
368 + caq(33,nn), caq(34,nn)
371 read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn),
372 + cag(4,nn), cag(5,nn), cag(6,nn), cag(7,nn), cag(8,nn),
373 + cag(9,nn), cag(10,nn), cag(11,nn), cag(12,nn),
375 + cag(14,nn), cag(15,nn), cag(16,nn), cag(17,nn),
377 + cag(19,nn), cag(20,nn), cag(21,nn), cag(22,nn),
379 + cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn),
381 + cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn),
382 + cag(33,nn), cag(34,nn)
386 OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91)
388 read (21,*) rrr(nn), daq(nn)
391 read (21,*) rrrg(nn), dag(nn)
396 90 PRINT*, 'input - output error'
397 91 PRINT*, 'input - output error #2'
403 =======================================================================
405 Adapted to ROOT macro by A. Dainese - 13/07/2003
406 Ported to class by C. Loizides - 12/02/2004
407 New version for extended R values added - 06/03/2004
410 Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
411 Double_t &continuous,Double_t &discrete) const
413 // Calculate Multiple Scattering approx.
414 // weights for given parton type,
415 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
421 //read-in data before first call
423 Error("CalcMult","Tables are not loaded.");
427 Error("CalcMult","Tables are not loaded for Multiple Scattering.");
431 Double_t rrin = rrrr;
432 Double_t xxin = xxxx;
434 if(xxin>fxx[260]) return -1;
435 Int_t nxlow = (Int_t)(xxin/0.01) + 1;
436 Int_t nxhigh = nxlow + 1;
437 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.01;
438 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01;
441 if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD
442 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
444 Int_t nrlow=0,nrhigh=0;
445 Double_t rrhigh=0,rrlow=0;
446 for(Int_t nr=1; nr<=34; nr++) {
447 if(rrin<frrr[nr-1]) {
450 rrhigh = frrr[nr-1-1];
460 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
461 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
464 rfraclow = TMath::Log2(rrhigh/rrin)/TMath::Log2(rrhigh/rrlow);
465 rfrachigh = TMath::Log2(rrin/rrlow)/TMath::Log2(rrhigh/rrlow);
467 if((ipart==1) && (rrin>=frrr[0]))
474 if((ipart==2) && (rrin>=frrrg[0]))
482 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
484 Double_t clow=0,chigh=0;
486 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
487 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
489 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
490 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
493 continuous = rfraclow*clow + rfrachigh*chigh;
494 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
495 //rfraclow,clow,rfrachigh,chigh,continuous);
498 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
500 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
506 Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *discall)
508 // read in tables for Single Hard Approx.
509 // path to continuum and to discrete part
511 fTablesLoaded = kFALSE;
515 snprintf(fname, 1024, "%s",gSystem->ExpandPathName(contall));
516 //PH ifstream fincont(fname);
517 fstream fincont(fname,ios::in);
518 #if defined(__HP_aCC) || defined(__DECCXX)
519 if(!fincont.rdbuf()->is_open()) return -1;
521 if(!fincont.is_open()) return -1;
525 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
526 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
527 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>
529 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>
531 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>
533 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>
542 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
543 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
544 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>
546 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>
548 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>
550 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>
558 snprintf(fname, 1024, "%s",gSystem->ExpandPathName(discall));
559 //PH ifstream findisc(fname);
560 fstream findisc(fname,ios::in);
561 #if defined(__HP_aCC) || defined(__DECCXX)
562 if(!findisc.rdbuf()->is_open()) return -1;
564 if(!findisc.is_open()) return -1;
568 while(findisc>>frrr[nn]>>fdaq[nn]) {
573 while(findisc>>frrrg[nn]>>fdag[nn]) {
579 fTablesLoaded = kTRUE;
584 C***************************************************************************
585 C Quenching Weights for Single Hard Scattering
590 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
592 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
595 C This package contains quenching weights for gluon radiation in the
596 C single hard scattering approximation.
598 C swqlin returns the quenching weight for a quark (ipart=1) or
599 C a gluon (ipart=2) traversing a medium with Debye screening mass mu and
600 C length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where
601 C wc=0.5*mu^2*L and w is the energy radiated. The output values are
602 C the continuous and discrete (prefactor of the delta function) parts
603 C of the quenching weights.
605 C In order to use this routine, the files contlin.all and disclin.all
606 C need to be in the working directory.
608 C An initialization of the tables is needed by doing call initlin before
611 C Please, send us any comment:
613 C urs.wiedemann@cern.ch
614 C carlos.salgado@cern.ch
617 C-------------------------------------------------------------------
620 SUBROUTINE swqlin(ipart,rrrr,xxxx,continuous,discrete)
622 REAL*8 xx(400), dalq(30), calq(30,261), rrr(30)
623 COMMON /datalinqua/ xx, dalq, calq, rrr
625 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
626 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
628 REAL*8 rrrr,xxxx, continuous, discrete
630 INTEGER nrlow, nrhigh, nxlow, nxhigh
631 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
632 REAL*8 xfraclow, xfrachigh
638 nxlow = int(xxin/0.038) + 1
640 xfraclow = (xx(nxhigh)-xxin)/0.038
641 xfrachigh = (xxin - xx(nxlow))/0.038
644 if (rrin.lt.rrr(nr)) then
656 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
657 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
660 clow = xfraclow*calq(nrlow,nxlow)+xfrachigh*calq(nrlow,nxhigh)
661 chigh = xfraclow*calq(nrhigh,nxlow)+xfrachigh*calq(nrhigh,nxhigh)
663 clow = xfraclow*calg(nrlow,nxlow)+xfrachigh*calg(nrlow,nxhigh)
664 chigh = xfraclow*calg(nrhigh,nxlow)+xfrachigh*calg(nrhigh,nxhigh)
667 continuous = rfraclow*clow + rfrachigh*chigh
670 discrete = rfraclow*dalq(nrlow) + rfrachigh*dalq(nrhigh)
672 discrete = rfraclow*dalg(nrlow) + rfrachigh*dalg(nrhigh)
678 REAL*8 xxlq(400), dalq(30), calq(30,261), rrr(30)
679 COMMON /datalinqua/ xxlq, dalq, calq, rrr
681 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
682 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
684 OPEN(UNIT=20,FILE='contlin.all',STATUS='OLD',ERR=90)
686 read (20,*) xxlq(nn), calq(1,nn), calq(2,nn), calq(3,nn),
687 + calq(4,nn), calq(5,nn), calq(6,nn), calq(7,nn), calq(8,nn),
688 + calq(9,nn), calq(10,nn), calq(11,nn), calq(12,nn),
690 + calq(14,nn), calq(15,nn), calq(16,nn), calq(17,nn),
692 + calq(19,nn), calq(20,nn), calq(21,nn), calq(22,nn),
694 + calq(24,nn), calq(25,nn), calq(26,nn), calq(27,nn),
696 + calq(29,nn), calq(30,nn)
699 read (20,*) xxlg(nn), calg(1,nn), calg(2,nn), calg(3,nn),
700 + calg(4,nn), calg(5,nn), calg(6,nn), calg(7,nn), calg(8,nn),
701 + calg(9,nn), calg(10,nn), calg(11,nn), calg(12,nn),
703 + calg(14,nn), calg(15,nn), calg(16,nn), calg(17,nn),
705 + calg(19,nn), calg(20,nn), calg(21,nn), calg(22,nn),
707 + calg(24,nn), calg(25,nn), calg(26,nn), calg(27,nn),
709 + calg(29,nn), calg(30,nn)
713 OPEN(UNIT=21,FILE='disclin.all',STATUS='OLD',ERR=91)
715 read (21,*) rrr(nn), dalq(nn)
718 read (21,*) rrrlg(nn), dalg(nn)
723 90 PRINT*, 'input - output error'
724 91 PRINT*, 'input - output error #2'
729 =======================================================================
731 Ported to class by C. Loizides - 17/02/2004
735 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx,
736 Double_t &continuous,Double_t &discrete) const
738 // calculate Single Hard approx.
739 // weights for given parton type,
740 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
742 // read-in data before first call
744 Error("CalcSingleHard","Tables are not loaded.");
748 Error("CalcSingleHard","Tables are not loaded for Single Hard Scattering.");
752 Double_t rrin = rrrr;
753 Double_t xxin = xxxx;
755 Int_t nxlow = (Int_t)(xxin/0.038) + 1;
756 Int_t nxhigh = nxlow + 1;
757 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.038;
758 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.038;
761 if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD
762 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
764 Int_t nrlow=0,nrhigh=0;
765 Double_t rrhigh=0,rrlow=0;
766 for(Int_t nr=1; nr<=30; nr++) {
767 if(rrin<frrr[nr-1]) {
770 rrhigh = frrr[nr-1-1];
780 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
781 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
783 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
785 Double_t clow=0,chigh=0;
787 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
788 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
790 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
791 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
794 continuous = rfraclow*clow + rfrachigh*chigh;
795 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
796 // rfraclow,clow,rfrachigh,chigh,continuous);
799 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
801 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
807 Int_t AliQuenchingWeights::CalcMult(Int_t ipart,
808 Double_t w,Double_t qtransp,Double_t length,
809 Double_t &continuous,Double_t &discrete) const
811 // Calculate Multiple Scattering approx.
812 // weights for given parton type,
813 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
815 Double_t wc=CalcWC(qtransp,length);
816 Double_t rrrr=CalcR(wc,length);
818 return CalcMult(ipart,rrrr,xxxx,continuous,discrete);
821 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
822 Double_t w,Double_t mu,Double_t length,
823 Double_t &continuous,Double_t &discrete) const
825 // calculate Single Hard approx.
826 // weights for given parton type,
827 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
829 Double_t wcbar=CalcWCbar(mu,length);
830 Double_t rrrr=CalcR(wcbar,length);
831 Double_t xxxx=w/wcbar;
832 return CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
835 Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const
837 //calculate r value and
838 //check if it is less then maximum
840 Double_t r = wc*l*fgkConvFmToInvGeV;
842 Warning("CalcR","Value of r = %.2f; should be less than %.2f", r, fgkRMax);
848 Double_t AliQuenchingWeights::CalcRk(Double_t k, Double_t I0, Double_t I1) const
850 //calculate R value and
851 //check if it is less then maximum
853 Double_t r = fgkRMax-1;
857 Warning("CalcRk","Value of r = %.2f; should be less than %.2f",r,fgkRMax);
863 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const
865 // return DeltaE for MS or SH scattering
866 // for given parton type, length and energy
867 // Dependant on ECM (energy constraint method)
868 // e is used to determine where to set bins to zero.
871 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
874 if((ipart<1) || (ipart>2)) {
875 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
879 Int_t l=GetIndex(length);
882 if(fECMethod==kReweight){
886 ret=fHistos[ipart-1][l-1]->GetRandom();
888 Warning("GetELossRandom",
889 "Stopping reweighting; maximum loss assigned after 1e6 trials.");
896 Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
901 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const
903 //return quenched parton energy
904 //for given parton type, length and energy
906 Double_t loss=GetELossRandom(ipart,length,e);
910 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const
912 // return DeltaE for MS or SH scattering
913 // for given parton type, length distribution and energy
914 // Dependant on ECM (energy constraint method)
915 // e is used to determine where to set bins to zero.
918 Warning("GetELossRandom","Pointer to length distribution is NULL.");
921 Double_t ell=hell->GetRandom();
922 return GetELossRandom(ipart,ell,e);
925 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e) const
927 //return quenched parton energy
928 //for given parton type, length distribution and energy
930 Double_t loss=GetELossRandom(ipart,hell,e);
934 Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
936 // return DeltaE for new dynamic version
937 // for given parton type, I0 and I1 value and energy
938 // Dependant on ECM (energy constraint method)
939 // e is used to determine where to set bins to zero.
941 // read-in data before first call
943 Fatal("GetELossRandomK","Tables are not loaded.");
946 if((ipart<1) || (ipart>2)) {
947 Fatal("GetELossRandomK","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
951 Double_t r=CalcRk(I0,I1);
953 Fatal("GetELossRandomK","R should not be negative");
956 Double_t wc=CalcWCk(I1);
958 Fatal("GetELossRandomK","wc should be greater than zero");
961 if(SampleEnergyLoss(ipart,r)!=0){
962 Fatal("GetELossRandomK","Could not sample energy loss");
966 if(fECMethod==kReweight){
970 ret=fHisto->GetRandom();
972 Warning("GetELossRandomK",
973 "Stopping reweighting; maximum loss assigned after 1e6 trials.");
981 Double_t ret=fHisto->GetRandom()*wc;
986 Double_t AliQuenchingWeights::CalcQuenchedEnergyK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
988 //return quenched parton energy
989 //for given parton type, I0 and I1 value and energy
991 Double_t loss=GetELossRandomK(ipart,I0,I1,e);
995 Double_t AliQuenchingWeights::GetELossRandomKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
997 // return DeltaE for new dynamic version
998 // for given parton type, I0 and I1 value and energy
999 // Dependant on ECM (energy constraint method)
1000 // e is used to determine where to set bins to zero.
1001 // method is optimized and should only be used if
1002 // all parameters are well within the bounds.
1003 // read-in data tables before first call
1005 Double_t r=CalcRk(I0,I1);
1010 Double_t wc=CalcWCk(I1);
1015 return GetELossRandomKFastR(ipart,r,wc,e);
1018 Double_t AliQuenchingWeights::GetELossRandomKFastR(Int_t ipart, Double_t r, Double_t wc, Double_t e)
1020 // return DeltaE for new dynamic version
1021 // for given parton type, R and wc value and energy
1022 // Dependant on ECM (energy constraint method)
1023 // e is used to determine where to set bins to zero.
1024 // method is optimized and should only be used if
1025 // all parameters are well within the bounds.
1026 // read-in data tables before first call
1032 Double_t discrete=0.;
1033 Double_t continuous=0.;
1035 Double_t xxxx = fHisto->GetBinCenter(bin);
1037 CalcMult(ipart,r,xxxx,continuous,discrete);
1039 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1042 return 0.; //no energy loss
1044 if (!fHisto) Init();
1046 fHisto->SetBinContent(bin,continuous);
1047 Int_t kbinmax=fHisto->FindBin(e/wc);
1048 if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
1049 if(kbinmax==1) return e; //maximum energy loss
1052 for(bin=2; bin<=kbinmax; bin++) {
1053 xxxx = fHisto->GetBinCenter(bin);
1054 CalcMult(ipart,r,xxxx,continuous,discrete);
1055 fHisto->SetBinContent(bin,continuous);
1058 for(bin=2; bin<=kbinmax; bin++) {
1059 xxxx = fHisto->GetBinCenter(bin);
1060 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1061 fHisto->SetBinContent(bin,continuous);
1065 if(fECMethod==kReweight){
1066 fHisto->SetBinContent(kbinmax+1,0);
1067 fHisto->Fill(0.,discrete*fgkBins/fgkMaxBin);
1068 } else if (fECMethod==kReweightCont) {
1069 fHisto->SetBinContent(kbinmax+1,0);
1070 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1071 fHisto->Scale(1./kdelta*(1-discrete));
1072 fHisto->Fill(0.,discrete);
1074 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1075 Double_t val=discrete*fgkBins/fgkMaxBin;
1076 fHisto->Fill(0.,val);
1077 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1079 for(bin=kbinmax+2; bin<=fgkBins; bin++) {
1080 fHisto->SetBinContent(bin,0);
1082 //cout << kbinmax << " " << discrete << " " << fHisto->Integral() << endl;
1083 Double_t ret=fHisto->GetRandom()*wc;
1088 Double_t AliQuenchingWeights::CalcQuenchedEnergyKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
1090 //return quenched parton energy (fast method)
1091 //for given parton type, I0 and I1 value and energy
1093 Double_t loss=GetELossRandomKFast(ipart,I0,I1,e);
1097 Double_t AliQuenchingWeights::GetDiscreteWeight(Int_t ipart, Double_t I0, Double_t I1)
1099 // return discrete weight
1101 Double_t r=CalcRk(I0,I1);
1105 return GetDiscreteWeightR(ipart,r);
1108 Double_t AliQuenchingWeights::GetDiscreteWeightR(Int_t ipart, Double_t r)
1110 // return discrete weight
1116 Double_t discrete=0.;
1117 Double_t continuous=0.;
1119 Double_t xxxx = fHisto->GetBinCenter(bin);
1121 CalcMult(ipart,r,xxxx,continuous,discrete);
1123 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1127 void AliQuenchingWeights::GetZeroLossProb(Double_t &p,Double_t &prw,Double_t &prwcont,
1128 Int_t ipart,Double_t I0,Double_t I1,Double_t e)
1130 //calculate the probabilty that there is no energy
1131 //loss for different ways of energy constraint
1132 p=1.;prw=1.;prwcont=1.;
1133 Double_t r=CalcRk(I0,I1);
1137 Double_t wc=CalcWCk(I1);
1141 GetZeroLossProbR(p,prw,prwcont,ipart,r,wc,e);
1144 void AliQuenchingWeights::GetZeroLossProbR(Double_t &p,Double_t &prw,Double_t &prwcont,
1145 Int_t ipart, Double_t r,Double_t wc,Double_t e)
1147 //calculate the probabilty that there is no energy
1148 //loss for different ways of energy constraint
1153 Double_t discrete=0.;
1154 Double_t continuous=0.;
1155 if (!fHisto) Init();
1156 Int_t kbinmax=fHisto->FindBin(e/wc);
1157 if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
1159 for(Int_t bin=1; bin<=kbinmax; bin++) {
1160 Double_t xxxx = fHisto->GetBinCenter(bin);
1161 CalcMult(ipart,r,xxxx,continuous,discrete);
1162 fHisto->SetBinContent(bin,continuous);
1165 for(Int_t bin=1; bin<=kbinmax; bin++) {
1166 Double_t xxxx = fHisto->GetBinCenter(bin);
1167 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1168 fHisto->SetBinContent(bin,continuous);
1172 //non-reweighted P(Delta E = 0)
1173 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1174 Double_t val=discrete*fgkBins/fgkMaxBin;
1175 fHisto->Fill(0.,val);
1176 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1177 Double_t hint=fHisto->Integral(1,kbinmax+1);
1178 p=fHisto->GetBinContent(1)/hint;
1181 hint=fHisto->Integral(1,kbinmax);
1182 prw=fHisto->GetBinContent(1)/hint;
1184 Double_t xxxx = fHisto->GetBinCenter(1);
1185 CalcMult(ipart,r,xxxx,continuous,discrete);
1186 fHisto->SetBinContent(1,continuous);
1187 hint=fHisto->Integral(1,kbinmax);
1188 fHisto->Scale(1./hint*(1-discrete));
1189 fHisto->Fill(0.,discrete);
1190 prwcont=fHisto->GetBinContent(1);
1193 Int_t AliQuenchingWeights::SampleEnergyLoss()
1195 // Has to be called to fill the histograms
1197 // For stored values fQTransport loop over
1198 // particle type and length = 1 to fMaxLength (fm)
1199 // to fill energy loss histos
1201 // Take histogram of continuous weights
1202 // Take discrete_weight
1203 // If discrete_weight > 1, put all channels to 0, except channel 1
1204 // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral
1206 // read-in data before first call
1208 Error("SampleEnergyLoss","Tables are not loaded.");
1213 Int_t lmax=CalcLengthMax(fQTransport);
1214 if(fLengthMax>lmax){
1215 Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
1219 Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested.");
1223 fHistos=new TH1F**[2];
1224 fHistos[0]=new TH1F*[4*fLengthMax];
1225 fHistos[1]=new TH1F*[4*fLengthMax];
1226 fLengthMaxOld=fLengthMax; //remember old value in case
1227 //user wants to reset
1230 Char_t meddesc[100];
1232 medvalue=(Int_t)(fQTransport*1000.);
1233 snprintf(meddesc, 100, "MS");
1235 medvalue=(Int_t)(fMu*1000.);
1236 snprintf(meddesc, 100, "SH");
1239 for(Int_t ipart=1;ipart<=2;ipart++){
1240 for(Int_t l=1;l<=4*fLengthMax;l++){
1242 snprintf(hname, 100, "hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
1244 Double_t wc = CalcWC(len);
1245 fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgkBins,0.,fgkMaxBin*wc);
1246 fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
1247 fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
1248 fHistos[ipart-1][l-1]->SetLineColor(4);
1250 Double_t rrrr = CalcR(wc,len);
1251 Double_t discrete=0.;
1252 // loop on histogram channels
1253 for(Int_t bin=1; bin<=fgkBins; bin++) {
1254 Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
1255 Double_t continuous;
1257 CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1259 CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1260 fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
1262 // add discrete part to distribution
1264 for(Int_t bin=2;bin<=fgkBins;bin++)
1265 fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
1267 Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1268 fHistos[ipart-1][l-1]->Fill(0.,val);
1270 Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1271 fHistos[ipart-1][l-1]->Scale(1./hint);
1277 Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t r)
1279 // Sample energy loss directly for one particle type
1280 // choses R (safe it and keep it until next call of function)
1282 // read-in data before first call
1284 Error("SampleEnergyLoss","Tables are not loaded.");
1288 Double_t discrete=0.;
1289 Double_t continuous=0;;
1291 if (!fHisto) Init();
1292 Double_t xxxx = fHisto->GetBinCenter(bin);
1294 CalcMult(ipart,r,xxxx,continuous,discrete);
1296 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1299 fHisto->SetBinContent(1,1.);
1300 for(bin=2;bin<=fgkBins;bin++)
1301 fHisto->SetBinContent(bin,0.);
1305 fHisto->SetBinContent(bin,continuous);
1306 for(bin=2; bin<=fgkBins; bin++) {
1307 xxxx = fHisto->GetBinCenter(bin);
1309 CalcMult(ipart,r,xxxx,continuous,discrete);
1311 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1312 fHisto->SetBinContent(bin,continuous);
1315 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
1316 fHisto->Fill(0.,val);
1317 Double_t hint=fHisto->Integral(1,fgkBins);
1319 fHisto->Scale(1./hint);
1321 //cout << discrete << " " << hint << " " << continuous << endl;
1327 const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) const
1329 //return quenching histograms
1330 //for ipart and length
1333 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
1336 if((ipart<1) || (ipart>2)) {
1337 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
1341 Int_t l=GetIndex(length);
1343 return fHistos[ipart-1][l-1];
1346 TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const
1348 // ipart = 1 for quark, 2 for gluon
1349 // medval a) qtransp = transport coefficient (GeV^2/fm)
1350 // b) mu = Debye mass (GeV)
1351 // length = path length in medium (fm)
1352 // Get from SW tables:
1353 // - continuous weight, as a function of dE/wc
1356 Char_t meddesc[100];
1358 wc=CalcWC(medval,length);
1359 snprintf(meddesc, 100, "MS");
1361 wc=CalcWCbar(medval,length);
1362 snprintf(meddesc, 100, "SH");
1366 snprintf(hname, 100, "hContQWHisto_%s_%d_%d_%d",meddesc,ipart,
1367 (Int_t)(medval*1000.),(Int_t)length);
1369 TH1F *hist = new TH1F("hist",hname,fgkBins,0.,fgkMaxBin*wc);
1370 hist->SetXTitle("#Delta E [GeV]");
1371 hist->SetYTitle("p(#Delta E)");
1372 hist->SetLineColor(4);
1374 Double_t rrrr = CalcR(wc,length);
1375 //loop on histogram channels
1376 for(Int_t bin=1; bin<=fgkBins; bin++) {
1377 Double_t xxxx = hist->GetBinCenter(bin)/wc;
1378 Double_t continuous,discrete;
1380 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1381 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1386 hist->SetBinContent(bin,continuous);
1391 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const
1393 // ipart = 1 for quark, 2 for gluon
1394 // medval a) qtransp = transport coefficient (GeV^2/fm)
1395 // b) mu = Debye mass (GeV)
1396 // length = path length in medium (fm)
1397 // Get from SW tables:
1398 // - continuous weight, as a function of dE/wc
1401 Char_t meddesc[100];
1403 wc=CalcWC(medval,length);
1404 snprintf(meddesc, 100, "MS");
1406 wc=CalcWCbar(medval,length);
1407 snprintf(meddesc, 100, "SH");
1411 snprintf(hname, 100, "hContQWHistox_%s_%d_%d_%d",meddesc,ipart,
1412 (Int_t)(medval*1000.),(Int_t)length);
1414 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
1415 histx->SetXTitle("x = #Delta E/#omega_{c}");
1417 histx->SetYTitle("p(#Delta E/#omega_{c})");
1419 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1420 histx->SetLineColor(4);
1422 Double_t rrrr = CalcR(wc,length);
1423 //loop on histogram channels
1424 for(Int_t bin=1; bin<=fgkBins; bin++) {
1425 Double_t xxxx = histx->GetBinCenter(bin);
1426 Double_t continuous,discrete;
1428 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1429 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1434 histx->SetBinContent(bin,continuous);
1439 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t r) const
1441 // compute P(E) distribution for
1442 // given ipart = 1 for quark, 2 for gluon
1445 Char_t meddesc[100];
1447 snprintf(meddesc, 100, "MS");
1449 snprintf(meddesc, 100, "SH");
1453 snprintf(hname, 100, "hQWHistox_%s_%d_%.2f",meddesc,ipart,r);
1454 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
1455 histx->SetXTitle("x = #Delta E/#omega_{c}");
1457 histx->SetYTitle("p(#Delta E/#omega_{c})");
1459 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1460 histx->SetLineColor(4);
1463 Double_t continuous=0.,discrete=0.;
1464 //loop on histogram channels
1465 for(Int_t bin=1; bin<=fgkBins; bin++) {
1466 Double_t xxxx = histx->GetBinCenter(bin);
1468 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1469 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1474 histx->SetBinContent(bin,continuous);
1477 //add discrete part to distribution
1479 for(Int_t bin=2;bin<=fgkBins;bin++)
1480 histx->SetBinContent(bin,0.);
1482 Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins);
1483 histx->Fill(0.,val);
1485 Double_t hint=histx->Integral(1,fgkBins);
1486 if(hint!=0) histx->Scale(1./hint);
1491 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
1493 // compute energy loss histogram for
1494 // parton type, medium value, length and energy
1496 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1498 dummy->SetQTransport(medval);
1501 dummy->SetMu(medval);
1502 dummy->InitSingleHard();
1504 dummy->SampleEnergyLoss();
1509 snprintf(name, 100, "Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1510 snprintf(hname,100, "hLossQuarks");
1512 snprintf(name, 100, "Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1513 snprintf(hname, 100, "hLossGluons");
1516 TH1F *h = new TH1F(hname,name,250,0,250);
1517 for(Int_t i=0;i<100000;i++){
1518 //if(i % 1000 == 0) cout << "." << flush;
1519 Double_t loss=dummy->GetELossRandom(ipart,l,e);
1527 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
1529 // compute energy loss histogram for
1530 // parton type, medium value,
1531 // length distribution and energy
1533 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1535 dummy->SetQTransport(medval);
1538 dummy->SetMu(medval);
1539 dummy->InitSingleHard();
1541 dummy->SampleEnergyLoss();
1546 snprintf(name, 100, "Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1547 snprintf(hname,100, "hLossQuarks");
1549 snprintf(name,100, "Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1550 snprintf(hname, 100, "hLossGluons");
1553 TH1F *h = new TH1F(hname,name,250,0,250);
1554 for(Int_t i=0;i<100000;i++){
1555 //if(i % 1000 == 0) cout << "." << flush;
1556 Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
1564 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t r) const
1566 // compute energy loss histogram for
1567 // parton type and given R
1569 TH1F *dummy = ComputeQWHistoX(ipart,r);
1570 if(!dummy) return 0;
1573 snprintf(hname, 100, "hELossHistox_%d_%.2f",ipart,r);
1574 TH1F *histx = new TH1F("histxr",hname,fgkBins,0.,fgkMaxBin);
1575 for(Int_t i=0;i<100000;i++){
1576 //if(i % 1000 == 0) cout << "." << flush;
1577 Double_t loss=dummy->GetRandom();
1584 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,Double_t l) const
1586 // compute average energy loss for
1587 // parton type, medium value, length and energy
1589 TH1F *dummy = ComputeELossHisto(ipart,medval,l);
1590 if(!dummy) return 0;
1591 Double_t ret=dummy->GetMean();
1596 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEll) const
1598 // compute average energy loss for
1599 // parton type, medium value,
1600 // length distribution and energy
1602 TH1F *dummy = ComputeELossHisto(ipart,medval,hEll);
1603 if(!dummy) return 0;
1604 Double_t ret=dummy->GetMean();
1609 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t r) const
1611 // compute average energy loss over wc
1612 // for parton type and given R
1614 TH1F *dummy = ComputeELossHisto(ipart,r);
1615 if(!dummy) return 0;
1616 Double_t ret=dummy->GetMean();
1621 void AliQuenchingWeights::PlotDiscreteWeights(Double_t len,Double_t qm) const
1623 // plot discrete weights for given length
1627 c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
1629 c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
1632 TH2F *hframe = new TH2F("hdisc","",2,0,qm+.1,2,0,1.25);
1633 hframe->SetStats(0);
1635 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1637 hframe->SetXTitle("#mu [GeV]");
1638 //hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
1639 hframe->SetYTitle("p_{0} (discrete weight)");
1642 Int_t points=(Int_t)qm*4;
1643 TGraph *gq=new TGraph(points);
1646 for(Double_t q=0.05;q<=qm+.05;q+=0.25){
1648 CalcMult(1,1.0,q,len,cont,disc);
1649 gq->SetPoint(i,q,disc);i++;
1652 for(Double_t m=0.05;m<=qm+.05;m+=0.25){
1654 CalcSingleHard(1,1.0,m,len,cont, disc);
1655 gq->SetPoint(i,m,disc);i++;
1658 gq->SetMarkerStyle(20);
1659 gq->SetMarkerColor(1);
1660 gq->SetLineStyle(1);
1661 gq->SetLineColor(1);
1664 TGraph *gg=new TGraph(points);
1667 for(Double_t q=0.05;q<=qm+.05;q+=0.25){
1669 CalcMult(2,1.0,q,len,cont,disc);
1670 gg->SetPoint(i,q,disc);i++;
1673 for(Double_t m=0.05;m<=qm+.05;m+=0.25){
1675 CalcSingleHard(2,1.0,m,len,cont,disc);
1676 gg->SetPoint(i,m,disc);i++;
1679 gg->SetMarkerStyle(24);
1680 gg->SetMarkerColor(2);
1681 gg->SetLineStyle(2);
1682 gg->SetLineColor(2);
1685 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1686 l1a->SetFillStyle(0);
1687 l1a->SetBorderSize(0);
1689 snprintf(label, 100, "L = %.1f fm",len);
1690 l1a->AddEntry(gq,label,"");
1691 l1a->AddEntry(gq,"quark projectile","l");
1692 l1a->AddEntry(gg,"gluon projectile","l");
1698 void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
1700 // plot continous weights for
1701 // given parton type and length
1708 snprintf(title, 1024, "Cont. Weight for Multiple Scattering - Quarks");
1710 snprintf(title, 1024, "Cont. Weight for Multiple Scattering - Gluons");
1712 medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1713 snprintf(name, 1024, "ccont-ms-%d",itype);
1716 snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Quarks");
1718 snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Gluons");
1720 medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1721 snprintf(name, 1024, "ccont-ms-%d",itype);
1724 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1726 TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1728 h1->SetTitle(title);
1730 h1->SetLineColor(1);
1732 TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1734 h2->SetLineColor(2);
1735 h2->DrawCopy("SAME");
1736 TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1738 h3->SetLineColor(3);
1739 h3->DrawCopy("SAME");
1741 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1742 l1a->SetFillStyle(0);
1743 l1a->SetBorderSize(0);
1745 snprintf(label, 100, "L = %.1f fm",ell);
1746 l1a->AddEntry(h1,label,"");
1748 snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1749 l1a->AddEntry(h1,label,"pl");
1750 snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1751 l1a->AddEntry(h2,label,"pl");
1752 snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1753 l1a->AddEntry(h3, label,"pl");
1755 snprintf(label, 100, "#mu = %.1f GeV",medvals[0]);
1756 l1a->AddEntry(h1,label,"pl");
1757 snprintf(label, 100, "#mu = %.1f GeV",medvals[1]);
1758 l1a->AddEntry(h2,label,"pl");
1759 snprintf(label, 100, "#mu = %.1f GeV",medvals[2]);
1760 l1a->AddEntry(h3,label,"pl");
1767 void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
1769 // plot continous weights for
1770 // given parton type and medium value
1776 snprintf(title,1024, "Cont. Weight for Multiple Scattering - Quarks");
1778 snprintf(title,1024, "Cont. Weight for Multiple Scattering - Gluons");
1780 snprintf(name,1024, "ccont2-ms-%d",itype);
1783 snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Quarks");
1785 snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Gluons");
1787 snprintf(name, 1024, "ccont2-sh-%d",itype);
1789 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1791 TH1F *h1=ComputeQWHisto(itype,medval,8);
1793 h1->SetTitle(title);
1795 h1->SetLineColor(1);
1797 TH1F *h2=ComputeQWHisto(itype,medval,5);
1799 h2->SetLineColor(2);
1800 h2->DrawCopy("SAME");
1801 TH1F *h3=ComputeQWHisto(itype,medval,2);
1803 h3->SetLineColor(3);
1804 h3->DrawCopy("SAME");
1806 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1807 l1a->SetFillStyle(0);
1808 l1a->SetBorderSize(0);
1811 snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medval);
1813 snprintf(label, 100, "#mu = %.1f GeV",medval);
1815 l1a->AddEntry(h1,label,"");
1816 l1a->AddEntry(h1,"L = 8 fm","pl");
1817 l1a->AddEntry(h2,"L = 5 fm","pl");
1818 l1a->AddEntry(h3,"L = 2 fm","pl");
1824 void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t qm,Double_t e) const
1826 // plot average energy loss for given length
1827 // and parton energy
1830 Error("PlotAvgELoss","Tables are not loaded.");
1837 snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
1838 snprintf(name, 1024, "cavgelossms");
1840 snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering");
1841 snprintf(name, 1024, "cavgelosssh");
1844 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1846 TH2F *hframe = new TH2F("avgloss","",2,0,qm+.1,2,0,100);
1847 hframe->SetStats(0);
1849 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1851 hframe->SetXTitle("#mu [GeV]");
1852 hframe->SetYTitle("<E_{loss}> [GeV]");
1855 TGraph *gq=new TGraph(20);
1857 for(Double_t v=0.05;v<=qm+.05;v+=0.25){
1858 TH1F *dummy=ComputeELossHisto(1,v,len,e);
1859 Double_t avgloss=dummy->GetMean();
1860 gq->SetPoint(i,v,avgloss);i++;
1863 gq->SetMarkerStyle(21);
1866 Int_t points=(Int_t)qm*4;
1867 TGraph *gg=new TGraph(points);
1869 for(Double_t v=0.05;v<=qm+.05;v+=0.25){
1870 TH1F *dummy=ComputeELossHisto(2,v,len,e);
1871 Double_t avgloss=dummy->GetMean();
1872 gg->SetPoint(i,v,avgloss);i++;
1875 gg->SetMarkerStyle(20);
1876 gg->SetMarkerColor(2);
1879 TGraph *gratio=new TGraph(points);
1880 for(i=0;i<points;i++){
1882 gg->GetPoint(i,x,y);
1883 gq->GetPoint(i,x2,y2);
1885 gratio->SetPoint(i,x,y/y2*10/2.25);
1886 else gratio->SetPoint(i,x,0);
1888 gratio->SetLineStyle(4);
1890 TLegend *l1a = new TLegend(0.15,0.60,0.50,0.90);
1891 l1a->SetFillStyle(0);
1892 l1a->SetBorderSize(0);
1894 snprintf(label, 100, "L = %.1f fm",len);
1895 l1a->AddEntry(gq,label,"");
1896 l1a->AddEntry(gq,"quark projectile","pl");
1897 l1a->AddEntry(gg,"gluon projectile","pl");
1898 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1904 void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
1906 // plot average energy loss for given
1907 // length distribution and parton energy
1910 Error("PlotAvgELossVs","Tables are not loaded.");
1917 snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
1918 snprintf(name, 1024, "cavgelossms2");
1920 snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering");
1921 snprintf(name, 1024, "cavgelosssh2");
1924 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1926 TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1927 hframe->SetStats(0);
1929 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1931 hframe->SetXTitle("#mu [GeV]");
1932 hframe->SetYTitle("<E_{loss}> [GeV]");
1935 TGraph *gq=new TGraph(20);
1937 for(Double_t v=0.05;v<=5.05;v+=0.25){
1938 TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1939 Double_t avgloss=dummy->GetMean();
1940 gq->SetPoint(i,v,avgloss);i++;
1943 gq->SetMarkerStyle(20);
1946 TGraph *gg=new TGraph(20);
1948 for(Double_t v=0.05;v<=5.05;v+=0.25){
1949 TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1950 Double_t avgloss=dummy->GetMean();
1951 gg->SetPoint(i,v,avgloss);i++;
1954 gg->SetMarkerStyle(24);
1957 TGraph *gratio=new TGraph(20);
1960 gg->GetPoint(i,x,y);
1961 gq->GetPoint(i,x2,y2);
1963 gratio->SetPoint(i,x,y/y2*10/2.25);
1964 else gratio->SetPoint(i,x,0);
1966 gratio->SetLineStyle(4);
1969 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1970 l1a->SetFillStyle(0);
1971 l1a->SetBorderSize(0);
1973 snprintf(label, 100, "<L> = %.2f fm",hEll->GetMean());
1974 l1a->AddEntry(gq,label,"");
1975 l1a->AddEntry(gq,"quark","pl");
1976 l1a->AddEntry(gg,"gluon","pl");
1977 //l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1983 void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const
1985 // plot average energy loss versus ell
1988 Error("PlotAvgELossVsEll","Tables are not loaded.");
1996 snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
1997 snprintf(name, 1024, "cavgelosslms");
2000 snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering");
2001 snprintf(name, 1024, "cavgelosslsh");
2005 TCanvas *c = new TCanvas(name,title,0,0,600,400);
2007 TH2F *hframe = new TH2F("avglossell",title,2,0,fLengthMax,2,0,250);
2008 hframe->SetStats(0);
2009 hframe->SetXTitle("length [fm]");
2010 hframe->SetYTitle("<E_{loss}> [GeV]");
2013 TGraph *gq=new TGraph((Int_t)fLengthMax*4);
2015 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
2016 TH1F *dummy=ComputeELossHisto(1,medval,v,e);
2017 Double_t avgloss=dummy->GetMean();
2018 gq->SetPoint(i,v,avgloss);i++;
2021 gq->SetMarkerStyle(20);
2024 TGraph *gg=new TGraph((Int_t)fLengthMax*4);
2026 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
2027 TH1F *dummy=ComputeELossHisto(2,medval,v,e);
2028 Double_t avgloss=dummy->GetMean();
2029 gg->SetPoint(i,v,avgloss);i++;
2032 gg->SetMarkerStyle(24);
2035 TGraph *gratio=new TGraph((Int_t)fLengthMax*4);
2036 for(i=0;i<=(Int_t)fLengthMax*4;i++){
2038 gg->GetPoint(i,x,y);
2039 gq->GetPoint(i,x2,y2);
2041 gratio->SetPoint(i,x,y/y2*100/2.25);
2042 else gratio->SetPoint(i,x,0);
2044 gratio->SetLineStyle(1);
2045 gratio->SetLineWidth(2);
2047 TLegend *l1a = new TLegend(0.15,0.65,.60,0.85);
2048 l1a->SetFillStyle(0);
2049 l1a->SetBorderSize(0);
2052 snprintf(label, 100, "#hat{q} = %.2f GeV^{2}/fm",medval);
2054 snprintf(label, 100, "#mu = %.2f GeV",medval);
2055 l1a->AddEntry(gq,label,"");
2056 l1a->AddEntry(gq,"quark","pl");
2057 l1a->AddEntry(gg,"gluon","pl");
2058 l1a->AddEntry(gratio,"gluon/quark/2.25*100","pl");
2061 TF1 *f=new TF1("f","100",0,fLengthMax);
2068 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
2070 // plot relative energy loss for given
2071 // length and parton energy versus pt
2074 Error("PlotAvgELossVsPt","Tables are not loaded.");
2081 snprintf(title, 1024, "Relative Energy Loss for Multiple Scattering");
2082 snprintf(name, 1024, "cavgelossvsptms");
2084 snprintf(title, 1024, "Relative Energy Loss for Single Hard Scattering");
2085 snprintf(name, 1024, "cavgelossvsptsh");
2088 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2090 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2091 hframe->SetStats(0);
2092 hframe->SetXTitle("p_{T} [GeV]");
2093 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2096 TGraph *gq=new TGraph(40);
2098 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2099 TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
2100 Double_t avgloss=dummy->GetMean();
2101 gq->SetPoint(i,pt,avgloss/pt);i++;
2104 gq->SetMarkerStyle(20);
2107 TGraph *gg=new TGraph(40);
2109 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2110 TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
2111 Double_t avgloss=dummy->GetMean();
2112 gg->SetPoint(i,pt,avgloss/pt);i++;
2115 gg->SetMarkerStyle(24);
2118 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2119 l1a->SetFillStyle(0);
2120 l1a->SetBorderSize(0);
2122 snprintf(label, 100, "L = %.1f fm",len);
2123 l1a->AddEntry(gq,label,"");
2124 l1a->AddEntry(gq,"quark","pl");
2125 l1a->AddEntry(gg,"gluon","pl");
2131 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
2133 // plot relative energy loss for given
2134 // length distribution and parton energy versus pt
2137 Error("PlotAvgELossVsPt","Tables are not loaded.");
2144 snprintf(title, 1024, "Relative Energy Loss for Multiple Scattering");
2145 snprintf(name, 1024, "cavgelossvsptms2");
2147 snprintf(title, 1024, "Relative Energy Loss for Single Hard Scattering");
2148 snprintf(name, 1024, "cavgelossvsptsh2");
2150 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2152 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2153 hframe->SetStats(0);
2154 hframe->SetXTitle("p_{T} [GeV]");
2155 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2158 TGraph *gq=new TGraph(40);
2160 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2161 TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
2162 Double_t avgloss=dummy->GetMean();
2163 gq->SetPoint(i,pt,avgloss/pt);i++;
2166 gq->SetMarkerStyle(20);
2169 TGraph *gg=new TGraph(40);
2171 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2172 TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
2173 Double_t avgloss=dummy->GetMean();
2174 gg->SetPoint(i,pt,avgloss/pt);i++;
2177 gg->SetMarkerStyle(24);
2180 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2181 l1a->SetFillStyle(0);
2182 l1a->SetBorderSize(0);
2184 snprintf(label, 100, "<L> = %.2f fm",hEll->GetMean());
2185 l1a->AddEntry(gq,label,"");
2186 l1a->AddEntry(gq,"quark","pl");
2187 l1a->AddEntry(gg,"gluon","pl");
2193 Int_t AliQuenchingWeights::GetIndex(Double_t len) const
2195 //get the index according to length
2196 if(len>fLengthMax) len=fLengthMax;
2198 Int_t l=Int_t(len/0.25);
2199 if((len-l*0.25)>0.125) l++;