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 **************************************************************************/
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
23 // C.A.Salgado and U.A.Wiedemann, Phys.Rev.D68 (2003) 014008 [hep-ph/0302184]
24 // A.Dainese, Eur.Phys.J.C, in press, [nucl-ex/0312005]
27 // Origin: C. Loizides constantinos.loizides@cern.ch
28 // A. Dainese andrea.dainese@pd.infn.it
30 //=================== Added by C. Loizides 27/03/04 ===========================
32 // Added support for k-Quenching, where wc=I1*k and R=2I1^2/I0*k
33 // (see the AliFastGlauber class for definition of I0/I1)
34 //-----------------------------------------------------------------------------
36 #include <Riostream.h>
45 #include "AliQuenchingWeights.h"
47 ClassImp(AliQuenchingWeights)
49 // conversion from fm to GeV^-1: 1 fm = fmGeV GeV^-1
50 const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197;
53 const Double_t AliQuenchingWeights::fgkRMax = 1.e6;
56 const Int_t AliQuenchingWeights::fgkBins = 1300;
57 const Double_t AliQuenchingWeights::fgkMaxBin = 1.3;
59 // counter for histogram labels
60 Int_t AliQuenchingWeights::fgCounter = 0;
63 AliQuenchingWeights::AliQuenchingWeights()
78 fInstanceNumber=fgCounter++;
80 sprintf(name,"hhistoqw_%d",fInstanceNumber);
81 fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin);
82 for(Int_t bin=1;bin<=fgkBins;bin++)
83 fHisto->SetBinContent(bin,0.);
86 AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a)
94 fMultSoft=a.GetMultSoft();;
97 fQTransport=a.GetQTransport();
98 fECMethod=(kECMethod)a.GetECMethod();
99 fLengthMax=a.GetLengthMax();
100 fInstanceNumber=fgCounter++;
102 sprintf(name,"hhistoqw_%d",fInstanceNumber);
103 fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin);
104 for(Int_t bin=1;bin<=fgkBins;bin++)
105 fHisto->SetBinContent(bin,0.);
107 //Missing in the class is the pathname
108 //to the tables, can be added if needed
111 AliQuenchingWeights::~AliQuenchingWeights()
117 void AliQuenchingWeights::Reset()
119 //reset tables if there were used
122 for(Int_t l=0;l<4*fLengthMaxOld;l++){
123 delete fHistos[0][l];
124 delete fHistos[1][l];
131 void AliQuenchingWeights::SetECMethod(kECMethod type)
133 //set energy constraint method
136 if(fECMethod==kDefault)
137 Info("SetECMethod","Energy Constraint Method set to DEFAULT:\nIf (sampled energy loss > parton energy) then sampled energy loss = parton energy.");
139 Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy.");
142 Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
144 // read in tables for multiple scattering approximation
145 // path to continuum and to discrete part
147 fTablesLoaded = kFALSE;
151 sprintf(fname,"%s",gSystem->ExpandPathName(contall));
152 //PH ifstream fincont(fname);
153 fstream fincont(fname,ios::in);
154 #if defined(__HP_aCC) || defined(__DECCXX)
155 if(!fincont.rdbuf()->is_open()) return -1;
157 if(!fincont.is_open()) return -1;
161 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
162 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
163 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>fcaq[13][nn]>>
164 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>fcaq[18][nn]>>
165 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>fcaq[23][nn]>>
166 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>fcaq[28][nn]>>
167 fcaq[29][nn]>>fcaq[30][nn]>>fcaq[31][nn]>>fcaq[32][nn]>>fcaq[33][nn])
174 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
175 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
176 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>fcag[13][nn]>>
177 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>fcag[18][nn]>>
178 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>fcag[23][nn]>>
179 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>fcag[28][nn]>>
180 fcag[29][nn]>>fcag[30][nn]>>fcag[31][nn]>>fcag[32][nn]>>fcag[33][nn])
187 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
188 //PH ifstream findisc(fname);
189 fstream findisc(fname,ios::in);
190 #if defined(__HP_aCC) || defined(__DECCXX)
191 if(!findisc.rdbuf()->is_open()) return -1;
193 if(!findisc.is_open()) return -1;
197 while(findisc>>frrr[nn]>>fdaq[nn]) {
202 while(findisc>>frrrg[nn]>>fdag[nn]) {
207 fTablesLoaded = kTRUE;
212 C***************************************************************************
213 C Quenching Weights for Multiple Soft Scattering
218 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
220 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
223 C This package contains quenching weights for gluon radiation in the
224 C multiple soft scattering approximation.
226 C swqmult returns the quenching weight for a quark (ipart=1) or
227 C a gluon (ipart=2) traversing a medium with transport coeficient q and
228 C length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where
229 C wc=0.5*q*L^2 and w is the energy radiated. The output values are
230 C the continuous and discrete (prefactor of the delta function) parts
231 C of the quenching weights.
233 C In order to use this routine, the files cont.all and disc.all need to be
234 C in the working directory.
236 C An initialization of the tables is needed by doing call initmult before
239 C Please, send us any comment:
241 C urs.wiedemann@cern.ch
242 C carlos.salgado@cern.ch
245 C-------------------------------------------------------------------
247 SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete)
249 REAL*8 xx(400), daq(34), caq(34,261), rrr(34)
250 COMMON /dataqua/ xx, daq, caq, rrr
252 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
253 COMMON /dataglu/ xxg, dag, cag, rrrg
255 REAL*8 rrrr,xxxx, continuous, discrete
257 INTEGER nrlow, nrhigh, nxlow, nxhigh
258 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
259 REAL*8 xfraclow, xfrachigh
270 if (rrin.lt.rrr(nr)) then
282 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
283 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
284 if (rrin.gt.10000d0) then
285 rfraclow = dlog(rrhigh/rrin)/dlog(rrhigh/rrlow)
286 rfrachigh = dlog(rrin/rrlow)/dlog(rrhigh/rrlow)
289 if (ipart.eq.1.and.rrin.ge.rrr(1)) then
296 if (ipart.ne.1.and.rrin.ge.rrrg(1)) then
303 if (xxxx.ge.xx(261)) go to 245
305 nxlow = int(xxin/0.01) + 1
307 xfraclow = (xx(nxhigh)-xxin)/0.01
308 xfrachigh = (xxin - xx(nxlow))/0.01
311 clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh)
312 chigh = xfraclow*caq(nrhigh,nxlow)+xfrachigh*caq(nrhigh,nxhigh)
314 clow = xfraclow*cag(nrlow,nxlow)+xfrachigh*cag(nrlow,nxhigh)
315 chigh = xfraclow*cag(nrhigh,nxlow)+xfrachigh*cag(nrhigh,nxhigh)
318 continuous = rfraclow*clow + rfrachigh*chigh
323 discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh)
325 discrete = rfraclow*dag(nrlow) + rfrachigh*dag(nrhigh)
331 REAL*8 xxq(400), daq(34), caq(34,261), rrr(34)
332 COMMON /dataqua/ xxq, daq, caq, rrr
334 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
335 COMMON /dataglu/ xxg, dag, cag, rrrg
337 OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90)
339 read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn),
340 + caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn),
341 + caq(9,nn), caq(10,nn), caq(11,nn), caq(12,nn),
343 + caq(14,nn), caq(15,nn), caq(16,nn), caq(17,nn),
345 + caq(19,nn), caq(20,nn), caq(21,nn), caq(22,nn),
347 + caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn),
349 + caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn),
350 + caq(33,nn), caq(34,nn)
353 read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn),
354 + cag(4,nn), cag(5,nn), cag(6,nn), cag(7,nn), cag(8,nn),
355 + cag(9,nn), cag(10,nn), cag(11,nn), cag(12,nn),
357 + cag(14,nn), cag(15,nn), cag(16,nn), cag(17,nn),
359 + cag(19,nn), cag(20,nn), cag(21,nn), cag(22,nn),
361 + cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn),
363 + cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn),
364 + cag(33,nn), cag(34,nn)
368 OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91)
370 read (21,*) rrr(nn), daq(nn)
373 read (21,*) rrrg(nn), dag(nn)
378 90 PRINT*, 'input - output error'
379 91 PRINT*, 'input - output error #2'
385 =======================================================================
387 Adapted to ROOT macro by A. Dainese - 13/07/2003
388 Ported to class by C. Loizides - 12/02/2004
389 New version for extended R values added - 06/03/2004
392 Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
393 Double_t &continuous,Double_t &discrete) const
395 // Calculate Multiple Scattering approx.
396 // weights for given parton type,
397 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
403 //read-in data before first call
405 Error("CalcMult","Tables are not loaded.");
409 Error("CalcMult","Tables are not loaded for Multiple Scattering.");
413 Double_t rrin = rrrr;
414 Double_t xxin = xxxx;
416 if(xxin>fxx[260]) return -1;
417 Int_t nxlow = (Int_t)(xxin/0.01) + 1;
418 Int_t nxhigh = nxlow + 1;
419 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.01;
420 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01;
423 if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD
424 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
426 Int_t nrlow=0,nrhigh=0;
427 Double_t rrhigh=0,rrlow=0;
428 for(Int_t nr=1; nr<=34; nr++) {
429 if(rrin<frrr[nr-1]) {
432 rrhigh = frrr[nr-1-1];
442 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
443 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
446 rfraclow = TMath::Log2(rrhigh/rrin)/TMath::Log2(rrhigh/rrlow);
447 rfrachigh = TMath::Log2(rrin/rrlow)/TMath::Log2(rrhigh/rrlow);
449 if((ipart==1) && (rrin>=frrr[0]))
456 if((ipart==2) && (rrin>=frrrg[0]))
464 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
466 Double_t clow=0,chigh=0;
468 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
469 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
471 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
472 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
475 continuous = rfraclow*clow + rfrachigh*chigh;
476 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
477 //rfraclow,clow,rfrachigh,chigh,continuous);
480 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
482 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
488 Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *discall)
490 // read in tables for Single Hard Approx.
491 // path to continuum and to discrete part
493 fTablesLoaded = kFALSE;
497 sprintf(fname,"%s",gSystem->ExpandPathName(contall));
498 //PH ifstream fincont(fname);
499 fstream fincont(fname,ios::in);
500 #if defined(__HP_aCC) || defined(__DECCXX)
501 if(!fincont.rdbuf()->is_open()) return -1;
503 if(!fincont.is_open()) return -1;
507 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
508 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
509 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>
511 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>
513 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>
515 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>
524 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
525 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
526 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>
528 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>
530 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>
532 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>
540 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
541 //PH ifstream findisc(fname);
542 fstream findisc(fname,ios::in);
543 #if defined(__HP_aCC) || defined(__DECCXX)
544 if(!findisc.rdbuf()->is_open()) return -1;
546 if(!findisc.is_open()) return -1;
550 while(findisc>>frrr[nn]>>fdaq[nn]) {
555 while(findisc>>frrrg[nn]>>fdag[nn]) {
561 fTablesLoaded = kTRUE;
566 C***************************************************************************
567 C Quenching Weights for Single Hard Scattering
572 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
574 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
577 C This package contains quenching weights for gluon radiation in the
578 C single hard scattering approximation.
580 C swqlin returns the quenching weight for a quark (ipart=1) or
581 C a gluon (ipart=2) traversing a medium with Debye screening mass mu and
582 C length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where
583 C wc=0.5*mu^2*L and w is the energy radiated. The output values are
584 C the continuous and discrete (prefactor of the delta function) parts
585 C of the quenching weights.
587 C In order to use this routine, the files contlin.all and disclin.all
588 C need to be in the working directory.
590 C An initialization of the tables is needed by doing call initlin before
593 C Please, send us any comment:
595 C urs.wiedemann@cern.ch
596 C carlos.salgado@cern.ch
599 C-------------------------------------------------------------------
602 SUBROUTINE swqlin(ipart,rrrr,xxxx,continuous,discrete)
604 REAL*8 xx(400), dalq(30), calq(30,261), rrr(30)
605 COMMON /datalinqua/ xx, dalq, calq, rrr
607 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
608 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
610 REAL*8 rrrr,xxxx, continuous, discrete
612 INTEGER nrlow, nrhigh, nxlow, nxhigh
613 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
614 REAL*8 xfraclow, xfrachigh
620 nxlow = int(xxin/0.038) + 1
622 xfraclow = (xx(nxhigh)-xxin)/0.038
623 xfrachigh = (xxin - xx(nxlow))/0.038
626 if (rrin.lt.rrr(nr)) then
638 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
639 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
642 clow = xfraclow*calq(nrlow,nxlow)+xfrachigh*calq(nrlow,nxhigh)
643 chigh = xfraclow*calq(nrhigh,nxlow)+xfrachigh*calq(nrhigh,nxhigh)
645 clow = xfraclow*calg(nrlow,nxlow)+xfrachigh*calg(nrlow,nxhigh)
646 chigh = xfraclow*calg(nrhigh,nxlow)+xfrachigh*calg(nrhigh,nxhigh)
649 continuous = rfraclow*clow + rfrachigh*chigh
652 discrete = rfraclow*dalq(nrlow) + rfrachigh*dalq(nrhigh)
654 discrete = rfraclow*dalg(nrlow) + rfrachigh*dalg(nrhigh)
660 REAL*8 xxlq(400), dalq(30), calq(30,261), rrr(30)
661 COMMON /datalinqua/ xxlq, dalq, calq, rrr
663 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
664 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
666 OPEN(UNIT=20,FILE='contlin.all',STATUS='OLD',ERR=90)
668 read (20,*) xxlq(nn), calq(1,nn), calq(2,nn), calq(3,nn),
669 + calq(4,nn), calq(5,nn), calq(6,nn), calq(7,nn), calq(8,nn),
670 + calq(9,nn), calq(10,nn), calq(11,nn), calq(12,nn),
672 + calq(14,nn), calq(15,nn), calq(16,nn), calq(17,nn),
674 + calq(19,nn), calq(20,nn), calq(21,nn), calq(22,nn),
676 + calq(24,nn), calq(25,nn), calq(26,nn), calq(27,nn),
678 + calq(29,nn), calq(30,nn)
681 read (20,*) xxlg(nn), calg(1,nn), calg(2,nn), calg(3,nn),
682 + calg(4,nn), calg(5,nn), calg(6,nn), calg(7,nn), calg(8,nn),
683 + calg(9,nn), calg(10,nn), calg(11,nn), calg(12,nn),
685 + calg(14,nn), calg(15,nn), calg(16,nn), calg(17,nn),
687 + calg(19,nn), calg(20,nn), calg(21,nn), calg(22,nn),
689 + calg(24,nn), calg(25,nn), calg(26,nn), calg(27,nn),
691 + calg(29,nn), calg(30,nn)
695 OPEN(UNIT=21,FILE='disclin.all',STATUS='OLD',ERR=91)
697 read (21,*) rrr(nn), dalq(nn)
700 read (21,*) rrrlg(nn), dalg(nn)
705 90 PRINT*, 'input - output error'
706 91 PRINT*, 'input - output error #2'
711 =======================================================================
713 Ported to class by C. Loizides - 17/02/2004
717 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx,
718 Double_t &continuous,Double_t &discrete) const
720 // calculate Single Hard approx.
721 // weights for given parton type,
722 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
724 // read-in data before first call
726 Error("CalcSingleHard","Tables are not loaded.");
730 Error("CalcSingleHard","Tables are not loaded for Single Hard Scattering.");
734 Double_t rrin = rrrr;
735 Double_t xxin = xxxx;
737 Int_t nxlow = (Int_t)(xxin/0.038) + 1;
738 Int_t nxhigh = nxlow + 1;
739 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.038;
740 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.038;
743 if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD
744 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
746 Int_t nrlow=0,nrhigh=0;
747 Double_t rrhigh=0,rrlow=0;
748 for(Int_t nr=1; nr<=30; nr++) {
749 if(rrin<frrr[nr-1]) {
752 rrhigh = frrr[nr-1-1];
762 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
763 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
765 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
767 Double_t clow=0,chigh=0;
769 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
770 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
772 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
773 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
776 continuous = rfraclow*clow + rfrachigh*chigh;
777 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
778 // rfraclow,clow,rfrachigh,chigh,continuous);
781 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
783 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
789 Int_t AliQuenchingWeights::CalcMult(Int_t ipart,
790 Double_t w,Double_t qtransp,Double_t length,
791 Double_t &continuous,Double_t &discrete) const
793 // Calculate Multiple Scattering approx.
794 // weights for given parton type,
795 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
797 Double_t wc=CalcWC(qtransp,length);
798 Double_t rrrr=CalcR(wc,length);
800 return CalcMult(ipart,rrrr,xxxx,continuous,discrete);
803 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
804 Double_t w,Double_t mu,Double_t length,
805 Double_t &continuous,Double_t &discrete) const
807 // calculate Single Hard approx.
808 // weights for given parton type,
809 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
811 Double_t wcbar=CalcWCbar(mu,length);
812 Double_t rrrr=CalcR(wcbar,length);
813 Double_t xxxx=w/wcbar;
814 return CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
817 Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const
819 //calculate R value and
820 //check if it is less then maximum
822 Double_t R = wc*l*fgkConvFmToInvGeV;
824 Warning("CalcR","Value of R = %.2f; should be less than %.2f",R,fgkRMax);
830 Double_t AliQuenchingWeights::CalcRk(Double_t k, Double_t I0, Double_t I1) const
832 //calculate R value and
833 //check if it is less then maximum
835 Double_t R = fgkRMax;
839 Warning("CalcRk","Value of R = %.2f; should be less than %.2f",R,fgkRMax);
845 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const
847 // return DeltaE for MS or SH scattering
848 // for given parton type, length and energy
849 // Dependant on ECM (energy constraint method)
850 // e is used to determine where to set bins to zero.
853 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
856 if((ipart<1) || (ipart>2)) {
857 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
861 Int_t l=GetIndex(length);
864 if(fECMethod==kReweight){
868 ret=fHistos[ipart-1][l-1]->GetRandom();
870 Warning("GetELossRandom",
871 "Aborted reweighting; maximum loss assigned after 1e6 trials.");
878 Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
883 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const
885 //return quenched parton energy
886 //for given parton type, length and energy
888 Double_t loss=GetELossRandom(ipart,length,e);
892 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const
894 // return DeltaE for MS or SH scattering
895 // for given parton type, length distribution and energy
896 // Dependant on ECM (energy constraint method)
897 // e is used to determine where to set bins to zero.
900 Warning("GetELossRandom","Pointer to length distribution is NULL.");
903 Double_t ell=hell->GetRandom();
904 return GetELossRandom(ipart,ell,e);
907 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e) const
909 //return quenched parton energy
910 //for given parton type, length distribution and energy
912 Double_t loss=GetELossRandom(ipart,hell,e);
916 Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
918 // return DeltaE for new dynamic version
919 // for given parton type, I0 and I1 value and energy
920 // Dependant on ECM (energy constraint method)
921 // e is used to determine where to set bins to zero.
923 // read-in data before first call
925 Fatal("GetELossRandomK","Tables are not loaded.");
928 if((ipart<1) || (ipart>2)) {
929 Fatal("GetELossRandomK","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
933 Double_t R=CalcRk(I0,I1);
935 Fatal("GetELossRandomK","R should not be negative");
938 Double_t wc=CalcWCk(I1);
940 Fatal("GetELossRandomK","wc should be greater than zero");
943 if(SampleEnergyLoss(ipart,R)!=0){
944 Fatal("GetELossRandomK","Could not sample energy loss");
948 if(fECMethod==kReweight){
952 ret=fHisto->GetRandom();
954 Warning("GetELossRandomK",
955 "Aborted reweighting; maximum loss assigned after 1e6 trials.");
963 Double_t ret=fHisto->GetRandom()*wc;
968 Double_t AliQuenchingWeights::CalcQuenchedEnergyK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
970 //return quenched parton energy
971 //for given parton type, I0 and I1 value and energy
973 Double_t loss=GetELossRandomK(ipart,I0,I1,e);
977 Double_t AliQuenchingWeights::GetELossRandomKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
979 // return DeltaE for new dynamic version
980 // for given parton type, I0 and I1 value and energy
981 // Dependant on ECM (energy constraint method)
982 // e is used to determine where to set bins to zero.
983 // method is optimized and should only be used if
984 // all parameters are well within the bounds.
985 // read-in data tables before first call
987 Double_t R=CalcRk(I0,I1);
992 Double_t wc=CalcWCk(I1);
997 Double_t discrete=0.;
998 Double_t continuous=0;
1000 Double_t xxxx = fHisto->GetBinCenter(bin);
1002 CalcMult(ipart,R,xxxx,continuous,discrete);
1004 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1006 if(discrete>0.999) {
1010 fHisto->SetBinContent(bin,continuous);
1011 #if 1 /* fast new version*/
1012 const Int_t kbinmax=fHisto->FindBin(e/wc);
1014 for(Int_t bin=2; bin<=kbinmax; bin++) {
1015 xxxx = fHisto->GetBinCenter(bin);
1016 CalcMult(ipart,R,xxxx,continuous,discrete);
1017 fHisto->SetBinContent(bin,continuous);
1020 for(Int_t bin=2; bin<=kbinmax; bin++) {
1021 xxxx = fHisto->GetBinCenter(bin);
1022 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1023 fHisto->SetBinContent(bin,continuous);
1027 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1029 Double_t val=discrete*fgkBins/fgkMaxBin;
1030 fHisto->Fill(0.,val);
1032 if(fECMethod==kReweight)
1033 fHisto->SetBinContent(kbinmax+1,0);
1035 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1037 for(Int_t bin=kbinmax+2; bin<=fgkBins; bin++) {
1038 fHisto->SetBinContent(bin,0);
1043 for(Int_t bin=2; bin<=fgkBins; bin++) {
1044 xxxx = fHisto->GetBinCenter(bin);
1045 CalcMult(ipart,R,xxxx,continuous,discrete);
1046 fHisto->SetBinContent(bin,continuous);
1049 for(Int_t bin=2; bin<=fgkBins; bin++) {
1050 xxxx = fHisto->GetBinCenter(bin);
1051 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1052 fHisto->SetBinContent(bin,continuous);
1056 // add discrete part to distribution
1057 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
1058 fHisto->Fill(0.,val);
1060 if(fECMethod==kReweight) {
1061 for(Int_t bin=fHisto->FindBin(e/wc)+1; bin<=fgkBins; bin++) {
1062 fHisto->SetBinContent(bin,0);
1067 Double_t ret=fHisto->GetRandom()*wc;
1072 Double_t AliQuenchingWeights::CalcQuenchedEnergyKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
1074 //return quenched parton energy (fast method)
1075 //for given parton type, I0 and I1 value and energy
1077 Double_t loss=GetELossRandomKFast(ipart,I0,I1,e);
1081 Int_t AliQuenchingWeights::SampleEnergyLoss()
1083 // Has to be called to fill the histograms
1085 // For stored values fQTransport loop over
1086 // particle type and length = 1 to fMaxLength (fm)
1087 // to fill energy loss histos
1089 // Take histogram of continuous weights
1090 // Take discrete_weight
1091 // If discrete_weight > 1, put all channels to 0, except channel 1
1092 // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral
1094 // read-in data before first call
1096 Error("SampleEnergyLoss","Tables are not loaded.");
1101 Int_t lmax=CalcLengthMax(fQTransport);
1102 if(fLengthMax>lmax){
1103 Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
1107 Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested.");
1111 fHistos=new TH1F**[2];
1112 fHistos[0]=new TH1F*[4*fLengthMax];
1113 fHistos[1]=new TH1F*[4*fLengthMax];
1114 fLengthMaxOld=fLengthMax; //remember old value in case
1115 //user wants to reset
1118 Char_t meddesc[100];
1120 medvalue=(Int_t)(fQTransport*1000.);
1121 sprintf(meddesc,"MS");
1123 medvalue=(Int_t)(fMu*1000.);
1124 sprintf(meddesc,"SH");
1127 for(Int_t ipart=1;ipart<=2;ipart++){
1128 for(Int_t l=1;l<=4*fLengthMax;l++){
1130 sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
1132 Double_t wc = CalcWC(len);
1133 fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgkBins,0.,fgkMaxBin*wc);
1134 fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
1135 fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
1136 fHistos[ipart-1][l-1]->SetLineColor(4);
1138 Double_t rrrr = CalcR(wc,len);
1139 Double_t discrete=0.;
1140 // loop on histogram channels
1141 for(Int_t bin=1; bin<=fgkBins; bin++) {
1142 Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
1143 Double_t continuous;
1145 CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1147 CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1148 fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
1150 // add discrete part to distribution
1152 for(Int_t bin=2;bin<=fgkBins;bin++)
1153 fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
1155 Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1156 fHistos[ipart-1][l-1]->Fill(0.,val);
1158 Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1159 fHistos[ipart-1][l-1]->Scale(1./hint);
1165 Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t R)
1167 // Sample energy loss directly for one particle type
1168 // choses R (safe it and keep it until next call of function)
1170 // read-in data before first call
1172 Error("SampleEnergyLoss","Tables are not loaded.");
1176 Double_t discrete=0.;
1177 Double_t continuous=0;;
1179 Double_t xxxx = fHisto->GetBinCenter(bin);
1181 CalcMult(ipart,R,xxxx,continuous,discrete);
1183 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1186 fHisto->SetBinContent(1,1.);
1187 for(Int_t bin=2;bin<=fgkBins;bin++)
1188 fHisto->SetBinContent(bin,0.);
1192 fHisto->SetBinContent(bin,continuous);
1193 for(Int_t bin=2; bin<=fgkBins; bin++) {
1194 xxxx = fHisto->GetBinCenter(bin);
1196 CalcMult(ipart,R,xxxx,continuous,discrete);
1198 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1199 fHisto->SetBinContent(bin,continuous);
1202 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
1203 fHisto->Fill(0.,val);
1204 Double_t hint=fHisto->Integral(1,fgkBins);
1206 fHisto->Scale(1./hint);
1208 //cout << discrete << " " << hint << " " << continuous << endl;
1214 const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) const
1216 //return quenching histograms
1217 //for ipart and length
1220 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
1223 if((ipart<1) || (ipart>2)) {
1224 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
1228 Int_t l=GetIndex(length);
1230 return fHistos[ipart-1][l-1];
1233 TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const
1235 // ipart = 1 for quark, 2 for gluon
1236 // medval a) qtransp = transport coefficient (GeV^2/fm)
1237 // b) mu = Debye mass (GeV)
1238 // length = path length in medium (fm)
1239 // Get from SW tables:
1240 // - continuous weight, as a function of dE/wc
1243 Char_t meddesc[100];
1245 wc=CalcWC(medval,length);
1246 sprintf(meddesc,"MS");
1248 wc=CalcWCbar(medval,length);
1249 sprintf(meddesc,"SH");
1253 sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart,
1254 (Int_t)(medval*1000.),(Int_t)length);
1256 TH1F *hist = new TH1F("hist",hname,fgkBins,0.,fgkMaxBin*wc);
1257 hist->SetXTitle("#Delta E [GeV]");
1258 hist->SetYTitle("p(#Delta E)");
1259 hist->SetLineColor(4);
1261 Double_t rrrr = CalcR(wc,length);
1262 //loop on histogram channels
1263 for(Int_t bin=1; bin<=fgkBins; bin++) {
1264 Double_t xxxx = hist->GetBinCenter(bin)/wc;
1265 Double_t continuous,discrete;
1267 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1268 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1273 hist->SetBinContent(bin,continuous);
1278 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const
1280 // ipart = 1 for quark, 2 for gluon
1281 // medval a) qtransp = transport coefficient (GeV^2/fm)
1282 // b) mu = Debye mass (GeV)
1283 // length = path length in medium (fm)
1284 // Get from SW tables:
1285 // - continuous weight, as a function of dE/wc
1288 Char_t meddesc[100];
1290 wc=CalcWC(medval,length);
1291 sprintf(meddesc,"MS");
1293 wc=CalcWCbar(medval,length);
1294 sprintf(meddesc,"SH");
1298 sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart,
1299 (Int_t)(medval*1000.),(Int_t)length);
1301 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
1302 histx->SetXTitle("x = #Delta E/#omega_{c}");
1304 histx->SetYTitle("p(#Delta E/#omega_{c})");
1306 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1307 histx->SetLineColor(4);
1309 Double_t rrrr = CalcR(wc,length);
1310 //loop on histogram channels
1311 for(Int_t bin=1; bin<=fgkBins; bin++) {
1312 Double_t xxxx = histx->GetBinCenter(bin);
1313 Double_t continuous,discrete;
1315 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1316 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1321 histx->SetBinContent(bin,continuous);
1326 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t R) const
1328 // compute P(E) distribution for
1329 // given ipart = 1 for quark, 2 for gluon
1332 Char_t meddesc[100];
1334 sprintf(meddesc,"MS");
1336 sprintf(meddesc,"SH");
1340 sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,R);
1341 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
1342 histx->SetXTitle("x = #Delta E/#omega_{c}");
1344 histx->SetYTitle("p(#Delta E/#omega_{c})");
1346 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1347 histx->SetLineColor(4);
1350 Double_t continuous=0.,discrete=0.;
1351 //loop on histogram channels
1352 for(Int_t bin=1; bin<=fgkBins; bin++) {
1353 Double_t xxxx = histx->GetBinCenter(bin);
1355 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1356 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1361 histx->SetBinContent(bin,continuous);
1364 //add discrete part to distribution
1366 for(Int_t bin=2;bin<=fgkBins;bin++)
1367 histx->SetBinContent(bin,0.);
1369 Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins);
1370 histx->Fill(0.,val);
1372 Double_t hint=histx->Integral(1,fgkBins);
1373 if(hint!=0) histx->Scale(1./hint);
1378 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
1380 // compute energy loss histogram for
1381 // parton type, medium value, length and energy
1383 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1385 dummy->SetQTransport(medval);
1388 dummy->SetMu(medval);
1389 dummy->InitSingleHard();
1391 dummy->SampleEnergyLoss();
1396 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1397 sprintf(hname,"hLossQuarks");
1399 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1400 sprintf(hname,"hLossGluons");
1403 TH1F *h = new TH1F(hname,name,250,0,250);
1404 for(Int_t i=0;i<100000;i++){
1405 //if(i % 1000 == 0) cout << "." << flush;
1406 Double_t loss=dummy->GetELossRandom(ipart,l,e);
1414 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
1416 // compute energy loss histogram for
1417 // parton type, medium value,
1418 // length distribution and energy
1420 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1422 dummy->SetQTransport(medval);
1425 dummy->SetMu(medval);
1426 dummy->InitSingleHard();
1428 dummy->SampleEnergyLoss();
1433 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1434 sprintf(hname,"hLossQuarks");
1436 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1437 sprintf(hname,"hLossGluons");
1440 TH1F *h = new TH1F(hname,name,250,0,250);
1441 for(Int_t i=0;i<100000;i++){
1442 //if(i % 1000 == 0) cout << "." << flush;
1443 Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
1451 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t R) const
1453 // compute energy loss histogram for
1454 // parton type and given R
1456 TH1F *dummy = ComputeQWHistoX(ipart,R);
1457 if(!dummy) return 0;
1460 sprintf(hname,"hELossHistox_%d_%.2f",ipart,R);
1461 TH1F *histx = new TH1F("histxr",hname,fgkBins,0.,fgkMaxBin);
1462 for(Int_t i=0;i<100000;i++){
1463 //if(i % 1000 == 0) cout << "." << flush;
1464 Double_t loss=dummy->GetRandom();
1471 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,Double_t l) const
1473 // compute average energy loss for
1474 // parton type, medium value, length and energy
1476 TH1F *dummy = ComputeELossHisto(ipart,medval,l);
1477 if(!dummy) return 0;
1478 Double_t ret=dummy->GetMean();
1483 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEll) const
1485 // compute average energy loss for
1486 // parton type, medium value,
1487 // length distribution and energy
1489 TH1F *dummy = ComputeELossHisto(ipart,medval,hEll);
1490 if(!dummy) return 0;
1491 Double_t ret=dummy->GetMean();
1496 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t R) const
1498 // compute average energy loss over wc
1499 // for parton type and given R
1501 TH1F *dummy = ComputeELossHisto(ipart,R);
1502 if(!dummy) return 0;
1503 Double_t ret=dummy->GetMean();
1508 void AliQuenchingWeights::PlotDiscreteWeights(Double_t len) const
1510 // plot discrete weights for given length
1514 c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
1516 c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
1519 TH2F *hframe = new TH2F("hdisc","",2,0,5.1,2,0,1);
1520 hframe->SetStats(0);
1522 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1524 hframe->SetXTitle("#mu [GeV]");
1525 hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
1528 TGraph *gq=new TGraph(20);
1531 for(Double_t q=0.05;q<=5.05;q+=0.25){
1533 CalcMult(1,1.0,q,len,cont,disc);
1534 gq->SetPoint(i,q,disc);i++;
1537 for(Double_t m=0.05;m<=5.05;m+=0.25){
1539 CalcSingleHard(1,1.0,m,len,cont, disc);
1540 gq->SetPoint(i,m,disc);i++;
1543 gq->SetMarkerStyle(20);
1546 TGraph *gg=new TGraph(20);
1549 for(Double_t q=0.05;q<=5.05;q+=0.25){
1551 CalcMult(2,1.0,q,len,cont,disc);
1552 gg->SetPoint(i,q,disc);i++;
1555 for(Double_t m=0.05;m<=5.05;m+=0.25){
1557 CalcSingleHard(2,1.0,m,len,cont,disc);
1558 gg->SetPoint(i,m,disc);i++;
1561 gg->SetMarkerStyle(24);
1564 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1565 l1a->SetFillStyle(0);
1566 l1a->SetBorderSize(0);
1568 sprintf(label,"L = %.1f fm",len);
1569 l1a->AddEntry(gq,label,"");
1570 l1a->AddEntry(gq,"quark","pl");
1571 l1a->AddEntry(gg,"gluon","pl");
1577 void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
1579 // plot continous weights for
1580 // given parton type and length
1587 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1589 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1591 medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1592 sprintf(name,"ccont-ms-%d",itype);
1595 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1597 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1599 medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1600 sprintf(name,"ccont-ms-%d",itype);
1603 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1605 TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1607 h1->SetTitle(title);
1609 h1->SetLineColor(1);
1611 TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1613 h2->SetLineColor(2);
1614 h2->DrawCopy("SAME");
1615 TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1617 h3->SetLineColor(3);
1618 h3->DrawCopy("SAME");
1620 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1621 l1a->SetFillStyle(0);
1622 l1a->SetBorderSize(0);
1624 sprintf(label,"L = %.1f fm",ell);
1625 l1a->AddEntry(h1,label,"");
1627 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1628 l1a->AddEntry(h1,label,"pl");
1629 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1630 l1a->AddEntry(h2,label,"pl");
1631 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1632 l1a->AddEntry(h3,label,"pl");
1634 sprintf(label,"#mu = %.1f GeV",medvals[0]);
1635 l1a->AddEntry(h1,label,"pl");
1636 sprintf(label,"#mu = %.1f GeV",medvals[1]);
1637 l1a->AddEntry(h2,label,"pl");
1638 sprintf(label,"#mu = %.1f GeV",medvals[2]);
1639 l1a->AddEntry(h3,label,"pl");
1646 void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
1648 // plot continous weights for
1649 // given parton type and medium value
1655 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1657 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1659 sprintf(name,"ccont2-ms-%d",itype);
1662 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1664 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1666 sprintf(name,"ccont2-sh-%d",itype);
1668 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1670 TH1F *h1=ComputeQWHisto(itype,medval,8);
1672 h1->SetTitle(title);
1674 h1->SetLineColor(1);
1676 TH1F *h2=ComputeQWHisto(itype,medval,5);
1678 h2->SetLineColor(2);
1679 h2->DrawCopy("SAME");
1680 TH1F *h3=ComputeQWHisto(itype,medval,2);
1682 h3->SetLineColor(3);
1683 h3->DrawCopy("SAME");
1685 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1686 l1a->SetFillStyle(0);
1687 l1a->SetBorderSize(0);
1690 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval);
1692 sprintf(label,"#mu = %.1f GeV",medval);
1694 l1a->AddEntry(h1,label,"");
1695 l1a->AddEntry(h1,"L = 8 fm","pl");
1696 l1a->AddEntry(h2,"L = 5 fm","pl");
1697 l1a->AddEntry(h3,"L = 2 fm","pl");
1703 void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t e) const
1705 // plot average energy loss for given length
1706 // and parton energy
1709 Error("PlotAvgELoss","Tables are not loaded.");
1716 sprintf(title,"Average Energy Loss for Multiple Scattering");
1717 sprintf(name,"cavgelossms");
1719 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1720 sprintf(name,"cavgelosssh");
1723 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1725 TH2F *hframe = new TH2F("avgloss",title,2,0,5.1,2,0,100);
1726 hframe->SetStats(0);
1728 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1730 hframe->SetXTitle("#mu [GeV]");
1731 hframe->SetYTitle("<E_{loss}> [GeV]");
1734 TGraph *gq=new TGraph(20);
1736 for(Double_t v=0.05;v<=5.05;v+=0.25){
1737 TH1F *dummy=ComputeELossHisto(1,v,len,e);
1738 Double_t avgloss=dummy->GetMean();
1739 gq->SetPoint(i,v,avgloss);i++;
1742 gq->SetMarkerStyle(20);
1745 TGraph *gg=new TGraph(20);
1747 for(Double_t v=0.05;v<=5.05;v+=0.25){
1748 TH1F *dummy=ComputeELossHisto(2,v,len,e);
1749 Double_t avgloss=dummy->GetMean();
1750 gg->SetPoint(i,v,avgloss);i++;
1753 gg->SetMarkerStyle(24);
1756 TGraph *gratio=new TGraph(20);
1757 for(Int_t i=0;i<20;i++){
1759 gg->GetPoint(i,x,y);
1760 gq->GetPoint(i,x2,y2);
1762 gratio->SetPoint(i,x,y/y2*10/2.25);
1763 else gratio->SetPoint(i,x,0);
1765 gratio->SetLineStyle(4);
1767 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1768 l1a->SetFillStyle(0);
1769 l1a->SetBorderSize(0);
1771 sprintf(label,"L = %.1f fm",len);
1772 l1a->AddEntry(gq,label,"");
1773 l1a->AddEntry(gq,"quark","pl");
1774 l1a->AddEntry(gg,"gluon","pl");
1775 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1781 void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
1783 // plot average energy loss for given
1784 // length distribution and parton energy
1787 Error("PlotAvgELossVs","Tables are not loaded.");
1794 sprintf(title,"Average Energy Loss for Multiple Scattering");
1795 sprintf(name,"cavgelossms2");
1797 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1798 sprintf(name,"cavgelosssh2");
1801 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1803 TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1804 hframe->SetStats(0);
1806 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1808 hframe->SetXTitle("#mu [GeV]");
1809 hframe->SetYTitle("<E_{loss}> [GeV]");
1812 TGraph *gq=new TGraph(20);
1814 for(Double_t v=0.05;v<=5.05;v+=0.25){
1815 TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1816 Double_t avgloss=dummy->GetMean();
1817 gq->SetPoint(i,v,avgloss);i++;
1820 gq->SetMarkerStyle(20);
1823 TGraph *gg=new TGraph(20);
1825 for(Double_t v=0.05;v<=5.05;v+=0.25){
1826 TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1827 Double_t avgloss=dummy->GetMean();
1828 gg->SetPoint(i,v,avgloss);i++;
1831 gg->SetMarkerStyle(24);
1834 TGraph *gratio=new TGraph(20);
1835 for(Int_t i=0;i<20;i++){
1837 gg->GetPoint(i,x,y);
1838 gq->GetPoint(i,x2,y2);
1840 gratio->SetPoint(i,x,y/y2*10/2.25);
1841 else gratio->SetPoint(i,x,0);
1843 gratio->SetLineStyle(4);
1846 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1847 l1a->SetFillStyle(0);
1848 l1a->SetBorderSize(0);
1850 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1851 l1a->AddEntry(gq,label,"");
1852 l1a->AddEntry(gq,"quark","pl");
1853 l1a->AddEntry(gg,"gluon","pl");
1854 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1860 void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const
1862 // plot average energy loss versus ell
1865 Error("PlotAvgELossVsEll","Tables are not loaded.");
1873 sprintf(title,"Average Energy Loss for Multiple Scattering");
1874 sprintf(name,"cavgelosslms");
1877 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1878 sprintf(name,"cavgelosslsh");
1882 TCanvas *c = new TCanvas(name,title,0,0,600,400);
1884 TH2F *hframe = new TH2F("avglossell",title,2,0,fLengthMax,2,0,250);
1885 hframe->SetStats(0);
1886 hframe->SetXTitle("length [fm]");
1887 hframe->SetYTitle("<E_{loss}> [GeV]");
1890 TGraph *gq=new TGraph((Int_t)fLengthMax*4);
1892 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
1893 TH1F *dummy=ComputeELossHisto(1,medval,v,e);
1894 Double_t avgloss=dummy->GetMean();
1895 gq->SetPoint(i,v,avgloss);i++;
1898 gq->SetMarkerStyle(20);
1901 TGraph *gg=new TGraph((Int_t)fLengthMax*4);
1903 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
1904 TH1F *dummy=ComputeELossHisto(2,medval,v,e);
1905 Double_t avgloss=dummy->GetMean();
1906 gg->SetPoint(i,v,avgloss);i++;
1909 gg->SetMarkerStyle(24);
1912 TGraph *gratio=new TGraph((Int_t)fLengthMax*4);
1913 for(Int_t i=0;i<=(Int_t)fLengthMax*4;i++){
1915 gg->GetPoint(i,x,y);
1916 gq->GetPoint(i,x2,y2);
1918 gratio->SetPoint(i,x,y/y2*100/2.25);
1919 else gratio->SetPoint(i,x,0);
1921 gratio->SetLineStyle(1);
1922 gratio->SetLineWidth(2);
1924 TLegend *l1a = new TLegend(0.15,0.65,.60,0.85);
1925 l1a->SetFillStyle(0);
1926 l1a->SetBorderSize(0);
1929 sprintf(label,"#hat{q} = %.2f GeV^{2}/fm",medval);
1931 sprintf(label,"#mu = %.2f GeV",medval);
1932 l1a->AddEntry(gq,label,"");
1933 l1a->AddEntry(gq,"quark","pl");
1934 l1a->AddEntry(gg,"gluon","pl");
1935 l1a->AddEntry(gratio,"gluon/quark/2.25*100","pl");
1938 TF1 *f=new TF1("f","100",0,fLengthMax);
1945 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
1947 // plot relative energy loss for given
1948 // length and parton energy versus pt
1951 Error("PlotAvgELossVsPt","Tables are not loaded.");
1958 sprintf(title,"Relative Energy Loss for Multiple Scattering");
1959 sprintf(name,"cavgelossvsptms");
1961 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
1962 sprintf(name,"cavgelossvsptsh");
1965 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1967 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
1968 hframe->SetStats(0);
1969 hframe->SetXTitle("p_{T} [GeV]");
1970 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
1973 TGraph *gq=new TGraph(40);
1975 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1976 TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
1977 Double_t avgloss=dummy->GetMean();
1978 gq->SetPoint(i,pt,avgloss/pt);i++;
1981 gq->SetMarkerStyle(20);
1984 TGraph *gg=new TGraph(40);
1986 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1987 TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
1988 Double_t avgloss=dummy->GetMean();
1989 gg->SetPoint(i,pt,avgloss/pt);i++;
1992 gg->SetMarkerStyle(24);
1995 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1996 l1a->SetFillStyle(0);
1997 l1a->SetBorderSize(0);
1999 sprintf(label,"L = %.1f fm",len);
2000 l1a->AddEntry(gq,label,"");
2001 l1a->AddEntry(gq,"quark","pl");
2002 l1a->AddEntry(gg,"gluon","pl");
2008 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
2010 // plot relative energy loss for given
2011 // length distribution and parton energy versus pt
2014 Error("PlotAvgELossVsPt","Tables are not loaded.");
2021 sprintf(title,"Relative Energy Loss for Multiple Scattering");
2022 sprintf(name,"cavgelossvsptms2");
2024 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
2025 sprintf(name,"cavgelossvsptsh2");
2027 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2029 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2030 hframe->SetStats(0);
2031 hframe->SetXTitle("p_{T} [GeV]");
2032 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2035 TGraph *gq=new TGraph(40);
2037 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2038 TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
2039 Double_t avgloss=dummy->GetMean();
2040 gq->SetPoint(i,pt,avgloss/pt);i++;
2043 gq->SetMarkerStyle(20);
2046 TGraph *gg=new TGraph(40);
2048 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2049 TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
2050 Double_t avgloss=dummy->GetMean();
2051 gg->SetPoint(i,pt,avgloss/pt);i++;
2054 gg->SetMarkerStyle(24);
2057 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2058 l1a->SetFillStyle(0);
2059 l1a->SetBorderSize(0);
2061 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
2062 l1a->AddEntry(gq,label,"");
2063 l1a->AddEntry(gq,"quark","pl");
2064 l1a->AddEntry(gg,"gluon","pl");
2070 Int_t AliQuenchingWeights::GetIndex(Double_t len) const
2072 if(len>fLengthMax) len=fLengthMax;
2074 Int_t l=Int_t(len/0.25);
2075 if((len-l*0.5)>0.125) l++;