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::fgBins = 1300;
57 const Double_t AliQuenchingWeights::fgMaxBin = 1.3;
59 // counter for histogram labels
60 Int_t AliQuenchingWeights::fgCounter = 0;
63 AliQuenchingWeights::AliQuenchingWeights()
74 fECMethod=kReweight; //this is to force printout
78 fInstanceNumber=fgCounter++;
80 sprintf(name,"hhistoqw_%d",fInstanceNumber);
81 fHisto = new TH1F(name,"",fgBins,0.,fgMaxBin);
82 for(Int_t bin=1;bin<=fgBins;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,"",fgBins,0.,fgMaxBin);
104 for(Int_t bin=1;bin<=fgBins;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
135 if(fECMethod==type) return;
137 if(fECMethod==kDefault)
138 Info("SetECMethod","Energy Constraint Method set to DEFAULT:\nIf (sampled energy loss > parton energy) then sampled energy loss = parton energy.");
140 Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy.");
143 Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
145 // read in tables for multiple scattering approximation
146 // path to continuum and to discrete part
148 fTablesLoaded = kFALSE;
152 sprintf(fname,"%s",gSystem->ExpandPathName(contall));
153 //PH ifstream fincont(fname);
154 fstream fincont(fname,ios::in);
155 #if defined(__HP_aCC) || defined(__DECCXX)
156 if(!fincont.rdbuf()->is_open()) return -1;
158 if(!fincont.is_open()) return -1;
162 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
163 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
164 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>fcaq[13][nn]>>
165 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>fcaq[18][nn]>>
166 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>fcaq[23][nn]>>
167 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>fcaq[28][nn]>>
168 fcaq[29][nn]>>fcaq[30][nn]>>fcaq[31][nn]>>fcaq[32][nn]>>fcaq[33][nn])
175 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
176 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
177 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>fcag[13][nn]>>
178 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>fcag[18][nn]>>
179 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>fcag[23][nn]>>
180 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>fcag[28][nn]>>
181 fcag[29][nn]>>fcag[30][nn]>>fcag[31][nn]>>fcag[32][nn]>>fcag[33][nn])
188 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
189 //PH ifstream findisc(fname);
190 fstream findisc(fname,ios::in);
191 #if defined(__HP_aCC) || defined(__DECCXX)
192 if(!findisc.rdbuf()->is_open()) return -1;
194 if(!findisc.is_open()) return -1;
198 while(findisc>>frrr[nn]>>fdaq[nn]) {
203 while(findisc>>frrrg[nn]>>fdag[nn]) {
208 fTablesLoaded = kTRUE;
213 C***************************************************************************
214 C Quenching Weights for Multiple Soft Scattering
219 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
221 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
224 C This package contains quenching weights for gluon radiation in the
225 C multiple soft scattering approximation.
227 C swqmult returns the quenching weight for a quark (ipart=1) or
228 C a gluon (ipart=2) traversing a medium with transport coeficient q and
229 C length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where
230 C wc=0.5*q*L^2 and w is the energy radiated. The output values are
231 C the continuous and discrete (prefactor of the delta function) parts
232 C of the quenching weights.
234 C In order to use this routine, the files cont.all and disc.all need to be
235 C in the working directory.
237 C An initialization of the tables is needed by doing call initmult before
240 C Please, send us any comment:
242 C urs.wiedemann@cern.ch
243 C carlos.salgado@cern.ch
246 C-------------------------------------------------------------------
248 SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete)
250 REAL*8 xx(400), daq(34), caq(34,261), rrr(34)
251 COMMON /dataqua/ xx, daq, caq, rrr
253 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
254 COMMON /dataglu/ xxg, dag, cag, rrrg
256 REAL*8 rrrr,xxxx, continuous, discrete
258 INTEGER nrlow, nrhigh, nxlow, nxhigh
259 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
260 REAL*8 xfraclow, xfrachigh
271 if (rrin.lt.rrr(nr)) then
283 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
284 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
285 if (rrin.gt.10000d0) then
286 rfraclow = dlog(rrhigh/rrin)/dlog(rrhigh/rrlow)
287 rfrachigh = dlog(rrin/rrlow)/dlog(rrhigh/rrlow)
290 if (ipart.eq.1.and.rrin.ge.rrr(1)) then
297 if (ipart.ne.1.and.rrin.ge.rrrg(1)) then
304 if (xxxx.ge.xx(261)) go to 245
306 nxlow = int(xxin/0.01) + 1
308 xfraclow = (xx(nxhigh)-xxin)/0.01
309 xfrachigh = (xxin - xx(nxlow))/0.01
312 clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh)
313 chigh = xfraclow*caq(nrhigh,nxlow)+xfrachigh*caq(nrhigh,nxhigh)
315 clow = xfraclow*cag(nrlow,nxlow)+xfrachigh*cag(nrlow,nxhigh)
316 chigh = xfraclow*cag(nrhigh,nxlow)+xfrachigh*cag(nrhigh,nxhigh)
319 continuous = rfraclow*clow + rfrachigh*chigh
324 discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh)
326 discrete = rfraclow*dag(nrlow) + rfrachigh*dag(nrhigh)
332 REAL*8 xxq(400), daq(34), caq(34,261), rrr(34)
333 COMMON /dataqua/ xxq, daq, caq, rrr
335 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
336 COMMON /dataglu/ xxg, dag, cag, rrrg
338 OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90)
340 read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn),
341 + caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn),
342 + caq(9,nn), caq(10,nn), caq(11,nn), caq(12,nn),
344 + caq(14,nn), caq(15,nn), caq(16,nn), caq(17,nn),
346 + caq(19,nn), caq(20,nn), caq(21,nn), caq(22,nn),
348 + caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn),
350 + caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn),
351 + caq(33,nn), caq(34,nn)
354 read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn),
355 + cag(4,nn), cag(5,nn), cag(6,nn), cag(7,nn), cag(8,nn),
356 + cag(9,nn), cag(10,nn), cag(11,nn), cag(12,nn),
358 + cag(14,nn), cag(15,nn), cag(16,nn), cag(17,nn),
360 + cag(19,nn), cag(20,nn), cag(21,nn), cag(22,nn),
362 + cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn),
364 + cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn),
365 + cag(33,nn), cag(34,nn)
369 OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91)
371 read (21,*) rrr(nn), daq(nn)
374 read (21,*) rrrg(nn), dag(nn)
379 90 PRINT*, 'input - output error'
380 91 PRINT*, 'input - output error #2'
386 =======================================================================
388 Adapted to ROOT macro by A. Dainese - 13/07/2003
389 Ported to class by C. Loizides - 12/02/2004
390 New version for extended R values added - 06/03/2004
393 Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
394 Double_t &continuous,Double_t &discrete) const
396 // Calculate Multiple Scattering approx.
397 // weights for given parton type,
398 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
404 //read-in data before first call
406 Error("CalcMult","Tables are not loaded.");
410 Error("CalcMult","Tables are not loaded for Multiple Scattering.");
414 Double_t rrin = rrrr;
415 Double_t xxin = xxxx;
417 if(xxin>fxx[260]) return -1;
418 Int_t nxlow = (Int_t)(xxin/0.01) + 1;
419 Int_t nxhigh = nxlow + 1;
420 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.01;
421 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01;
424 if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD
425 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
427 Int_t nrlow=0,nrhigh=0;
428 Double_t rrhigh=0,rrlow=0;
429 for(Int_t nr=1; nr<=34; nr++) {
430 if(rrin<frrr[nr-1]) {
433 rrhigh = frrr[nr-1-1];
443 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
444 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
447 rfraclow = TMath::Log2(rrhigh/rrin)/TMath::Log2(rrhigh/rrlow);
448 rfrachigh = TMath::Log2(rrin/rrlow)/TMath::Log2(rrhigh/rrlow);
450 if((ipart==1) && (rrin>=frrr[0]))
457 if((ipart==2) && (rrin>=frrrg[0]))
465 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
467 Double_t clow=0,chigh=0;
469 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
470 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
472 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
473 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
476 continuous = rfraclow*clow + rfrachigh*chigh;
477 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
478 //rfraclow,clow,rfrachigh,chigh,continuous);
481 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
483 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
489 Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *discall)
491 // read in tables for Single Hard Approx.
492 // path to continuum and to discrete part
494 fTablesLoaded = kFALSE;
498 sprintf(fname,"%s",gSystem->ExpandPathName(contall));
499 //PH ifstream fincont(fname);
500 fstream fincont(fname,ios::in);
501 #if defined(__HP_aCC) || defined(__DECCXX)
502 if(!fincont.rdbuf()->is_open()) return -1;
504 if(!fincont.is_open()) return -1;
508 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
509 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
510 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>
512 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>
514 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>
516 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>
525 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
526 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
527 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>
529 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>
531 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>
533 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>
541 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
542 //PH ifstream findisc(fname);
543 fstream findisc(fname,ios::in);
544 #if defined(__HP_aCC) || defined(__DECCXX)
545 if(!findisc.rdbuf()->is_open()) return -1;
547 if(!findisc.is_open()) return -1;
551 while(findisc>>frrr[nn]>>fdaq[nn]) {
556 while(findisc>>frrrg[nn]>>fdag[nn]) {
562 fTablesLoaded = kTRUE;
567 C***************************************************************************
568 C Quenching Weights for Single Hard Scattering
573 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
575 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
578 C This package contains quenching weights for gluon radiation in the
579 C single hard scattering approximation.
581 C swqlin returns the quenching weight for a quark (ipart=1) or
582 C a gluon (ipart=2) traversing a medium with Debye screening mass mu and
583 C length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where
584 C wc=0.5*mu^2*L and w is the energy radiated. The output values are
585 C the continuous and discrete (prefactor of the delta function) parts
586 C of the quenching weights.
588 C In order to use this routine, the files contlin.all and disclin.all
589 C need to be in the working directory.
591 C An initialization of the tables is needed by doing call initlin before
594 C Please, send us any comment:
596 C urs.wiedemann@cern.ch
597 C carlos.salgado@cern.ch
600 C-------------------------------------------------------------------
603 SUBROUTINE swqlin(ipart,rrrr,xxxx,continuous,discrete)
605 REAL*8 xx(400), dalq(30), calq(30,261), rrr(30)
606 COMMON /datalinqua/ xx, dalq, calq, rrr
608 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
609 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
611 REAL*8 rrrr,xxxx, continuous, discrete
613 INTEGER nrlow, nrhigh, nxlow, nxhigh
614 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
615 REAL*8 xfraclow, xfrachigh
621 nxlow = int(xxin/0.038) + 1
623 xfraclow = (xx(nxhigh)-xxin)/0.038
624 xfrachigh = (xxin - xx(nxlow))/0.038
627 if (rrin.lt.rrr(nr)) then
639 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
640 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
643 clow = xfraclow*calq(nrlow,nxlow)+xfrachigh*calq(nrlow,nxhigh)
644 chigh = xfraclow*calq(nrhigh,nxlow)+xfrachigh*calq(nrhigh,nxhigh)
646 clow = xfraclow*calg(nrlow,nxlow)+xfrachigh*calg(nrlow,nxhigh)
647 chigh = xfraclow*calg(nrhigh,nxlow)+xfrachigh*calg(nrhigh,nxhigh)
650 continuous = rfraclow*clow + rfrachigh*chigh
653 discrete = rfraclow*dalq(nrlow) + rfrachigh*dalq(nrhigh)
655 discrete = rfraclow*dalg(nrlow) + rfrachigh*dalg(nrhigh)
661 REAL*8 xxlq(400), dalq(30), calq(30,261), rrr(30)
662 COMMON /datalinqua/ xxlq, dalq, calq, rrr
664 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
665 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
667 OPEN(UNIT=20,FILE='contlin.all',STATUS='OLD',ERR=90)
669 read (20,*) xxlq(nn), calq(1,nn), calq(2,nn), calq(3,nn),
670 + calq(4,nn), calq(5,nn), calq(6,nn), calq(7,nn), calq(8,nn),
671 + calq(9,nn), calq(10,nn), calq(11,nn), calq(12,nn),
673 + calq(14,nn), calq(15,nn), calq(16,nn), calq(17,nn),
675 + calq(19,nn), calq(20,nn), calq(21,nn), calq(22,nn),
677 + calq(24,nn), calq(25,nn), calq(26,nn), calq(27,nn),
679 + calq(29,nn), calq(30,nn)
682 read (20,*) xxlg(nn), calg(1,nn), calg(2,nn), calg(3,nn),
683 + calg(4,nn), calg(5,nn), calg(6,nn), calg(7,nn), calg(8,nn),
684 + calg(9,nn), calg(10,nn), calg(11,nn), calg(12,nn),
686 + calg(14,nn), calg(15,nn), calg(16,nn), calg(17,nn),
688 + calg(19,nn), calg(20,nn), calg(21,nn), calg(22,nn),
690 + calg(24,nn), calg(25,nn), calg(26,nn), calg(27,nn),
692 + calg(29,nn), calg(30,nn)
696 OPEN(UNIT=21,FILE='disclin.all',STATUS='OLD',ERR=91)
698 read (21,*) rrr(nn), dalq(nn)
701 read (21,*) rrrlg(nn), dalg(nn)
706 90 PRINT*, 'input - output error'
707 91 PRINT*, 'input - output error #2'
712 =======================================================================
714 Ported to class by C. Loizides - 17/02/2004
718 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx,
719 Double_t &continuous,Double_t &discrete) const
721 // calculate Single Hard approx.
722 // weights for given parton type,
723 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
725 // read-in data before first call
727 Error("CalcSingleHard","Tables are not loaded.");
731 Error("CalcSingleHard","Tables are not loaded for Single Hard Scattering.");
735 Double_t rrin = rrrr;
736 Double_t xxin = xxxx;
738 Int_t nxlow = (Int_t)(xxin/0.038) + 1;
739 Int_t nxhigh = nxlow + 1;
740 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.038;
741 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.038;
744 if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD
745 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
747 Int_t nrlow=0,nrhigh=0;
748 Double_t rrhigh=0,rrlow=0;
749 for(Int_t nr=1; nr<=30; nr++) {
750 if(rrin<frrr[nr-1]) {
753 rrhigh = frrr[nr-1-1];
763 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
764 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
766 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
768 Double_t clow=0,chigh=0;
770 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
771 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
773 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
774 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
777 continuous = rfraclow*clow + rfrachigh*chigh;
778 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
779 // rfraclow,clow,rfrachigh,chigh,continuous);
782 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
784 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
790 Int_t AliQuenchingWeights::CalcMult(Int_t ipart,
791 Double_t w,Double_t qtransp,Double_t length,
792 Double_t &continuous,Double_t &discrete) const
794 // Calculate Multiple Scattering approx.
795 // weights for given parton type,
796 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
798 Double_t wc=CalcWC(qtransp,length);
799 Double_t rrrr=CalcR(wc,length);
801 return CalcMult(ipart,rrrr,xxxx,continuous,discrete);
804 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
805 Double_t w,Double_t mu,Double_t length,
806 Double_t &continuous,Double_t &discrete) const
808 // calculate Single Hard approx.
809 // weights for given parton type,
810 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
812 Double_t wcbar=CalcWCbar(mu,length);
813 Double_t rrrr=CalcR(wcbar,length);
814 Double_t xxxx=w/wcbar;
815 return CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
818 Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const
820 //calculate R value and
821 //check if it is less then maximum
823 Double_t R = wc*l*fgkConvFmToInvGeV;
825 Warning("CalcR","Value of R = %.2f; should be less than %.2f",R,fgkRMax);
831 Double_t AliQuenchingWeights::CalcRk(Double_t k, Double_t I0, Double_t I1) const
833 //calculate R value and
834 //check if it is less then maximum
836 Double_t R = fgkRMax;
840 Warning("CalcRk","Value of R = %.2f; should be less than %.2f",R,fgkRMax);
846 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const
848 // return DeltaE for MS or SH scattering
849 // for given parton type, length and energy
850 // Dependant on ECM (energy constraint method)
851 // e is used to determine where to set bins to zero.
854 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
857 if((ipart<1) || (ipart>2)) {
858 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
862 Int_t l=GetIndex(length);
865 if(fECMethod==kReweight){
869 ret=fHistos[ipart-1][l-1]->GetRandom();
871 Warning("GetELossRandomK",
872 "Aborted reweighting; maximum loss assigned after 1e5 trials.");
879 Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
884 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const
886 //return quenched parton energy
887 //for given parton type, length and energy
889 Double_t loss=GetELossRandom(ipart,length,e);
893 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const
895 // return DeltaE for MS or SH scattering
896 // for given parton type, length distribution and energy
897 // Dependant on ECM (energy constraint method)
898 // e is used to determine where to set bins to zero.
901 Warning("GetELossRandom","Pointer to length distribution is NULL.");
904 Double_t ell=hell->GetRandom();
905 return GetELossRandom(ipart,ell,e);
908 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e) const
910 //return quenched parton energy
911 //for given parton type, length distribution and energy
913 Double_t loss=GetELossRandom(ipart,hell,e);
917 Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
919 // return DeltaE for new dynamic version
920 // for given parton type, I0 and I1 value and energy
921 // Dependant on ECM (energy constraint method)
922 // e is used to determine where to set bins to zero.
924 // read-in data before first call
926 Fatal("GetELossRandomK","Tables are not loaded.");
929 if((ipart<1) || (ipart>2)) {
930 Fatal("GetELossRandomK","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
934 Double_t R=CalcRk(I0,I1);
936 Fatal("GetELossRandomK","R should not be negative");
939 Double_t wc=CalcWCk(I1);
941 Fatal("GetELossRandomK","wc should be greater than zero");
944 if(SampleEnergyLoss(ipart,R)!=0){
945 cout << ipart << " " << R << endl;
946 Fatal("GetELossRandomK","Could not sample energy loss");
950 if(fECMethod==kReweight){
954 ret=fHisto->GetRandom();
956 Warning("GetELossRandomK",
957 "Aborted reweighting; maximum loss assigned after 1e5 trials.");
965 Double_t ret=fHisto->GetRandom()*wc;
970 Double_t AliQuenchingWeights::CalcQuenchedEnergyK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
972 //return quenched parton energy
973 //for given parton type, I0 and I1 value and energy
975 Double_t loss=GetELossRandomK(ipart,I0,I1,e);
979 Int_t AliQuenchingWeights::SampleEnergyLoss()
981 // Has to be called to fill the histograms
983 // For stored values fQTransport loop over
984 // particle type and length = 1 to fMaxLength (fm)
985 // to fill energy loss histos
987 // Take histogram of continuous weights
988 // Take discrete_weight
989 // If discrete_weight > 1, put all channels to 0, except channel 1
990 // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral
992 // read-in data before first call
994 Error("SampleEnergyLoss","Tables are not loaded.");
999 Int_t lmax=CalcLengthMax(fQTransport);
1000 if(fLengthMax>lmax){
1001 Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
1005 Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested.");
1009 fHistos=new TH1F**[2];
1010 fHistos[0]=new TH1F*[4*fLengthMax];
1011 fHistos[1]=new TH1F*[4*fLengthMax];
1012 fLengthMaxOld=fLengthMax; //remember old value in case
1013 //user wants to reset
1016 Char_t meddesc[100];
1018 medvalue=(Int_t)(fQTransport*1000.);
1019 sprintf(meddesc,"MS");
1021 medvalue=(Int_t)(fMu*1000.);
1022 sprintf(meddesc,"SH");
1025 for(Int_t ipart=1;ipart<=2;ipart++){
1026 for(Int_t l=1;l<=4*fLengthMax;l++){
1028 sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
1030 Double_t wc = CalcWC(len);
1031 fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgBins,0.,fgMaxBin*wc);
1032 fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
1033 fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
1034 fHistos[ipart-1][l-1]->SetLineColor(4);
1036 Double_t rrrr = CalcR(wc,len);
1037 Double_t discrete=0.;
1038 // loop on histogram channels
1039 for(Int_t bin=1; bin<=fgBins; bin++) {
1040 Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
1041 Double_t continuous;
1043 CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1045 CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1046 fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
1048 // add discrete part to distribution
1050 for(Int_t bin=2;bin<=fgBins;bin++)
1051 fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
1053 Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgBins);
1054 fHistos[ipart-1][l-1]->Fill(0.,val);
1056 Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgBins);
1057 fHistos[ipart-1][l-1]->Scale(1./hint);
1063 Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t R)
1065 // Sample energy loss directly for one particle type
1066 // choses R (safe it and keep it until next call of function)
1068 // read-in data before first call
1070 Error("SampleEnergyLoss","Tables are not loaded.");
1074 Double_t discrete=0.;
1075 Double_t continuous=0;;
1077 Double_t xxxx = fHisto->GetBinCenter(bin);
1079 CalcMult(ipart,R,xxxx,continuous,discrete);
1081 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1084 fHisto->SetBinContent(1,1.);
1085 for(Int_t bin=2;bin<=fgBins;bin++)
1086 fHisto->SetBinContent(bin,0.);
1090 fHisto->SetBinContent(bin,continuous);
1091 for(Int_t bin=2; bin<=fgBins; bin++) {
1092 xxxx = fHisto->GetBinCenter(bin);
1094 CalcMult(ipart,R,xxxx,continuous,discrete);
1096 CalcSingleHard(ipart,R,xxxx,continuous,discrete);
1097 fHisto->SetBinContent(bin,continuous);
1099 // add discrete part to distribution
1100 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgBins);
1101 fHisto->Fill(0.,val);
1103 Double_t hint=fHisto->Integral(1,fgBins);
1105 fHisto->Scale(1./hint);
1107 cout << discrete << " " << hint << " " << continuous << endl;
1113 const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) const
1115 //return quenching histograms
1116 //for ipart and length
1119 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
1122 if((ipart<1) || (ipart>2)) {
1123 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
1127 Int_t l=GetIndex(length);
1129 return fHistos[ipart-1][l-1];
1132 TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const
1134 // ipart = 1 for quark, 2 for gluon
1135 // medval a) qtransp = transport coefficient (GeV^2/fm)
1136 // b) mu = Debye mass (GeV)
1137 // length = path length in medium (fm)
1138 // Get from SW tables:
1139 // - continuous weight, as a function of dE/wc
1142 Char_t meddesc[100];
1144 wc=CalcWC(medval,length);
1145 sprintf(meddesc,"MS");
1147 wc=CalcWCbar(medval,length);
1148 sprintf(meddesc,"SH");
1152 sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart,
1153 (Int_t)(medval*1000.),(Int_t)length);
1155 TH1F *hist = new TH1F("hist",hname,fgBins,0.,fgMaxBin*wc);
1156 hist->SetXTitle("#Delta E [GeV]");
1157 hist->SetYTitle("p(#Delta E)");
1158 hist->SetLineColor(4);
1160 Double_t rrrr = CalcR(wc,length);
1161 //loop on histogram channels
1162 for(Int_t bin=1; bin<=fgBins; bin++) {
1163 Double_t xxxx = hist->GetBinCenter(bin)/wc;
1164 Double_t continuous,discrete;
1166 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1167 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1172 hist->SetBinContent(bin,continuous);
1177 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const
1179 // ipart = 1 for quark, 2 for gluon
1180 // medval a) qtransp = transport coefficient (GeV^2/fm)
1181 // b) mu = Debye mass (GeV)
1182 // length = path length in medium (fm)
1183 // Get from SW tables:
1184 // - continuous weight, as a function of dE/wc
1187 Char_t meddesc[100];
1189 wc=CalcWC(medval,length);
1190 sprintf(meddesc,"MS");
1192 wc=CalcWCbar(medval,length);
1193 sprintf(meddesc,"SH");
1197 sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart,
1198 (Int_t)(medval*1000.),(Int_t)length);
1200 TH1F *histx = new TH1F("histx",hname,fgBins,0.,fgMaxBin);
1201 histx->SetXTitle("x = #Delta E/#omega_{c}");
1203 histx->SetYTitle("p(#Delta E/#omega_{c})");
1205 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1206 histx->SetLineColor(4);
1208 Double_t rrrr = CalcR(wc,length);
1209 //loop on histogram channels
1210 for(Int_t bin=1; bin<=fgBins; bin++) {
1211 Double_t xxxx = histx->GetBinCenter(bin);
1212 Double_t continuous,discrete;
1214 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1215 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1220 histx->SetBinContent(bin,continuous);
1225 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t R) const
1227 // compute P(E) distribution for
1228 // given ipart = 1 for quark, 2 for gluon
1231 Char_t meddesc[100];
1233 sprintf(meddesc,"MS");
1235 sprintf(meddesc,"SH");
1239 sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,R);
1240 TH1F *histx = new TH1F("histx",hname,fgBins,0.,fgMaxBin);
1241 histx->SetXTitle("x = #Delta E/#omega_{c}");
1243 histx->SetYTitle("p(#Delta E/#omega_{c})");
1245 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1246 histx->SetLineColor(4);
1249 Double_t continuous=0.,discrete=0.;
1250 //loop on histogram channels
1251 for(Int_t bin=1; bin<=fgBins; bin++) {
1252 Double_t xxxx = histx->GetBinCenter(bin);
1254 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1255 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1260 histx->SetBinContent(bin,continuous);
1263 //add discrete part to distribution
1265 for(Int_t bin=2;bin<=fgBins;bin++)
1266 histx->SetBinContent(bin,0.);
1268 Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgBins);
1269 histx->Fill(0.,val);
1271 Double_t hint=histx->Integral(1,fgBins);
1272 if(hint!=0) histx->Scale(1./hint);
1277 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
1279 // compute energy loss histogram for
1280 // parton type, medium value, length and energy
1282 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1284 dummy->SetQTransport(medval);
1287 dummy->SetMu(medval);
1288 dummy->InitSingleHard();
1290 dummy->SampleEnergyLoss();
1295 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1296 sprintf(hname,"hLossQuarks");
1298 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1299 sprintf(hname,"hLossGluons");
1302 TH1F *h = new TH1F(hname,name,250,0,250);
1303 for(Int_t i=0;i<100000;i++){
1304 //if(i % 1000 == 0) cout << "." << flush;
1305 Double_t loss=dummy->GetELossRandom(ipart,l,e);
1313 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
1315 // compute energy loss histogram for
1316 // parton type, medium value,
1317 // length distribution and energy
1319 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1321 dummy->SetQTransport(medval);
1324 dummy->SetMu(medval);
1325 dummy->InitSingleHard();
1327 dummy->SampleEnergyLoss();
1332 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1333 sprintf(hname,"hLossQuarks");
1335 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1336 sprintf(hname,"hLossGluons");
1339 TH1F *h = new TH1F(hname,name,250,0,250);
1340 for(Int_t i=0;i<100000;i++){
1341 //if(i % 1000 == 0) cout << "." << flush;
1342 Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
1350 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t R) const
1352 // compute energy loss histogram for
1353 // parton type and given R
1355 TH1F *dummy = ComputeQWHistoX(ipart,R);
1356 if(!dummy) return 0;
1359 sprintf(hname,"hELossHistox_%d_%.2f",ipart,R);
1360 TH1F *histx = new TH1F("histxr",hname,fgBins,0.,fgMaxBin);
1361 for(Int_t i=0;i<100000;i++){
1362 //if(i % 1000 == 0) cout << "." << flush;
1363 Double_t loss=dummy->GetRandom();
1370 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,Double_t l) const
1372 // compute average energy loss for
1373 // parton type, medium value, length and energy
1375 TH1F *dummy = ComputeELossHisto(ipart,medval,l);
1376 if(!dummy) return 0;
1377 Double_t ret=dummy->GetMean();
1382 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEll) const
1384 // compute average energy loss for
1385 // parton type, medium value,
1386 // length distribution and energy
1388 TH1F *dummy = ComputeELossHisto(ipart,medval,hEll);
1389 if(!dummy) return 0;
1390 Double_t ret=dummy->GetMean();
1395 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t R) const
1397 // compute average energy loss over wc
1398 // for parton type and given R
1400 TH1F *dummy = ComputeELossHisto(ipart,R);
1401 if(!dummy) return 0;
1402 Double_t ret=dummy->GetMean();
1407 void AliQuenchingWeights::PlotDiscreteWeights(Double_t len) const
1409 // plot discrete weights for given length
1413 c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
1415 c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
1418 TH2F *hframe = new TH2F("hdisc","",2,0,5.1,2,0,1);
1419 hframe->SetStats(0);
1421 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1423 hframe->SetXTitle("#mu [GeV]");
1424 hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
1427 TGraph *gq=new TGraph(20);
1430 for(Double_t q=0.05;q<=5.05;q+=0.25){
1432 CalcMult(1,1.0,q,len,cont,disc);
1433 gq->SetPoint(i,q,disc);i++;
1436 for(Double_t m=0.05;m<=5.05;m+=0.25){
1438 CalcSingleHard(1,1.0,m,len,cont, disc);
1439 gq->SetPoint(i,m,disc);i++;
1442 gq->SetMarkerStyle(20);
1445 TGraph *gg=new TGraph(20);
1448 for(Double_t q=0.05;q<=5.05;q+=0.25){
1450 CalcMult(2,1.0,q,len,cont,disc);
1451 gg->SetPoint(i,q,disc);i++;
1454 for(Double_t m=0.05;m<=5.05;m+=0.25){
1456 CalcSingleHard(2,1.0,m,len,cont,disc);
1457 gg->SetPoint(i,m,disc);i++;
1460 gg->SetMarkerStyle(24);
1463 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1464 l1a->SetFillStyle(0);
1465 l1a->SetBorderSize(0);
1467 sprintf(label,"L = %.1f fm",len);
1468 l1a->AddEntry(gq,label,"");
1469 l1a->AddEntry(gq,"quark","pl");
1470 l1a->AddEntry(gg,"gluon","pl");
1476 void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
1478 // plot continous weights for
1479 // given parton type and length
1486 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1488 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1490 medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1491 sprintf(name,"ccont-ms-%d",itype);
1494 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1496 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1498 medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1499 sprintf(name,"ccont-ms-%d",itype);
1502 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1504 TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1506 h1->SetTitle(title);
1508 h1->SetLineColor(1);
1510 TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1512 h2->SetLineColor(2);
1513 h2->DrawCopy("SAME");
1514 TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1516 h3->SetLineColor(3);
1517 h3->DrawCopy("SAME");
1519 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1520 l1a->SetFillStyle(0);
1521 l1a->SetBorderSize(0);
1523 sprintf(label,"L = %.1f fm",ell);
1524 l1a->AddEntry(h1,label,"");
1526 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1527 l1a->AddEntry(h1,label,"pl");
1528 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1529 l1a->AddEntry(h2,label,"pl");
1530 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1531 l1a->AddEntry(h3,label,"pl");
1533 sprintf(label,"#mu = %.1f GeV",medvals[0]);
1534 l1a->AddEntry(h1,label,"pl");
1535 sprintf(label,"#mu = %.1f GeV",medvals[1]);
1536 l1a->AddEntry(h2,label,"pl");
1537 sprintf(label,"#mu = %.1f GeV",medvals[2]);
1538 l1a->AddEntry(h3,label,"pl");
1545 void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
1547 // plot continous weights for
1548 // given parton type and medium value
1554 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1556 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1558 sprintf(name,"ccont2-ms-%d",itype);
1561 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1563 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1565 sprintf(name,"ccont2-sh-%d",itype);
1567 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1569 TH1F *h1=ComputeQWHisto(itype,medval,8);
1571 h1->SetTitle(title);
1573 h1->SetLineColor(1);
1575 TH1F *h2=ComputeQWHisto(itype,medval,5);
1577 h2->SetLineColor(2);
1578 h2->DrawCopy("SAME");
1579 TH1F *h3=ComputeQWHisto(itype,medval,2);
1581 h3->SetLineColor(3);
1582 h3->DrawCopy("SAME");
1584 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1585 l1a->SetFillStyle(0);
1586 l1a->SetBorderSize(0);
1589 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval);
1591 sprintf(label,"#mu = %.1f GeV",medval);
1593 l1a->AddEntry(h1,label,"");
1594 l1a->AddEntry(h1,"L = 8 fm","pl");
1595 l1a->AddEntry(h2,"L = 5 fm","pl");
1596 l1a->AddEntry(h3,"L = 2 fm","pl");
1602 void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t e) const
1604 // plot average energy loss for given length
1605 // and parton energy
1608 Error("PlotAvgELoss","Tables are not loaded.");
1615 sprintf(title,"Average Energy Loss for Multiple Scattering");
1616 sprintf(name,"cavgelossms");
1618 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1619 sprintf(name,"cavgelosssh");
1622 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1624 TH2F *hframe = new TH2F("avgloss",title,2,0,5.1,2,0,100);
1625 hframe->SetStats(0);
1627 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1629 hframe->SetXTitle("#mu [GeV]");
1630 hframe->SetYTitle("<E_{loss}> [GeV]");
1633 TGraph *gq=new TGraph(20);
1635 for(Double_t v=0.05;v<=5.05;v+=0.25){
1636 TH1F *dummy=ComputeELossHisto(1,v,len,e);
1637 Double_t avgloss=dummy->GetMean();
1638 gq->SetPoint(i,v,avgloss);i++;
1641 gq->SetMarkerStyle(20);
1644 TGraph *gg=new TGraph(20);
1646 for(Double_t v=0.05;v<=5.05;v+=0.25){
1647 TH1F *dummy=ComputeELossHisto(2,v,len,e);
1648 Double_t avgloss=dummy->GetMean();
1649 gg->SetPoint(i,v,avgloss);i++;
1652 gg->SetMarkerStyle(24);
1655 TGraph *gratio=new TGraph(20);
1656 for(Int_t i=0;i<20;i++){
1658 gg->GetPoint(i,x,y);
1659 gq->GetPoint(i,x2,y2);
1661 gratio->SetPoint(i,x,y/y2*10/2.25);
1662 else gratio->SetPoint(i,x,0);
1664 gratio->SetLineStyle(4);
1666 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1667 l1a->SetFillStyle(0);
1668 l1a->SetBorderSize(0);
1670 sprintf(label,"L = %.1f fm",len);
1671 l1a->AddEntry(gq,label,"");
1672 l1a->AddEntry(gq,"quark","pl");
1673 l1a->AddEntry(gg,"gluon","pl");
1674 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1680 void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
1682 // plot average energy loss for given
1683 // length distribution and parton energy
1686 Error("PlotAvgELossVs","Tables are not loaded.");
1693 sprintf(title,"Average Energy Loss for Multiple Scattering");
1694 sprintf(name,"cavgelossms2");
1696 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1697 sprintf(name,"cavgelosssh2");
1700 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1702 TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1703 hframe->SetStats(0);
1705 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1707 hframe->SetXTitle("#mu [GeV]");
1708 hframe->SetYTitle("<E_{loss}> [GeV]");
1711 TGraph *gq=new TGraph(20);
1713 for(Double_t v=0.05;v<=5.05;v+=0.25){
1714 TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1715 Double_t avgloss=dummy->GetMean();
1716 gq->SetPoint(i,v,avgloss);i++;
1719 gq->SetMarkerStyle(20);
1722 TGraph *gg=new TGraph(20);
1724 for(Double_t v=0.05;v<=5.05;v+=0.25){
1725 TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1726 Double_t avgloss=dummy->GetMean();
1727 gg->SetPoint(i,v,avgloss);i++;
1730 gg->SetMarkerStyle(24);
1733 TGraph *gratio=new TGraph(20);
1734 for(Int_t i=0;i<20;i++){
1736 gg->GetPoint(i,x,y);
1737 gq->GetPoint(i,x2,y2);
1739 gratio->SetPoint(i,x,y/y2*10/2.25);
1740 else gratio->SetPoint(i,x,0);
1742 gratio->SetLineStyle(4);
1745 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1746 l1a->SetFillStyle(0);
1747 l1a->SetBorderSize(0);
1749 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1750 l1a->AddEntry(gq,label,"");
1751 l1a->AddEntry(gq,"quark","pl");
1752 l1a->AddEntry(gg,"gluon","pl");
1753 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1759 void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const
1761 // plot average energy loss versus ell
1764 Error("PlotAvgELossVsEll","Tables are not loaded.");
1772 sprintf(title,"Average Energy Loss for Multiple Scattering");
1773 sprintf(name,"cavgelosslms");
1776 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1777 sprintf(name,"cavgelosslsh");
1781 TCanvas *c = new TCanvas(name,title,0,0,600,400);
1783 TH2F *hframe = new TH2F("avglossell",title,2,0,fLengthMax,2,0,250);
1784 hframe->SetStats(0);
1785 hframe->SetXTitle("length [fm]");
1786 hframe->SetYTitle("<E_{loss}> [GeV]");
1789 TGraph *gq=new TGraph((Int_t)fLengthMax*4);
1791 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
1792 TH1F *dummy=ComputeELossHisto(1,medval,v,e);
1793 Double_t avgloss=dummy->GetMean();
1794 gq->SetPoint(i,v,avgloss);i++;
1797 gq->SetMarkerStyle(20);
1800 TGraph *gg=new TGraph((Int_t)fLengthMax*4);
1802 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
1803 TH1F *dummy=ComputeELossHisto(2,medval,v,e);
1804 Double_t avgloss=dummy->GetMean();
1805 gg->SetPoint(i,v,avgloss);i++;
1808 gg->SetMarkerStyle(24);
1811 TGraph *gratio=new TGraph((Int_t)fLengthMax*4);
1812 for(Int_t i=0;i<=(Int_t)fLengthMax*4;i++){
1814 gg->GetPoint(i,x,y);
1815 gq->GetPoint(i,x2,y2);
1817 gratio->SetPoint(i,x,y/y2*100/2.25);
1818 else gratio->SetPoint(i,x,0);
1820 gratio->SetLineStyle(1);
1821 gratio->SetLineWidth(2);
1823 TLegend *l1a = new TLegend(0.15,0.65,.60,0.85);
1824 l1a->SetFillStyle(0);
1825 l1a->SetBorderSize(0);
1828 sprintf(label,"#hat{q} = %.2f GeV^{2}/fm",medval);
1830 sprintf(label,"#mu = %.2f GeV",medval);
1831 l1a->AddEntry(gq,label,"");
1832 l1a->AddEntry(gq,"quark","pl");
1833 l1a->AddEntry(gg,"gluon","pl");
1834 l1a->AddEntry(gratio,"gluon/quark/2.25*100","pl");
1837 TF1 *f=new TF1("f","100",0,fLengthMax);
1844 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
1846 // plot relative energy loss for given
1847 // length and parton energy versus pt
1850 Error("PlotAvgELossVsPt","Tables are not loaded.");
1857 sprintf(title,"Relative Energy Loss for Multiple Scattering");
1858 sprintf(name,"cavgelossvsptms");
1860 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
1861 sprintf(name,"cavgelossvsptsh");
1864 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1866 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
1867 hframe->SetStats(0);
1868 hframe->SetXTitle("p_{T} [GeV]");
1869 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
1872 TGraph *gq=new TGraph(40);
1874 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1875 TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
1876 Double_t avgloss=dummy->GetMean();
1877 gq->SetPoint(i,pt,avgloss/pt);i++;
1880 gq->SetMarkerStyle(20);
1883 TGraph *gg=new TGraph(40);
1885 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1886 TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
1887 Double_t avgloss=dummy->GetMean();
1888 gg->SetPoint(i,pt,avgloss/pt);i++;
1891 gg->SetMarkerStyle(24);
1894 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1895 l1a->SetFillStyle(0);
1896 l1a->SetBorderSize(0);
1898 sprintf(label,"L = %.1f fm",len);
1899 l1a->AddEntry(gq,label,"");
1900 l1a->AddEntry(gq,"quark","pl");
1901 l1a->AddEntry(gg,"gluon","pl");
1907 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
1909 // plot relative energy loss for given
1910 // length distribution and parton energy versus pt
1913 Error("PlotAvgELossVsPt","Tables are not loaded.");
1920 sprintf(title,"Relative Energy Loss for Multiple Scattering");
1921 sprintf(name,"cavgelossvsptms2");
1923 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
1924 sprintf(name,"cavgelossvsptsh2");
1926 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1928 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
1929 hframe->SetStats(0);
1930 hframe->SetXTitle("p_{T} [GeV]");
1931 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
1934 TGraph *gq=new TGraph(40);
1936 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1937 TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
1938 Double_t avgloss=dummy->GetMean();
1939 gq->SetPoint(i,pt,avgloss/pt);i++;
1942 gq->SetMarkerStyle(20);
1945 TGraph *gg=new TGraph(40);
1947 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
1948 TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
1949 Double_t avgloss=dummy->GetMean();
1950 gg->SetPoint(i,pt,avgloss/pt);i++;
1953 gg->SetMarkerStyle(24);
1956 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1957 l1a->SetFillStyle(0);
1958 l1a->SetBorderSize(0);
1960 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1961 l1a->AddEntry(gq,label,"");
1962 l1a->AddEntry(gq,"quark","pl");
1963 l1a->AddEntry(gg,"gluon","pl");
1969 Int_t AliQuenchingWeights::GetIndex(Double_t len) const
1971 if(len>fLengthMax) len=fLengthMax;
1973 Int_t l=Int_t(len/0.25);
1974 if((len-l*0.5)>0.125) l++;