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>
44 #include "AliQuenchingWeights.h"
46 ClassImp(AliQuenchingWeights)
48 // conversion from fm to GeV^-1: 1 fm = fmGeV GeV^-1
49 const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197;
52 const Double_t AliQuenchingWeights::fgkRMax = 1.e6;
55 const Int_t AliQuenchingWeights::fgkBins = 1300;
56 const Double_t AliQuenchingWeights::fgkMaxBin = 1.3;
58 // counter for histogram labels
59 Int_t AliQuenchingWeights::fgCounter = 0;
62 AliQuenchingWeights::AliQuenchingWeights()
76 fInstanceNumber=fgCounter++;
78 sprintf(name,"hhistoqw_%d",fInstanceNumber);
79 fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin);
80 for(Int_t bin=1;bin<=fgkBins;bin++)
81 fHisto->SetBinContent(bin,0.);
84 AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a)
92 fMultSoft=a.GetMultSoft();;
95 fQTransport=a.GetQTransport();
96 fECMethod=(kECMethod)a.GetECMethod();
97 fLengthMax=a.GetLengthMax();
98 fInstanceNumber=fgCounter++;
100 sprintf(name,"hhistoqw_%d",fInstanceNumber);
101 fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin);
102 for(Int_t bin=1;bin<=fgkBins;bin++)
103 fHisto->SetBinContent(bin,0.);
105 //Missing in the class is the pathname
106 //to the tables, can be added if needed
109 AliQuenchingWeights::~AliQuenchingWeights()
115 void AliQuenchingWeights::Reset()
117 //reset tables if there were used
120 for(Int_t l=0;l<4*fLengthMaxOld;l++){
121 delete fHistos[0][l];
122 delete fHistos[1][l];
129 void AliQuenchingWeights::SetECMethod(kECMethod type)
131 //set energy constraint method
134 if(fECMethod==kDefault)
135 Info("SetECMethod","Energy Constraint Method set to DEFAULT:\nIf (sampled energy loss > parton energy) then sampled energy loss = parton energy.");
136 else if(fECMethod==kReweight)
137 Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy.");
138 else Info("SetECMethod","Energy Constraint Method set to REWEIGHTCONT:\nRequire sampled energy loss <= parton energy (only implemented for FAST method.");
141 Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
143 // read in tables for multiple scattering approximation
144 // path to continuum and to discrete part
146 fTablesLoaded = kFALSE;
150 sprintf(fname,"%s",gSystem->ExpandPathName(contall));
151 //PH ifstream fincont(fname);
152 fstream fincont(fname,ios::in);
153 #if defined(__HP_aCC) || defined(__DECCXX)
154 if(!fincont.rdbuf()->is_open()) return -1;
156 if(!fincont.is_open()) return -1;
160 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
161 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
162 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>fcaq[13][nn]>>
163 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>fcaq[18][nn]>>
164 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>fcaq[23][nn]>>
165 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>fcaq[28][nn]>>
166 fcaq[29][nn]>>fcaq[30][nn]>>fcaq[31][nn]>>fcaq[32][nn]>>fcaq[33][nn])
173 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
174 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
175 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>fcag[13][nn]>>
176 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>fcag[18][nn]>>
177 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>fcag[23][nn]>>
178 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>fcag[28][nn]>>
179 fcag[29][nn]>>fcag[30][nn]>>fcag[31][nn]>>fcag[32][nn]>>fcag[33][nn])
186 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
187 //PH ifstream findisc(fname);
188 fstream findisc(fname,ios::in);
189 #if defined(__HP_aCC) || defined(__DECCXX)
190 if(!findisc.rdbuf()->is_open()) return -1;
192 if(!findisc.is_open()) return -1;
196 while(findisc>>frrr[nn]>>fdaq[nn]) {
201 while(findisc>>frrrg[nn]>>fdag[nn]) {
206 fTablesLoaded = kTRUE;
211 C***************************************************************************
212 C Quenching Weights for Multiple Soft Scattering
217 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
219 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
222 C This package contains quenching weights for gluon radiation in the
223 C multiple soft scattering approximation.
225 C swqmult returns the quenching weight for a quark (ipart=1) or
226 C a gluon (ipart=2) traversing a medium with transport coeficient q and
227 C length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where
228 C wc=0.5*q*L^2 and w is the energy radiated. The output values are
229 C the continuous and discrete (prefactor of the delta function) parts
230 C of the quenching weights.
232 C In order to use this routine, the files cont.all and disc.all need to be
233 C in the working directory.
235 C An initialization of the tables is needed by doing call initmult before
238 C Please, send us any comment:
240 C urs.wiedemann@cern.ch
241 C carlos.salgado@cern.ch
244 C-------------------------------------------------------------------
246 SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete)
248 REAL*8 xx(400), daq(34), caq(34,261), rrr(34)
249 COMMON /dataqua/ xx, daq, caq, rrr
251 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
252 COMMON /dataglu/ xxg, dag, cag, rrrg
254 REAL*8 rrrr,xxxx, continuous, discrete
256 INTEGER nrlow, nrhigh, nxlow, nxhigh
257 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
258 REAL*8 xfraclow, xfrachigh
269 if (rrin.lt.rrr(nr)) then
281 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
282 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
283 if (rrin.gt.10000d0) then
284 rfraclow = dlog(rrhigh/rrin)/dlog(rrhigh/rrlow)
285 rfrachigh = dlog(rrin/rrlow)/dlog(rrhigh/rrlow)
288 if (ipart.eq.1.and.rrin.ge.rrr(1)) then
295 if (ipart.ne.1.and.rrin.ge.rrrg(1)) then
302 if (xxxx.ge.xx(261)) go to 245
304 nxlow = int(xxin/0.01) + 1
306 xfraclow = (xx(nxhigh)-xxin)/0.01
307 xfrachigh = (xxin - xx(nxlow))/0.01
310 clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh)
311 chigh = xfraclow*caq(nrhigh,nxlow)+xfrachigh*caq(nrhigh,nxhigh)
313 clow = xfraclow*cag(nrlow,nxlow)+xfrachigh*cag(nrlow,nxhigh)
314 chigh = xfraclow*cag(nrhigh,nxlow)+xfrachigh*cag(nrhigh,nxhigh)
317 continuous = rfraclow*clow + rfrachigh*chigh
322 discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh)
324 discrete = rfraclow*dag(nrlow) + rfrachigh*dag(nrhigh)
330 REAL*8 xxq(400), daq(34), caq(34,261), rrr(34)
331 COMMON /dataqua/ xxq, daq, caq, rrr
333 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
334 COMMON /dataglu/ xxg, dag, cag, rrrg
336 OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90)
338 read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn),
339 + caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn),
340 + caq(9,nn), caq(10,nn), caq(11,nn), caq(12,nn),
342 + caq(14,nn), caq(15,nn), caq(16,nn), caq(17,nn),
344 + caq(19,nn), caq(20,nn), caq(21,nn), caq(22,nn),
346 + caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn),
348 + caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn),
349 + caq(33,nn), caq(34,nn)
352 read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn),
353 + cag(4,nn), cag(5,nn), cag(6,nn), cag(7,nn), cag(8,nn),
354 + cag(9,nn), cag(10,nn), cag(11,nn), cag(12,nn),
356 + cag(14,nn), cag(15,nn), cag(16,nn), cag(17,nn),
358 + cag(19,nn), cag(20,nn), cag(21,nn), cag(22,nn),
360 + cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn),
362 + cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn),
363 + cag(33,nn), cag(34,nn)
367 OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91)
369 read (21,*) rrr(nn), daq(nn)
372 read (21,*) rrrg(nn), dag(nn)
377 90 PRINT*, 'input - output error'
378 91 PRINT*, 'input - output error #2'
384 =======================================================================
386 Adapted to ROOT macro by A. Dainese - 13/07/2003
387 Ported to class by C. Loizides - 12/02/2004
388 New version for extended R values added - 06/03/2004
391 Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
392 Double_t &continuous,Double_t &discrete) const
394 // Calculate Multiple Scattering approx.
395 // weights for given parton type,
396 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
402 //read-in data before first call
404 Error("CalcMult","Tables are not loaded.");
408 Error("CalcMult","Tables are not loaded for Multiple Scattering.");
412 Double_t rrin = rrrr;
413 Double_t xxin = xxxx;
415 if(xxin>fxx[260]) return -1;
416 Int_t nxlow = (Int_t)(xxin/0.01) + 1;
417 Int_t nxhigh = nxlow + 1;
418 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.01;
419 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01;
422 if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD
423 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
425 Int_t nrlow=0,nrhigh=0;
426 Double_t rrhigh=0,rrlow=0;
427 for(Int_t nr=1; nr<=34; nr++) {
428 if(rrin<frrr[nr-1]) {
431 rrhigh = frrr[nr-1-1];
441 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
442 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
445 rfraclow = TMath::Log2(rrhigh/rrin)/TMath::Log2(rrhigh/rrlow);
446 rfrachigh = TMath::Log2(rrin/rrlow)/TMath::Log2(rrhigh/rrlow);
448 if((ipart==1) && (rrin>=frrr[0]))
455 if((ipart==2) && (rrin>=frrrg[0]))
463 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
465 Double_t clow=0,chigh=0;
467 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
468 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
470 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
471 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
474 continuous = rfraclow*clow + rfrachigh*chigh;
475 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
476 //rfraclow,clow,rfrachigh,chigh,continuous);
479 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
481 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
487 Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *discall)
489 // read in tables for Single Hard Approx.
490 // path to continuum and to discrete part
492 fTablesLoaded = kFALSE;
496 sprintf(fname,"%s",gSystem->ExpandPathName(contall));
497 //PH ifstream fincont(fname);
498 fstream fincont(fname,ios::in);
499 #if defined(__HP_aCC) || defined(__DECCXX)
500 if(!fincont.rdbuf()->is_open()) return -1;
502 if(!fincont.is_open()) return -1;
506 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
507 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
508 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>
510 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>
512 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>
514 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>
523 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
524 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
525 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>
527 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>
529 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>
531 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>
539 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
540 //PH ifstream findisc(fname);
541 fstream findisc(fname,ios::in);
542 #if defined(__HP_aCC) || defined(__DECCXX)
543 if(!findisc.rdbuf()->is_open()) return -1;
545 if(!findisc.is_open()) return -1;
549 while(findisc>>frrr[nn]>>fdaq[nn]) {
554 while(findisc>>frrrg[nn]>>fdag[nn]) {
560 fTablesLoaded = kTRUE;
565 C***************************************************************************
566 C Quenching Weights for Single Hard Scattering
571 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
573 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
576 C This package contains quenching weights for gluon radiation in the
577 C single hard scattering approximation.
579 C swqlin returns the quenching weight for a quark (ipart=1) or
580 C a gluon (ipart=2) traversing a medium with Debye screening mass mu and
581 C length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where
582 C wc=0.5*mu^2*L and w is the energy radiated. The output values are
583 C the continuous and discrete (prefactor of the delta function) parts
584 C of the quenching weights.
586 C In order to use this routine, the files contlin.all and disclin.all
587 C need to be in the working directory.
589 C An initialization of the tables is needed by doing call initlin before
592 C Please, send us any comment:
594 C urs.wiedemann@cern.ch
595 C carlos.salgado@cern.ch
598 C-------------------------------------------------------------------
601 SUBROUTINE swqlin(ipart,rrrr,xxxx,continuous,discrete)
603 REAL*8 xx(400), dalq(30), calq(30,261), rrr(30)
604 COMMON /datalinqua/ xx, dalq, calq, rrr
606 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
607 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
609 REAL*8 rrrr,xxxx, continuous, discrete
611 INTEGER nrlow, nrhigh, nxlow, nxhigh
612 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
613 REAL*8 xfraclow, xfrachigh
619 nxlow = int(xxin/0.038) + 1
621 xfraclow = (xx(nxhigh)-xxin)/0.038
622 xfrachigh = (xxin - xx(nxlow))/0.038
625 if (rrin.lt.rrr(nr)) then
637 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
638 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
641 clow = xfraclow*calq(nrlow,nxlow)+xfrachigh*calq(nrlow,nxhigh)
642 chigh = xfraclow*calq(nrhigh,nxlow)+xfrachigh*calq(nrhigh,nxhigh)
644 clow = xfraclow*calg(nrlow,nxlow)+xfrachigh*calg(nrlow,nxhigh)
645 chigh = xfraclow*calg(nrhigh,nxlow)+xfrachigh*calg(nrhigh,nxhigh)
648 continuous = rfraclow*clow + rfrachigh*chigh
651 discrete = rfraclow*dalq(nrlow) + rfrachigh*dalq(nrhigh)
653 discrete = rfraclow*dalg(nrlow) + rfrachigh*dalg(nrhigh)
659 REAL*8 xxlq(400), dalq(30), calq(30,261), rrr(30)
660 COMMON /datalinqua/ xxlq, dalq, calq, rrr
662 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
663 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
665 OPEN(UNIT=20,FILE='contlin.all',STATUS='OLD',ERR=90)
667 read (20,*) xxlq(nn), calq(1,nn), calq(2,nn), calq(3,nn),
668 + calq(4,nn), calq(5,nn), calq(6,nn), calq(7,nn), calq(8,nn),
669 + calq(9,nn), calq(10,nn), calq(11,nn), calq(12,nn),
671 + calq(14,nn), calq(15,nn), calq(16,nn), calq(17,nn),
673 + calq(19,nn), calq(20,nn), calq(21,nn), calq(22,nn),
675 + calq(24,nn), calq(25,nn), calq(26,nn), calq(27,nn),
677 + calq(29,nn), calq(30,nn)
680 read (20,*) xxlg(nn), calg(1,nn), calg(2,nn), calg(3,nn),
681 + calg(4,nn), calg(5,nn), calg(6,nn), calg(7,nn), calg(8,nn),
682 + calg(9,nn), calg(10,nn), calg(11,nn), calg(12,nn),
684 + calg(14,nn), calg(15,nn), calg(16,nn), calg(17,nn),
686 + calg(19,nn), calg(20,nn), calg(21,nn), calg(22,nn),
688 + calg(24,nn), calg(25,nn), calg(26,nn), calg(27,nn),
690 + calg(29,nn), calg(30,nn)
694 OPEN(UNIT=21,FILE='disclin.all',STATUS='OLD',ERR=91)
696 read (21,*) rrr(nn), dalq(nn)
699 read (21,*) rrrlg(nn), dalg(nn)
704 90 PRINT*, 'input - output error'
705 91 PRINT*, 'input - output error #2'
710 =======================================================================
712 Ported to class by C. Loizides - 17/02/2004
716 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx,
717 Double_t &continuous,Double_t &discrete) const
719 // calculate Single Hard approx.
720 // weights for given parton type,
721 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
723 // read-in data before first call
725 Error("CalcSingleHard","Tables are not loaded.");
729 Error("CalcSingleHard","Tables are not loaded for Single Hard Scattering.");
733 Double_t rrin = rrrr;
734 Double_t xxin = xxxx;
736 Int_t nxlow = (Int_t)(xxin/0.038) + 1;
737 Int_t nxhigh = nxlow + 1;
738 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.038;
739 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.038;
742 if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD
743 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
745 Int_t nrlow=0,nrhigh=0;
746 Double_t rrhigh=0,rrlow=0;
747 for(Int_t nr=1; nr<=30; nr++) {
748 if(rrin<frrr[nr-1]) {
751 rrhigh = frrr[nr-1-1];
761 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
762 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
764 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
766 Double_t clow=0,chigh=0;
768 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
769 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
771 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
772 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
775 continuous = rfraclow*clow + rfrachigh*chigh;
776 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
777 // rfraclow,clow,rfrachigh,chigh,continuous);
780 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
782 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
788 Int_t AliQuenchingWeights::CalcMult(Int_t ipart,
789 Double_t w,Double_t qtransp,Double_t length,
790 Double_t &continuous,Double_t &discrete) const
792 // Calculate Multiple Scattering approx.
793 // weights for given parton type,
794 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
796 Double_t wc=CalcWC(qtransp,length);
797 Double_t rrrr=CalcR(wc,length);
799 return CalcMult(ipart,rrrr,xxxx,continuous,discrete);
802 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
803 Double_t w,Double_t mu,Double_t length,
804 Double_t &continuous,Double_t &discrete) const
806 // calculate Single Hard approx.
807 // weights for given parton type,
808 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
810 Double_t wcbar=CalcWCbar(mu,length);
811 Double_t rrrr=CalcR(wcbar,length);
812 Double_t xxxx=w/wcbar;
813 return CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
816 Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const
818 //calculate r value and
819 //check if it is less then maximum
821 Double_t r = wc*l*fgkConvFmToInvGeV;
823 Warning("CalcR","Value of r = %.2f; should be less than %.2f", r, fgkRMax);
829 Double_t AliQuenchingWeights::CalcRk(Double_t k, Double_t I0, Double_t I1) const
831 //calculate R value and
832 //check if it is less then maximum
834 Double_t r = fgkRMax-1;
838 Warning("CalcRk","Value of r = %.2f; should be less than %.2f",r,fgkRMax);
844 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const
846 // return DeltaE for MS or SH scattering
847 // for given parton type, length and energy
848 // Dependant on ECM (energy constraint method)
849 // e is used to determine where to set bins to zero.
852 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
855 if((ipart<1) || (ipart>2)) {
856 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
860 Int_t l=GetIndex(length);
863 if(fECMethod==kReweight){
867 ret=fHistos[ipart-1][l-1]->GetRandom();
869 Warning("GetELossRandom",
870 "Aborted reweighting; maximum loss assigned after 1e6 trials.");
877 Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
882 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const
884 //return quenched parton energy
885 //for given parton type, length and energy
887 Double_t loss=GetELossRandom(ipart,length,e);
891 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const
893 // return DeltaE for MS or SH scattering
894 // for given parton type, length distribution and energy
895 // Dependant on ECM (energy constraint method)
896 // e is used to determine where to set bins to zero.
899 Warning("GetELossRandom","Pointer to length distribution is NULL.");
902 Double_t ell=hell->GetRandom();
903 return GetELossRandom(ipart,ell,e);
906 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e) const
908 //return quenched parton energy
909 //for given parton type, length distribution and energy
911 Double_t loss=GetELossRandom(ipart,hell,e);
915 Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
917 // return DeltaE for new dynamic version
918 // for given parton type, I0 and I1 value and energy
919 // Dependant on ECM (energy constraint method)
920 // e is used to determine where to set bins to zero.
922 // read-in data before first call
924 Fatal("GetELossRandomK","Tables are not loaded.");
927 if((ipart<1) || (ipart>2)) {
928 Fatal("GetELossRandomK","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
932 Double_t r=CalcRk(I0,I1);
934 Fatal("GetELossRandomK","R should not be negative");
937 Double_t wc=CalcWCk(I1);
939 Fatal("GetELossRandomK","wc should be greater than zero");
942 if(SampleEnergyLoss(ipart,r)!=0){
943 Fatal("GetELossRandomK","Could not sample energy loss");
947 if(fECMethod==kReweight){
951 ret=fHisto->GetRandom();
953 Warning("GetELossRandomK",
954 "Aborted reweighting; maximum loss assigned after 1e6 trials.");
962 Double_t ret=fHisto->GetRandom()*wc;
967 Double_t AliQuenchingWeights::CalcQuenchedEnergyK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
969 //return quenched parton energy
970 //for given parton type, I0 and I1 value and energy
972 Double_t loss=GetELossRandomK(ipart,I0,I1,e);
976 Double_t AliQuenchingWeights::GetELossRandomKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
978 // return DeltaE for new dynamic version
979 // for given parton type, I0 and I1 value and energy
980 // Dependant on ECM (energy constraint method)
981 // e is used to determine where to set bins to zero.
982 // method is optimized and should only be used if
983 // all parameters are well within the bounds.
984 // read-in data tables before first call
986 Double_t r=CalcRk(I0,I1);
991 Double_t wc=CalcWCk(I1);
996 return GetELossRandomKFastR(ipart,r,wc,e);
999 Double_t AliQuenchingWeights::GetELossRandomKFastR(Int_t ipart, Double_t r, Double_t wc, Double_t e)
1001 // return DeltaE for new dynamic version
1002 // for given parton type, R and wc value and energy
1003 // Dependant on ECM (energy constraint method)
1004 // e is used to determine where to set bins to zero.
1005 // method is optimized and should only be used if
1006 // all parameters are well within the bounds.
1007 // read-in data tables before first call
1013 Double_t discrete=0.;
1014 Double_t continuous=0.;
1016 Double_t xxxx = fHisto->GetBinCenter(bin);
1018 CalcMult(ipart,r,xxxx,continuous,discrete);
1020 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1023 return 0.; //no energy loss
1026 fHisto->SetBinContent(bin,continuous);
1027 Int_t kbinmax=fHisto->FindBin(e/wc);
1028 if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
1029 if(kbinmax==1) return e; //maximum energy loss
1032 for(Int_t bin=2; bin<=kbinmax; bin++) {
1033 xxxx = fHisto->GetBinCenter(bin);
1034 CalcMult(ipart,r,xxxx,continuous,discrete);
1035 fHisto->SetBinContent(bin,continuous);
1038 for(Int_t bin=2; bin<=kbinmax; bin++) {
1039 xxxx = fHisto->GetBinCenter(bin);
1040 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1041 fHisto->SetBinContent(bin,continuous);
1045 if(fECMethod==kReweight){
1046 fHisto->SetBinContent(kbinmax+1,0);
1047 fHisto->Fill(0.,discrete*fgkBins/fgkMaxBin);
1048 } else if (fECMethod==kReweightCont) {
1049 fHisto->SetBinContent(kbinmax+1,0);
1050 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1051 fHisto->Scale(1./kdelta*(1-discrete));
1052 fHisto->Fill(0.,discrete);
1054 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1055 Double_t val=discrete*fgkBins/fgkMaxBin;
1056 fHisto->Fill(0.,val);
1057 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1059 for(Int_t bin=kbinmax+2; bin<=fgkBins; bin++) {
1060 fHisto->SetBinContent(bin,0);
1062 //cout << kbinmax << " " << discrete << " " << fHisto->Integral() << endl;
1063 Double_t ret=fHisto->GetRandom()*wc;
1068 Double_t AliQuenchingWeights::CalcQuenchedEnergyKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
1070 //return quenched parton energy (fast method)
1071 //for given parton type, I0 and I1 value and energy
1073 Double_t loss=GetELossRandomKFast(ipart,I0,I1,e);
1077 Double_t AliQuenchingWeights::GetDiscreteWeight(Int_t ipart, Double_t I0, Double_t I1)
1079 // return discrete weight
1081 Double_t r=CalcRk(I0,I1);
1085 return GetDiscreteWeightR(ipart,r);
1088 Double_t AliQuenchingWeights::GetDiscreteWeightR(Int_t ipart, Double_t r)
1090 // return discrete weight
1096 Double_t discrete=0.;
1097 Double_t continuous=0.;
1099 Double_t xxxx = fHisto->GetBinCenter(bin);
1101 CalcMult(ipart,r,xxxx,continuous,discrete);
1103 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1107 void AliQuenchingWeights::GetZeroLossProb(Double_t &p,Double_t &prw,Double_t &prwcont,
1108 Int_t ipart,Double_t I0,Double_t I1,Double_t e)
1110 //calculate the probabilty that there is no energy
1111 //loss for different ways of energy constraint
1112 p=1.;prw=1.;prwcont=1.;
1113 Double_t r=CalcRk(I0,I1);
1117 Double_t wc=CalcWCk(I1);
1121 GetZeroLossProbR(p,prw,prwcont,ipart,r,wc,e);
1124 void AliQuenchingWeights::GetZeroLossProbR(Double_t &p,Double_t &prw,Double_t &prwcont,
1125 Int_t ipart, Double_t r,Double_t wc,Double_t e)
1127 //calculate the probabilty that there is no energy
1128 //loss for different ways of energy constraint
1133 Double_t discrete=0.;
1134 Double_t continuous=0.;
1136 Int_t kbinmax=fHisto->FindBin(e/wc);
1137 if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
1139 for(Int_t bin=1; bin<=kbinmax; bin++) {
1140 Double_t xxxx = fHisto->GetBinCenter(bin);
1141 CalcMult(ipart,r,xxxx,continuous,discrete);
1142 fHisto->SetBinContent(bin,continuous);
1145 for(Int_t bin=1; bin<=kbinmax; bin++) {
1146 Double_t xxxx = fHisto->GetBinCenter(bin);
1147 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1148 fHisto->SetBinContent(bin,continuous);
1152 //non-reweighted P(Delta E = 0)
1153 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1154 Double_t val=discrete*fgkBins/fgkMaxBin;
1155 fHisto->Fill(0.,val);
1156 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1157 Double_t hint=fHisto->Integral(1,kbinmax+1);
1158 p=fHisto->GetBinContent(1)/hint;
1161 hint=fHisto->Integral(1,kbinmax);
1162 prw=fHisto->GetBinContent(1)/hint;
1164 Double_t xxxx = fHisto->GetBinCenter(1);
1165 CalcMult(ipart,r,xxxx,continuous,discrete);
1166 fHisto->SetBinContent(1,continuous);
1167 hint=fHisto->Integral(1,kbinmax);
1168 fHisto->Scale(1./hint*(1-discrete));
1169 fHisto->Fill(0.,discrete);
1170 prwcont=fHisto->GetBinContent(1);
1173 Int_t AliQuenchingWeights::SampleEnergyLoss()
1175 // Has to be called to fill the histograms
1177 // For stored values fQTransport loop over
1178 // particle type and length = 1 to fMaxLength (fm)
1179 // to fill energy loss histos
1181 // Take histogram of continuous weights
1182 // Take discrete_weight
1183 // If discrete_weight > 1, put all channels to 0, except channel 1
1184 // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral
1186 // read-in data before first call
1188 Error("SampleEnergyLoss","Tables are not loaded.");
1193 Int_t lmax=CalcLengthMax(fQTransport);
1194 if(fLengthMax>lmax){
1195 Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
1199 Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested.");
1203 fHistos=new TH1F**[2];
1204 fHistos[0]=new TH1F*[4*fLengthMax];
1205 fHistos[1]=new TH1F*[4*fLengthMax];
1206 fLengthMaxOld=fLengthMax; //remember old value in case
1207 //user wants to reset
1210 Char_t meddesc[100];
1212 medvalue=(Int_t)(fQTransport*1000.);
1213 sprintf(meddesc,"MS");
1215 medvalue=(Int_t)(fMu*1000.);
1216 sprintf(meddesc,"SH");
1219 for(Int_t ipart=1;ipart<=2;ipart++){
1220 for(Int_t l=1;l<=4*fLengthMax;l++){
1222 sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
1224 Double_t wc = CalcWC(len);
1225 fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgkBins,0.,fgkMaxBin*wc);
1226 fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
1227 fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
1228 fHistos[ipart-1][l-1]->SetLineColor(4);
1230 Double_t rrrr = CalcR(wc,len);
1231 Double_t discrete=0.;
1232 // loop on histogram channels
1233 for(Int_t bin=1; bin<=fgkBins; bin++) {
1234 Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
1235 Double_t continuous;
1237 CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1239 CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1240 fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
1242 // add discrete part to distribution
1244 for(Int_t bin=2;bin<=fgkBins;bin++)
1245 fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
1247 Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1248 fHistos[ipart-1][l-1]->Fill(0.,val);
1250 Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1251 fHistos[ipart-1][l-1]->Scale(1./hint);
1257 Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t r)
1259 // Sample energy loss directly for one particle type
1260 // choses R (safe it and keep it until next call of function)
1262 // read-in data before first call
1264 Error("SampleEnergyLoss","Tables are not loaded.");
1268 Double_t discrete=0.;
1269 Double_t continuous=0;;
1271 Double_t xxxx = fHisto->GetBinCenter(bin);
1273 CalcMult(ipart,r,xxxx,continuous,discrete);
1275 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1278 fHisto->SetBinContent(1,1.);
1279 for(Int_t bin=2;bin<=fgkBins;bin++)
1280 fHisto->SetBinContent(bin,0.);
1284 fHisto->SetBinContent(bin,continuous);
1285 for(Int_t bin=2; bin<=fgkBins; bin++) {
1286 xxxx = fHisto->GetBinCenter(bin);
1288 CalcMult(ipart,r,xxxx,continuous,discrete);
1290 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1291 fHisto->SetBinContent(bin,continuous);
1294 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
1295 fHisto->Fill(0.,val);
1296 Double_t hint=fHisto->Integral(1,fgkBins);
1298 fHisto->Scale(1./hint);
1300 //cout << discrete << " " << hint << " " << continuous << endl;
1306 const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) const
1308 //return quenching histograms
1309 //for ipart and length
1312 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
1315 if((ipart<1) || (ipart>2)) {
1316 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
1320 Int_t l=GetIndex(length);
1322 return fHistos[ipart-1][l-1];
1325 TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const
1327 // ipart = 1 for quark, 2 for gluon
1328 // medval a) qtransp = transport coefficient (GeV^2/fm)
1329 // b) mu = Debye mass (GeV)
1330 // length = path length in medium (fm)
1331 // Get from SW tables:
1332 // - continuous weight, as a function of dE/wc
1335 Char_t meddesc[100];
1337 wc=CalcWC(medval,length);
1338 sprintf(meddesc,"MS");
1340 wc=CalcWCbar(medval,length);
1341 sprintf(meddesc,"SH");
1345 sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart,
1346 (Int_t)(medval*1000.),(Int_t)length);
1348 TH1F *hist = new TH1F("hist",hname,fgkBins,0.,fgkMaxBin*wc);
1349 hist->SetXTitle("#Delta E [GeV]");
1350 hist->SetYTitle("p(#Delta E)");
1351 hist->SetLineColor(4);
1353 Double_t rrrr = CalcR(wc,length);
1354 //loop on histogram channels
1355 for(Int_t bin=1; bin<=fgkBins; bin++) {
1356 Double_t xxxx = hist->GetBinCenter(bin)/wc;
1357 Double_t continuous,discrete;
1359 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1360 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1365 hist->SetBinContent(bin,continuous);
1370 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const
1372 // ipart = 1 for quark, 2 for gluon
1373 // medval a) qtransp = transport coefficient (GeV^2/fm)
1374 // b) mu = Debye mass (GeV)
1375 // length = path length in medium (fm)
1376 // Get from SW tables:
1377 // - continuous weight, as a function of dE/wc
1380 Char_t meddesc[100];
1382 wc=CalcWC(medval,length);
1383 sprintf(meddesc,"MS");
1385 wc=CalcWCbar(medval,length);
1386 sprintf(meddesc,"SH");
1390 sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart,
1391 (Int_t)(medval*1000.),(Int_t)length);
1393 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
1394 histx->SetXTitle("x = #Delta E/#omega_{c}");
1396 histx->SetYTitle("p(#Delta E/#omega_{c})");
1398 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1399 histx->SetLineColor(4);
1401 Double_t rrrr = CalcR(wc,length);
1402 //loop on histogram channels
1403 for(Int_t bin=1; bin<=fgkBins; bin++) {
1404 Double_t xxxx = histx->GetBinCenter(bin);
1405 Double_t continuous,discrete;
1407 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1408 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1413 histx->SetBinContent(bin,continuous);
1418 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t r) const
1420 // compute P(E) distribution for
1421 // given ipart = 1 for quark, 2 for gluon
1424 Char_t meddesc[100];
1426 sprintf(meddesc,"MS");
1428 sprintf(meddesc,"SH");
1432 sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,r);
1433 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
1434 histx->SetXTitle("x = #Delta E/#omega_{c}");
1436 histx->SetYTitle("p(#Delta E/#omega_{c})");
1438 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1439 histx->SetLineColor(4);
1442 Double_t continuous=0.,discrete=0.;
1443 //loop on histogram channels
1444 for(Int_t bin=1; bin<=fgkBins; bin++) {
1445 Double_t xxxx = histx->GetBinCenter(bin);
1447 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1448 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1453 histx->SetBinContent(bin,continuous);
1456 //add discrete part to distribution
1458 for(Int_t bin=2;bin<=fgkBins;bin++)
1459 histx->SetBinContent(bin,0.);
1461 Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins);
1462 histx->Fill(0.,val);
1464 Double_t hint=histx->Integral(1,fgkBins);
1465 if(hint!=0) histx->Scale(1./hint);
1470 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
1472 // compute energy loss histogram for
1473 // parton type, medium value, length and energy
1475 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1477 dummy->SetQTransport(medval);
1480 dummy->SetMu(medval);
1481 dummy->InitSingleHard();
1483 dummy->SampleEnergyLoss();
1488 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1489 sprintf(hname,"hLossQuarks");
1491 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1492 sprintf(hname,"hLossGluons");
1495 TH1F *h = new TH1F(hname,name,250,0,250);
1496 for(Int_t i=0;i<100000;i++){
1497 //if(i % 1000 == 0) cout << "." << flush;
1498 Double_t loss=dummy->GetELossRandom(ipart,l,e);
1506 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
1508 // compute energy loss histogram for
1509 // parton type, medium value,
1510 // length distribution and energy
1512 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1514 dummy->SetQTransport(medval);
1517 dummy->SetMu(medval);
1518 dummy->InitSingleHard();
1520 dummy->SampleEnergyLoss();
1525 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1526 sprintf(hname,"hLossQuarks");
1528 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1529 sprintf(hname,"hLossGluons");
1532 TH1F *h = new TH1F(hname,name,250,0,250);
1533 for(Int_t i=0;i<100000;i++){
1534 //if(i % 1000 == 0) cout << "." << flush;
1535 Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
1543 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t r) const
1545 // compute energy loss histogram for
1546 // parton type and given R
1548 TH1F *dummy = ComputeQWHistoX(ipart,r);
1549 if(!dummy) return 0;
1552 sprintf(hname,"hELossHistox_%d_%.2f",ipart,r);
1553 TH1F *histx = new TH1F("histxr",hname,fgkBins,0.,fgkMaxBin);
1554 for(Int_t i=0;i<100000;i++){
1555 //if(i % 1000 == 0) cout << "." << flush;
1556 Double_t loss=dummy->GetRandom();
1563 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,Double_t l) const
1565 // compute average energy loss for
1566 // parton type, medium value, length and energy
1568 TH1F *dummy = ComputeELossHisto(ipart,medval,l);
1569 if(!dummy) return 0;
1570 Double_t ret=dummy->GetMean();
1575 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEll) const
1577 // compute average energy loss for
1578 // parton type, medium value,
1579 // length distribution and energy
1581 TH1F *dummy = ComputeELossHisto(ipart,medval,hEll);
1582 if(!dummy) return 0;
1583 Double_t ret=dummy->GetMean();
1588 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t r) const
1590 // compute average energy loss over wc
1591 // for parton type and given R
1593 TH1F *dummy = ComputeELossHisto(ipart,r);
1594 if(!dummy) return 0;
1595 Double_t ret=dummy->GetMean();
1600 void AliQuenchingWeights::PlotDiscreteWeights(Double_t len,Double_t qm) const
1602 // plot discrete weights for given length
1606 c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
1608 c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
1611 TH2F *hframe = new TH2F("hdisc","",2,0,qm+.1,2,0,1.25);
1612 hframe->SetStats(0);
1614 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1616 hframe->SetXTitle("#mu [GeV]");
1617 //hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
1618 hframe->SetYTitle("p_{0} (discrete weight)");
1621 Int_t points=(Int_t)qm*4;
1622 TGraph *gq=new TGraph(points);
1625 for(Double_t q=0.05;q<=qm+.05;q+=0.25){
1627 CalcMult(1,1.0,q,len,cont,disc);
1628 gq->SetPoint(i,q,disc);i++;
1631 for(Double_t m=0.05;m<=qm+.05;m+=0.25){
1633 CalcSingleHard(1,1.0,m,len,cont, disc);
1634 gq->SetPoint(i,m,disc);i++;
1637 gq->SetMarkerStyle(20);
1638 gq->SetMarkerColor(1);
1639 gq->SetLineStyle(1);
1640 gq->SetLineColor(1);
1643 TGraph *gg=new TGraph(points);
1646 for(Double_t q=0.05;q<=qm+.05;q+=0.25){
1648 CalcMult(2,1.0,q,len,cont,disc);
1649 gg->SetPoint(i,q,disc);i++;
1652 for(Double_t m=0.05;m<=qm+.05;m+=0.25){
1654 CalcSingleHard(2,1.0,m,len,cont,disc);
1655 gg->SetPoint(i,m,disc);i++;
1658 gg->SetMarkerStyle(24);
1659 gg->SetMarkerColor(2);
1660 gg->SetLineStyle(2);
1661 gg->SetLineColor(2);
1664 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1665 l1a->SetFillStyle(0);
1666 l1a->SetBorderSize(0);
1668 sprintf(label,"L = %.1f fm",len);
1669 l1a->AddEntry(gq,label,"");
1670 l1a->AddEntry(gq,"quark projectile","l");
1671 l1a->AddEntry(gg,"gluon projectile","l");
1677 void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
1679 // plot continous weights for
1680 // given parton type and length
1687 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1689 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1691 medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1692 sprintf(name,"ccont-ms-%d",itype);
1695 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1697 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1699 medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1700 sprintf(name,"ccont-ms-%d",itype);
1703 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1705 TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1707 h1->SetTitle(title);
1709 h1->SetLineColor(1);
1711 TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1713 h2->SetLineColor(2);
1714 h2->DrawCopy("SAME");
1715 TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1717 h3->SetLineColor(3);
1718 h3->DrawCopy("SAME");
1720 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1721 l1a->SetFillStyle(0);
1722 l1a->SetBorderSize(0);
1724 sprintf(label,"L = %.1f fm",ell);
1725 l1a->AddEntry(h1,label,"");
1727 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1728 l1a->AddEntry(h1,label,"pl");
1729 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1730 l1a->AddEntry(h2,label,"pl");
1731 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1732 l1a->AddEntry(h3,label,"pl");
1734 sprintf(label,"#mu = %.1f GeV",medvals[0]);
1735 l1a->AddEntry(h1,label,"pl");
1736 sprintf(label,"#mu = %.1f GeV",medvals[1]);
1737 l1a->AddEntry(h2,label,"pl");
1738 sprintf(label,"#mu = %.1f GeV",medvals[2]);
1739 l1a->AddEntry(h3,label,"pl");
1746 void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
1748 // plot continous weights for
1749 // given parton type and medium value
1755 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1757 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1759 sprintf(name,"ccont2-ms-%d",itype);
1762 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1764 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1766 sprintf(name,"ccont2-sh-%d",itype);
1768 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1770 TH1F *h1=ComputeQWHisto(itype,medval,8);
1772 h1->SetTitle(title);
1774 h1->SetLineColor(1);
1776 TH1F *h2=ComputeQWHisto(itype,medval,5);
1778 h2->SetLineColor(2);
1779 h2->DrawCopy("SAME");
1780 TH1F *h3=ComputeQWHisto(itype,medval,2);
1782 h3->SetLineColor(3);
1783 h3->DrawCopy("SAME");
1785 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1786 l1a->SetFillStyle(0);
1787 l1a->SetBorderSize(0);
1790 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval);
1792 sprintf(label,"#mu = %.1f GeV",medval);
1794 l1a->AddEntry(h1,label,"");
1795 l1a->AddEntry(h1,"L = 8 fm","pl");
1796 l1a->AddEntry(h2,"L = 5 fm","pl");
1797 l1a->AddEntry(h3,"L = 2 fm","pl");
1803 void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t qm,Double_t e) const
1805 // plot average energy loss for given length
1806 // and parton energy
1809 Error("PlotAvgELoss","Tables are not loaded.");
1816 sprintf(title,"Average Energy Loss for Multiple Scattering");
1817 sprintf(name,"cavgelossms");
1819 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1820 sprintf(name,"cavgelosssh");
1823 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1825 TH2F *hframe = new TH2F("avgloss","",2,0,qm+.1,2,0,100);
1826 hframe->SetStats(0);
1828 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1830 hframe->SetXTitle("#mu [GeV]");
1831 hframe->SetYTitle("<E_{loss}> [GeV]");
1834 TGraph *gq=new TGraph(20);
1836 for(Double_t v=0.05;v<=qm+.05;v+=0.25){
1837 TH1F *dummy=ComputeELossHisto(1,v,len,e);
1838 Double_t avgloss=dummy->GetMean();
1839 gq->SetPoint(i,v,avgloss);i++;
1842 gq->SetMarkerStyle(21);
1845 Int_t points=(Int_t)qm*4;
1846 TGraph *gg=new TGraph(points);
1848 for(Double_t v=0.05;v<=qm+.05;v+=0.25){
1849 TH1F *dummy=ComputeELossHisto(2,v,len,e);
1850 Double_t avgloss=dummy->GetMean();
1851 gg->SetPoint(i,v,avgloss);i++;
1854 gg->SetMarkerStyle(20);
1855 gg->SetMarkerColor(2);
1858 TGraph *gratio=new TGraph(points);
1859 for(Int_t i=0;i<points;i++){
1861 gg->GetPoint(i,x,y);
1862 gq->GetPoint(i,x2,y2);
1864 gratio->SetPoint(i,x,y/y2*10/2.25);
1865 else gratio->SetPoint(i,x,0);
1867 gratio->SetLineStyle(4);
1869 TLegend *l1a = new TLegend(0.15,0.60,0.50,0.90);
1870 l1a->SetFillStyle(0);
1871 l1a->SetBorderSize(0);
1873 sprintf(label,"L = %.1f fm",len);
1874 l1a->AddEntry(gq,label,"");
1875 l1a->AddEntry(gq,"quark projectile","pl");
1876 l1a->AddEntry(gg,"gluon projectile","pl");
1877 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1883 void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
1885 // plot average energy loss for given
1886 // length distribution and parton energy
1889 Error("PlotAvgELossVs","Tables are not loaded.");
1896 sprintf(title,"Average Energy Loss for Multiple Scattering");
1897 sprintf(name,"cavgelossms2");
1899 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1900 sprintf(name,"cavgelosssh2");
1903 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1905 TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1906 hframe->SetStats(0);
1908 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1910 hframe->SetXTitle("#mu [GeV]");
1911 hframe->SetYTitle("<E_{loss}> [GeV]");
1914 TGraph *gq=new TGraph(20);
1916 for(Double_t v=0.05;v<=5.05;v+=0.25){
1917 TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1918 Double_t avgloss=dummy->GetMean();
1919 gq->SetPoint(i,v,avgloss);i++;
1922 gq->SetMarkerStyle(20);
1925 TGraph *gg=new TGraph(20);
1927 for(Double_t v=0.05;v<=5.05;v+=0.25){
1928 TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1929 Double_t avgloss=dummy->GetMean();
1930 gg->SetPoint(i,v,avgloss);i++;
1933 gg->SetMarkerStyle(24);
1936 TGraph *gratio=new TGraph(20);
1937 for(Int_t i=0;i<20;i++){
1939 gg->GetPoint(i,x,y);
1940 gq->GetPoint(i,x2,y2);
1942 gratio->SetPoint(i,x,y/y2*10/2.25);
1943 else gratio->SetPoint(i,x,0);
1945 gratio->SetLineStyle(4);
1948 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1949 l1a->SetFillStyle(0);
1950 l1a->SetBorderSize(0);
1952 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1953 l1a->AddEntry(gq,label,"");
1954 l1a->AddEntry(gq,"quark","pl");
1955 l1a->AddEntry(gg,"gluon","pl");
1956 //l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1962 void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const
1964 // plot average energy loss versus ell
1967 Error("PlotAvgELossVsEll","Tables are not loaded.");
1975 sprintf(title,"Average Energy Loss for Multiple Scattering");
1976 sprintf(name,"cavgelosslms");
1979 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1980 sprintf(name,"cavgelosslsh");
1984 TCanvas *c = new TCanvas(name,title,0,0,600,400);
1986 TH2F *hframe = new TH2F("avglossell",title,2,0,fLengthMax,2,0,250);
1987 hframe->SetStats(0);
1988 hframe->SetXTitle("length [fm]");
1989 hframe->SetYTitle("<E_{loss}> [GeV]");
1992 TGraph *gq=new TGraph((Int_t)fLengthMax*4);
1994 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
1995 TH1F *dummy=ComputeELossHisto(1,medval,v,e);
1996 Double_t avgloss=dummy->GetMean();
1997 gq->SetPoint(i,v,avgloss);i++;
2000 gq->SetMarkerStyle(20);
2003 TGraph *gg=new TGraph((Int_t)fLengthMax*4);
2005 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
2006 TH1F *dummy=ComputeELossHisto(2,medval,v,e);
2007 Double_t avgloss=dummy->GetMean();
2008 gg->SetPoint(i,v,avgloss);i++;
2011 gg->SetMarkerStyle(24);
2014 TGraph *gratio=new TGraph((Int_t)fLengthMax*4);
2015 for(Int_t i=0;i<=(Int_t)fLengthMax*4;i++){
2017 gg->GetPoint(i,x,y);
2018 gq->GetPoint(i,x2,y2);
2020 gratio->SetPoint(i,x,y/y2*100/2.25);
2021 else gratio->SetPoint(i,x,0);
2023 gratio->SetLineStyle(1);
2024 gratio->SetLineWidth(2);
2026 TLegend *l1a = new TLegend(0.15,0.65,.60,0.85);
2027 l1a->SetFillStyle(0);
2028 l1a->SetBorderSize(0);
2031 sprintf(label,"#hat{q} = %.2f GeV^{2}/fm",medval);
2033 sprintf(label,"#mu = %.2f GeV",medval);
2034 l1a->AddEntry(gq,label,"");
2035 l1a->AddEntry(gq,"quark","pl");
2036 l1a->AddEntry(gg,"gluon","pl");
2037 l1a->AddEntry(gratio,"gluon/quark/2.25*100","pl");
2040 TF1 *f=new TF1("f","100",0,fLengthMax);
2047 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
2049 // plot relative energy loss for given
2050 // length and parton energy versus pt
2053 Error("PlotAvgELossVsPt","Tables are not loaded.");
2060 sprintf(title,"Relative Energy Loss for Multiple Scattering");
2061 sprintf(name,"cavgelossvsptms");
2063 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
2064 sprintf(name,"cavgelossvsptsh");
2067 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2069 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2070 hframe->SetStats(0);
2071 hframe->SetXTitle("p_{T} [GeV]");
2072 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2075 TGraph *gq=new TGraph(40);
2077 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2078 TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
2079 Double_t avgloss=dummy->GetMean();
2080 gq->SetPoint(i,pt,avgloss/pt);i++;
2083 gq->SetMarkerStyle(20);
2086 TGraph *gg=new TGraph(40);
2088 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2089 TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
2090 Double_t avgloss=dummy->GetMean();
2091 gg->SetPoint(i,pt,avgloss/pt);i++;
2094 gg->SetMarkerStyle(24);
2097 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2098 l1a->SetFillStyle(0);
2099 l1a->SetBorderSize(0);
2101 sprintf(label,"L = %.1f fm",len);
2102 l1a->AddEntry(gq,label,"");
2103 l1a->AddEntry(gq,"quark","pl");
2104 l1a->AddEntry(gg,"gluon","pl");
2110 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
2112 // plot relative energy loss for given
2113 // length distribution and parton energy versus pt
2116 Error("PlotAvgELossVsPt","Tables are not loaded.");
2123 sprintf(title,"Relative Energy Loss for Multiple Scattering");
2124 sprintf(name,"cavgelossvsptms2");
2126 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
2127 sprintf(name,"cavgelossvsptsh2");
2129 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2131 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2132 hframe->SetStats(0);
2133 hframe->SetXTitle("p_{T} [GeV]");
2134 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2137 TGraph *gq=new TGraph(40);
2139 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2140 TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
2141 Double_t avgloss=dummy->GetMean();
2142 gq->SetPoint(i,pt,avgloss/pt);i++;
2145 gq->SetMarkerStyle(20);
2148 TGraph *gg=new TGraph(40);
2150 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2151 TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
2152 Double_t avgloss=dummy->GetMean();
2153 gg->SetPoint(i,pt,avgloss/pt);i++;
2156 gg->SetMarkerStyle(24);
2159 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2160 l1a->SetFillStyle(0);
2161 l1a->SetBorderSize(0);
2163 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
2164 l1a->AddEntry(gq,label,"");
2165 l1a->AddEntry(gq,"quark","pl");
2166 l1a->AddEntry(gg,"gluon","pl");
2172 Int_t AliQuenchingWeights::GetIndex(Double_t len) const
2174 //get the index according to length
2175 if(len>fLengthMax) len=fLengthMax;
2177 Int_t l=Int_t(len/0.25);
2178 if((len-l*0.25)>0.125) l++;