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()
64 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)
86 fInstanceNumber(fgCounter++),
100 fTablesLoaded=kFALSE;
103 fMultSoft=a.GetMultSoft();;
106 fQTransport=a.GetQTransport();
107 fECMethod=(kECMethod)a.GetECMethod();
108 fLengthMax=a.GetLengthMax();
109 fInstanceNumber=fgCounter++;
111 sprintf(name,"hhistoqw_%d",fInstanceNumber);
112 fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin);
113 for(Int_t bin=1;bin<=fgkBins;bin++)
114 fHisto->SetBinContent(bin,0.);
116 //Missing in the class is the pathname
117 //to the tables, can be added if needed
120 AliQuenchingWeights::~AliQuenchingWeights()
126 void AliQuenchingWeights::Reset()
128 //reset tables if there were used
131 for(Int_t l=0;l<4*fLengthMaxOld;l++){
132 delete fHistos[0][l];
133 delete fHistos[1][l];
140 void AliQuenchingWeights::SetECMethod(kECMethod type)
142 //set energy constraint method
145 if(fECMethod==kDefault)
146 Info("SetECMethod","Energy Constraint Method set to DEFAULT:\nIf (sampled energy loss > parton energy) then sampled energy loss = parton energy.");
147 else if(fECMethod==kReweight)
148 Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy.");
149 else Info("SetECMethod","Energy Constraint Method set to REWEIGHTCONT:\nRequire sampled energy loss <= parton energy (only implemented for FAST method.");
152 Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
154 // read in tables for multiple scattering approximation
155 // path to continuum and to discrete part
157 fTablesLoaded = kFALSE;
161 sprintf(fname,"%s",gSystem->ExpandPathName(contall));
162 //PH ifstream fincont(fname);
163 fstream fincont(fname,ios::in);
164 #if defined(__HP_aCC) || defined(__DECCXX)
165 if(!fincont.rdbuf()->is_open()) return -1;
167 if(!fincont.is_open()) return -1;
171 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
172 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
173 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>fcaq[13][nn]>>
174 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>fcaq[18][nn]>>
175 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>fcaq[23][nn]>>
176 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>fcaq[28][nn]>>
177 fcaq[29][nn]>>fcaq[30][nn]>>fcaq[31][nn]>>fcaq[32][nn]>>fcaq[33][nn])
184 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
185 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
186 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>fcag[13][nn]>>
187 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>fcag[18][nn]>>
188 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>fcag[23][nn]>>
189 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>fcag[28][nn]>>
190 fcag[29][nn]>>fcag[30][nn]>>fcag[31][nn]>>fcag[32][nn]>>fcag[33][nn])
197 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
198 //PH ifstream findisc(fname);
199 fstream findisc(fname,ios::in);
200 #if defined(__HP_aCC) || defined(__DECCXX)
201 if(!findisc.rdbuf()->is_open()) return -1;
203 if(!findisc.is_open()) return -1;
207 while(findisc>>frrr[nn]>>fdaq[nn]) {
212 while(findisc>>frrrg[nn]>>fdag[nn]) {
217 fTablesLoaded = kTRUE;
222 C***************************************************************************
223 C Quenching Weights for Multiple Soft Scattering
228 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
230 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
233 C This package contains quenching weights for gluon radiation in the
234 C multiple soft scattering approximation.
236 C swqmult returns the quenching weight for a quark (ipart=1) or
237 C a gluon (ipart=2) traversing a medium with transport coeficient q and
238 C length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where
239 C wc=0.5*q*L^2 and w is the energy radiated. The output values are
240 C the continuous and discrete (prefactor of the delta function) parts
241 C of the quenching weights.
243 C In order to use this routine, the files cont.all and disc.all need to be
244 C in the working directory.
246 C An initialization of the tables is needed by doing call initmult before
249 C Please, send us any comment:
251 C urs.wiedemann@cern.ch
252 C carlos.salgado@cern.ch
255 C-------------------------------------------------------------------
257 SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete)
259 REAL*8 xx(400), daq(34), caq(34,261), rrr(34)
260 COMMON /dataqua/ xx, daq, caq, rrr
262 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
263 COMMON /dataglu/ xxg, dag, cag, rrrg
265 REAL*8 rrrr,xxxx, continuous, discrete
267 INTEGER nrlow, nrhigh, nxlow, nxhigh
268 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
269 REAL*8 xfraclow, xfrachigh
280 if (rrin.lt.rrr(nr)) then
292 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
293 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
294 if (rrin.gt.10000d0) then
295 rfraclow = dlog(rrhigh/rrin)/dlog(rrhigh/rrlow)
296 rfrachigh = dlog(rrin/rrlow)/dlog(rrhigh/rrlow)
299 if (ipart.eq.1.and.rrin.ge.rrr(1)) then
306 if (ipart.ne.1.and.rrin.ge.rrrg(1)) then
313 if (xxxx.ge.xx(261)) go to 245
315 nxlow = int(xxin/0.01) + 1
317 xfraclow = (xx(nxhigh)-xxin)/0.01
318 xfrachigh = (xxin - xx(nxlow))/0.01
321 clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh)
322 chigh = xfraclow*caq(nrhigh,nxlow)+xfrachigh*caq(nrhigh,nxhigh)
324 clow = xfraclow*cag(nrlow,nxlow)+xfrachigh*cag(nrlow,nxhigh)
325 chigh = xfraclow*cag(nrhigh,nxlow)+xfrachigh*cag(nrhigh,nxhigh)
328 continuous = rfraclow*clow + rfrachigh*chigh
333 discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh)
335 discrete = rfraclow*dag(nrlow) + rfrachigh*dag(nrhigh)
341 REAL*8 xxq(400), daq(34), caq(34,261), rrr(34)
342 COMMON /dataqua/ xxq, daq, caq, rrr
344 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
345 COMMON /dataglu/ xxg, dag, cag, rrrg
347 OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90)
349 read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn),
350 + caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn),
351 + caq(9,nn), caq(10,nn), caq(11,nn), caq(12,nn),
353 + caq(14,nn), caq(15,nn), caq(16,nn), caq(17,nn),
355 + caq(19,nn), caq(20,nn), caq(21,nn), caq(22,nn),
357 + caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn),
359 + caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn),
360 + caq(33,nn), caq(34,nn)
363 read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn),
364 + cag(4,nn), cag(5,nn), cag(6,nn), cag(7,nn), cag(8,nn),
365 + cag(9,nn), cag(10,nn), cag(11,nn), cag(12,nn),
367 + cag(14,nn), cag(15,nn), cag(16,nn), cag(17,nn),
369 + cag(19,nn), cag(20,nn), cag(21,nn), cag(22,nn),
371 + cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn),
373 + cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn),
374 + cag(33,nn), cag(34,nn)
378 OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91)
380 read (21,*) rrr(nn), daq(nn)
383 read (21,*) rrrg(nn), dag(nn)
388 90 PRINT*, 'input - output error'
389 91 PRINT*, 'input - output error #2'
395 =======================================================================
397 Adapted to ROOT macro by A. Dainese - 13/07/2003
398 Ported to class by C. Loizides - 12/02/2004
399 New version for extended R values added - 06/03/2004
402 Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
403 Double_t &continuous,Double_t &discrete) const
405 // Calculate Multiple Scattering approx.
406 // weights for given parton type,
407 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
413 //read-in data before first call
415 Error("CalcMult","Tables are not loaded.");
419 Error("CalcMult","Tables are not loaded for Multiple Scattering.");
423 Double_t rrin = rrrr;
424 Double_t xxin = xxxx;
426 if(xxin>fxx[260]) return -1;
427 Int_t nxlow = (Int_t)(xxin/0.01) + 1;
428 Int_t nxhigh = nxlow + 1;
429 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.01;
430 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01;
433 if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD
434 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
436 Int_t nrlow=0,nrhigh=0;
437 Double_t rrhigh=0,rrlow=0;
438 for(Int_t nr=1; nr<=34; nr++) {
439 if(rrin<frrr[nr-1]) {
442 rrhigh = frrr[nr-1-1];
452 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
453 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
456 rfraclow = TMath::Log2(rrhigh/rrin)/TMath::Log2(rrhigh/rrlow);
457 rfrachigh = TMath::Log2(rrin/rrlow)/TMath::Log2(rrhigh/rrlow);
459 if((ipart==1) && (rrin>=frrr[0]))
466 if((ipart==2) && (rrin>=frrrg[0]))
474 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
476 Double_t clow=0,chigh=0;
478 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
479 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
481 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
482 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
485 continuous = rfraclow*clow + rfrachigh*chigh;
486 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
487 //rfraclow,clow,rfrachigh,chigh,continuous);
490 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
492 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
498 Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *discall)
500 // read in tables for Single Hard Approx.
501 // path to continuum and to discrete part
503 fTablesLoaded = kFALSE;
507 sprintf(fname,"%s",gSystem->ExpandPathName(contall));
508 //PH ifstream fincont(fname);
509 fstream fincont(fname,ios::in);
510 #if defined(__HP_aCC) || defined(__DECCXX)
511 if(!fincont.rdbuf()->is_open()) return -1;
513 if(!fincont.is_open()) return -1;
517 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
518 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
519 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>
521 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>
523 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>
525 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>
534 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
535 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
536 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>
538 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>
540 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>
542 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>
550 sprintf(fname,"%s",gSystem->ExpandPathName(discall));
551 //PH ifstream findisc(fname);
552 fstream findisc(fname,ios::in);
553 #if defined(__HP_aCC) || defined(__DECCXX)
554 if(!findisc.rdbuf()->is_open()) return -1;
556 if(!findisc.is_open()) return -1;
560 while(findisc>>frrr[nn]>>fdaq[nn]) {
565 while(findisc>>frrrg[nn]>>fdag[nn]) {
571 fTablesLoaded = kTRUE;
576 C***************************************************************************
577 C Quenching Weights for Single Hard Scattering
582 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
584 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
587 C This package contains quenching weights for gluon radiation in the
588 C single hard scattering approximation.
590 C swqlin returns the quenching weight for a quark (ipart=1) or
591 C a gluon (ipart=2) traversing a medium with Debye screening mass mu and
592 C length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where
593 C wc=0.5*mu^2*L and w is the energy radiated. The output values are
594 C the continuous and discrete (prefactor of the delta function) parts
595 C of the quenching weights.
597 C In order to use this routine, the files contlin.all and disclin.all
598 C need to be in the working directory.
600 C An initialization of the tables is needed by doing call initlin before
603 C Please, send us any comment:
605 C urs.wiedemann@cern.ch
606 C carlos.salgado@cern.ch
609 C-------------------------------------------------------------------
612 SUBROUTINE swqlin(ipart,rrrr,xxxx,continuous,discrete)
614 REAL*8 xx(400), dalq(30), calq(30,261), rrr(30)
615 COMMON /datalinqua/ xx, dalq, calq, rrr
617 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
618 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
620 REAL*8 rrrr,xxxx, continuous, discrete
622 INTEGER nrlow, nrhigh, nxlow, nxhigh
623 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
624 REAL*8 xfraclow, xfrachigh
630 nxlow = int(xxin/0.038) + 1
632 xfraclow = (xx(nxhigh)-xxin)/0.038
633 xfrachigh = (xxin - xx(nxlow))/0.038
636 if (rrin.lt.rrr(nr)) then
648 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
649 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
652 clow = xfraclow*calq(nrlow,nxlow)+xfrachigh*calq(nrlow,nxhigh)
653 chigh = xfraclow*calq(nrhigh,nxlow)+xfrachigh*calq(nrhigh,nxhigh)
655 clow = xfraclow*calg(nrlow,nxlow)+xfrachigh*calg(nrlow,nxhigh)
656 chigh = xfraclow*calg(nrhigh,nxlow)+xfrachigh*calg(nrhigh,nxhigh)
659 continuous = rfraclow*clow + rfrachigh*chigh
662 discrete = rfraclow*dalq(nrlow) + rfrachigh*dalq(nrhigh)
664 discrete = rfraclow*dalg(nrlow) + rfrachigh*dalg(nrhigh)
670 REAL*8 xxlq(400), dalq(30), calq(30,261), rrr(30)
671 COMMON /datalinqua/ xxlq, dalq, calq, rrr
673 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
674 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
676 OPEN(UNIT=20,FILE='contlin.all',STATUS='OLD',ERR=90)
678 read (20,*) xxlq(nn), calq(1,nn), calq(2,nn), calq(3,nn),
679 + calq(4,nn), calq(5,nn), calq(6,nn), calq(7,nn), calq(8,nn),
680 + calq(9,nn), calq(10,nn), calq(11,nn), calq(12,nn),
682 + calq(14,nn), calq(15,nn), calq(16,nn), calq(17,nn),
684 + calq(19,nn), calq(20,nn), calq(21,nn), calq(22,nn),
686 + calq(24,nn), calq(25,nn), calq(26,nn), calq(27,nn),
688 + calq(29,nn), calq(30,nn)
691 read (20,*) xxlg(nn), calg(1,nn), calg(2,nn), calg(3,nn),
692 + calg(4,nn), calg(5,nn), calg(6,nn), calg(7,nn), calg(8,nn),
693 + calg(9,nn), calg(10,nn), calg(11,nn), calg(12,nn),
695 + calg(14,nn), calg(15,nn), calg(16,nn), calg(17,nn),
697 + calg(19,nn), calg(20,nn), calg(21,nn), calg(22,nn),
699 + calg(24,nn), calg(25,nn), calg(26,nn), calg(27,nn),
701 + calg(29,nn), calg(30,nn)
705 OPEN(UNIT=21,FILE='disclin.all',STATUS='OLD',ERR=91)
707 read (21,*) rrr(nn), dalq(nn)
710 read (21,*) rrrlg(nn), dalg(nn)
715 90 PRINT*, 'input - output error'
716 91 PRINT*, 'input - output error #2'
721 =======================================================================
723 Ported to class by C. Loizides - 17/02/2004
727 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx,
728 Double_t &continuous,Double_t &discrete) const
730 // calculate Single Hard approx.
731 // weights for given parton type,
732 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
734 // read-in data before first call
736 Error("CalcSingleHard","Tables are not loaded.");
740 Error("CalcSingleHard","Tables are not loaded for Single Hard Scattering.");
744 Double_t rrin = rrrr;
745 Double_t xxin = xxxx;
747 Int_t nxlow = (Int_t)(xxin/0.038) + 1;
748 Int_t nxhigh = nxlow + 1;
749 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.038;
750 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.038;
753 if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD
754 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
756 Int_t nrlow=0,nrhigh=0;
757 Double_t rrhigh=0,rrlow=0;
758 for(Int_t nr=1; nr<=30; nr++) {
759 if(rrin<frrr[nr-1]) {
762 rrhigh = frrr[nr-1-1];
772 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
773 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
775 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
777 Double_t clow=0,chigh=0;
779 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
780 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
782 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
783 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
786 continuous = rfraclow*clow + rfrachigh*chigh;
787 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
788 // rfraclow,clow,rfrachigh,chigh,continuous);
791 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
793 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
799 Int_t AliQuenchingWeights::CalcMult(Int_t ipart,
800 Double_t w,Double_t qtransp,Double_t length,
801 Double_t &continuous,Double_t &discrete) const
803 // Calculate Multiple Scattering approx.
804 // weights for given parton type,
805 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
807 Double_t wc=CalcWC(qtransp,length);
808 Double_t rrrr=CalcR(wc,length);
810 return CalcMult(ipart,rrrr,xxxx,continuous,discrete);
813 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
814 Double_t w,Double_t mu,Double_t length,
815 Double_t &continuous,Double_t &discrete) const
817 // calculate Single Hard approx.
818 // weights for given parton type,
819 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
821 Double_t wcbar=CalcWCbar(mu,length);
822 Double_t rrrr=CalcR(wcbar,length);
823 Double_t xxxx=w/wcbar;
824 return CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
827 Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const
829 //calculate r value and
830 //check if it is less then maximum
832 Double_t r = wc*l*fgkConvFmToInvGeV;
834 Warning("CalcR","Value of r = %.2f; should be less than %.2f", r, fgkRMax);
840 Double_t AliQuenchingWeights::CalcRk(Double_t k, Double_t I0, Double_t I1) const
842 //calculate R value and
843 //check if it is less then maximum
845 Double_t r = fgkRMax-1;
849 Warning("CalcRk","Value of r = %.2f; should be less than %.2f",r,fgkRMax);
855 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const
857 // return DeltaE for MS or SH scattering
858 // for given parton type, length and energy
859 // Dependant on ECM (energy constraint method)
860 // e is used to determine where to set bins to zero.
863 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
866 if((ipart<1) || (ipart>2)) {
867 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
871 Int_t l=GetIndex(length);
874 if(fECMethod==kReweight){
878 ret=fHistos[ipart-1][l-1]->GetRandom();
880 Warning("GetELossRandom",
881 "Aborted reweighting; maximum loss assigned after 1e6 trials.");
888 Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
893 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const
895 //return quenched parton energy
896 //for given parton type, length and energy
898 Double_t loss=GetELossRandom(ipart,length,e);
902 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const
904 // return DeltaE for MS or SH scattering
905 // for given parton type, length distribution and energy
906 // Dependant on ECM (energy constraint method)
907 // e is used to determine where to set bins to zero.
910 Warning("GetELossRandom","Pointer to length distribution is NULL.");
913 Double_t ell=hell->GetRandom();
914 return GetELossRandom(ipart,ell,e);
917 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e) const
919 //return quenched parton energy
920 //for given parton type, length distribution and energy
922 Double_t loss=GetELossRandom(ipart,hell,e);
926 Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
928 // return DeltaE for new dynamic version
929 // for given parton type, I0 and I1 value and energy
930 // Dependant on ECM (energy constraint method)
931 // e is used to determine where to set bins to zero.
933 // read-in data before first call
935 Fatal("GetELossRandomK","Tables are not loaded.");
938 if((ipart<1) || (ipart>2)) {
939 Fatal("GetELossRandomK","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
943 Double_t r=CalcRk(I0,I1);
945 Fatal("GetELossRandomK","R should not be negative");
948 Double_t wc=CalcWCk(I1);
950 Fatal("GetELossRandomK","wc should be greater than zero");
953 if(SampleEnergyLoss(ipart,r)!=0){
954 Fatal("GetELossRandomK","Could not sample energy loss");
958 if(fECMethod==kReweight){
962 ret=fHisto->GetRandom();
964 Warning("GetELossRandomK",
965 "Aborted reweighting; maximum loss assigned after 1e6 trials.");
973 Double_t ret=fHisto->GetRandom()*wc;
978 Double_t AliQuenchingWeights::CalcQuenchedEnergyK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
980 //return quenched parton energy
981 //for given parton type, I0 and I1 value and energy
983 Double_t loss=GetELossRandomK(ipart,I0,I1,e);
987 Double_t AliQuenchingWeights::GetELossRandomKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
989 // return DeltaE for new dynamic version
990 // for given parton type, I0 and I1 value and energy
991 // Dependant on ECM (energy constraint method)
992 // e is used to determine where to set bins to zero.
993 // method is optimized and should only be used if
994 // all parameters are well within the bounds.
995 // read-in data tables before first call
997 Double_t r=CalcRk(I0,I1);
1002 Double_t wc=CalcWCk(I1);
1007 return GetELossRandomKFastR(ipart,r,wc,e);
1010 Double_t AliQuenchingWeights::GetELossRandomKFastR(Int_t ipart, Double_t r, Double_t wc, Double_t e)
1012 // return DeltaE for new dynamic version
1013 // for given parton type, R and wc value and energy
1014 // Dependant on ECM (energy constraint method)
1015 // e is used to determine where to set bins to zero.
1016 // method is optimized and should only be used if
1017 // all parameters are well within the bounds.
1018 // read-in data tables before first call
1024 Double_t discrete=0.;
1025 Double_t continuous=0.;
1027 Double_t xxxx = fHisto->GetBinCenter(bin);
1029 CalcMult(ipart,r,xxxx,continuous,discrete);
1031 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1034 return 0.; //no energy loss
1037 fHisto->SetBinContent(bin,continuous);
1038 Int_t kbinmax=fHisto->FindBin(e/wc);
1039 if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
1040 if(kbinmax==1) return e; //maximum energy loss
1043 for(bin=2; bin<=kbinmax; bin++) {
1044 xxxx = fHisto->GetBinCenter(bin);
1045 CalcMult(ipart,r,xxxx,continuous,discrete);
1046 fHisto->SetBinContent(bin,continuous);
1049 for(bin=2; bin<=kbinmax; bin++) {
1050 xxxx = fHisto->GetBinCenter(bin);
1051 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1052 fHisto->SetBinContent(bin,continuous);
1056 if(fECMethod==kReweight){
1057 fHisto->SetBinContent(kbinmax+1,0);
1058 fHisto->Fill(0.,discrete*fgkBins/fgkMaxBin);
1059 } else if (fECMethod==kReweightCont) {
1060 fHisto->SetBinContent(kbinmax+1,0);
1061 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1062 fHisto->Scale(1./kdelta*(1-discrete));
1063 fHisto->Fill(0.,discrete);
1065 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1066 Double_t val=discrete*fgkBins/fgkMaxBin;
1067 fHisto->Fill(0.,val);
1068 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1070 for(bin=kbinmax+2; bin<=fgkBins; bin++) {
1071 fHisto->SetBinContent(bin,0);
1073 //cout << kbinmax << " " << discrete << " " << fHisto->Integral() << endl;
1074 Double_t ret=fHisto->GetRandom()*wc;
1079 Double_t AliQuenchingWeights::CalcQuenchedEnergyKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
1081 //return quenched parton energy (fast method)
1082 //for given parton type, I0 and I1 value and energy
1084 Double_t loss=GetELossRandomKFast(ipart,I0,I1,e);
1088 Double_t AliQuenchingWeights::GetDiscreteWeight(Int_t ipart, Double_t I0, Double_t I1)
1090 // return discrete weight
1092 Double_t r=CalcRk(I0,I1);
1096 return GetDiscreteWeightR(ipart,r);
1099 Double_t AliQuenchingWeights::GetDiscreteWeightR(Int_t ipart, Double_t r)
1101 // return discrete weight
1107 Double_t discrete=0.;
1108 Double_t continuous=0.;
1110 Double_t xxxx = fHisto->GetBinCenter(bin);
1112 CalcMult(ipart,r,xxxx,continuous,discrete);
1114 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1118 void AliQuenchingWeights::GetZeroLossProb(Double_t &p,Double_t &prw,Double_t &prwcont,
1119 Int_t ipart,Double_t I0,Double_t I1,Double_t e)
1121 //calculate the probabilty that there is no energy
1122 //loss for different ways of energy constraint
1123 p=1.;prw=1.;prwcont=1.;
1124 Double_t r=CalcRk(I0,I1);
1128 Double_t wc=CalcWCk(I1);
1132 GetZeroLossProbR(p,prw,prwcont,ipart,r,wc,e);
1135 void AliQuenchingWeights::GetZeroLossProbR(Double_t &p,Double_t &prw,Double_t &prwcont,
1136 Int_t ipart, Double_t r,Double_t wc,Double_t e)
1138 //calculate the probabilty that there is no energy
1139 //loss for different ways of energy constraint
1144 Double_t discrete=0.;
1145 Double_t continuous=0.;
1147 Int_t kbinmax=fHisto->FindBin(e/wc);
1148 if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
1150 for(Int_t bin=1; bin<=kbinmax; bin++) {
1151 Double_t xxxx = fHisto->GetBinCenter(bin);
1152 CalcMult(ipart,r,xxxx,continuous,discrete);
1153 fHisto->SetBinContent(bin,continuous);
1156 for(Int_t bin=1; bin<=kbinmax; bin++) {
1157 Double_t xxxx = fHisto->GetBinCenter(bin);
1158 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1159 fHisto->SetBinContent(bin,continuous);
1163 //non-reweighted P(Delta E = 0)
1164 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1165 Double_t val=discrete*fgkBins/fgkMaxBin;
1166 fHisto->Fill(0.,val);
1167 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1168 Double_t hint=fHisto->Integral(1,kbinmax+1);
1169 p=fHisto->GetBinContent(1)/hint;
1172 hint=fHisto->Integral(1,kbinmax);
1173 prw=fHisto->GetBinContent(1)/hint;
1175 Double_t xxxx = fHisto->GetBinCenter(1);
1176 CalcMult(ipart,r,xxxx,continuous,discrete);
1177 fHisto->SetBinContent(1,continuous);
1178 hint=fHisto->Integral(1,kbinmax);
1179 fHisto->Scale(1./hint*(1-discrete));
1180 fHisto->Fill(0.,discrete);
1181 prwcont=fHisto->GetBinContent(1);
1184 Int_t AliQuenchingWeights::SampleEnergyLoss()
1186 // Has to be called to fill the histograms
1188 // For stored values fQTransport loop over
1189 // particle type and length = 1 to fMaxLength (fm)
1190 // to fill energy loss histos
1192 // Take histogram of continuous weights
1193 // Take discrete_weight
1194 // If discrete_weight > 1, put all channels to 0, except channel 1
1195 // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral
1197 // read-in data before first call
1199 Error("SampleEnergyLoss","Tables are not loaded.");
1204 Int_t lmax=CalcLengthMax(fQTransport);
1205 if(fLengthMax>lmax){
1206 Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
1210 Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested.");
1214 fHistos=new TH1F**[2];
1215 fHistos[0]=new TH1F*[4*fLengthMax];
1216 fHistos[1]=new TH1F*[4*fLengthMax];
1217 fLengthMaxOld=fLengthMax; //remember old value in case
1218 //user wants to reset
1221 Char_t meddesc[100];
1223 medvalue=(Int_t)(fQTransport*1000.);
1224 sprintf(meddesc,"MS");
1226 medvalue=(Int_t)(fMu*1000.);
1227 sprintf(meddesc,"SH");
1230 for(Int_t ipart=1;ipart<=2;ipart++){
1231 for(Int_t l=1;l<=4*fLengthMax;l++){
1233 sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
1235 Double_t wc = CalcWC(len);
1236 fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgkBins,0.,fgkMaxBin*wc);
1237 fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
1238 fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
1239 fHistos[ipart-1][l-1]->SetLineColor(4);
1241 Double_t rrrr = CalcR(wc,len);
1242 Double_t discrete=0.;
1243 // loop on histogram channels
1244 for(Int_t bin=1; bin<=fgkBins; bin++) {
1245 Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
1246 Double_t continuous;
1248 CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1250 CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1251 fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
1253 // add discrete part to distribution
1255 for(Int_t bin=2;bin<=fgkBins;bin++)
1256 fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
1258 Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1259 fHistos[ipart-1][l-1]->Fill(0.,val);
1261 Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1262 fHistos[ipart-1][l-1]->Scale(1./hint);
1268 Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t r)
1270 // Sample energy loss directly for one particle type
1271 // choses R (safe it and keep it until next call of function)
1273 // read-in data before first call
1275 Error("SampleEnergyLoss","Tables are not loaded.");
1279 Double_t discrete=0.;
1280 Double_t continuous=0;;
1282 Double_t xxxx = fHisto->GetBinCenter(bin);
1284 CalcMult(ipart,r,xxxx,continuous,discrete);
1286 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1289 fHisto->SetBinContent(1,1.);
1290 for(bin=2;bin<=fgkBins;bin++)
1291 fHisto->SetBinContent(bin,0.);
1295 fHisto->SetBinContent(bin,continuous);
1296 for(bin=2; bin<=fgkBins; bin++) {
1297 xxxx = fHisto->GetBinCenter(bin);
1299 CalcMult(ipart,r,xxxx,continuous,discrete);
1301 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1302 fHisto->SetBinContent(bin,continuous);
1305 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
1306 fHisto->Fill(0.,val);
1307 Double_t hint=fHisto->Integral(1,fgkBins);
1309 fHisto->Scale(1./hint);
1311 //cout << discrete << " " << hint << " " << continuous << endl;
1317 const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) const
1319 //return quenching histograms
1320 //for ipart and length
1323 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
1326 if((ipart<1) || (ipart>2)) {
1327 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
1331 Int_t l=GetIndex(length);
1333 return fHistos[ipart-1][l-1];
1336 TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const
1338 // ipart = 1 for quark, 2 for gluon
1339 // medval a) qtransp = transport coefficient (GeV^2/fm)
1340 // b) mu = Debye mass (GeV)
1341 // length = path length in medium (fm)
1342 // Get from SW tables:
1343 // - continuous weight, as a function of dE/wc
1346 Char_t meddesc[100];
1348 wc=CalcWC(medval,length);
1349 sprintf(meddesc,"MS");
1351 wc=CalcWCbar(medval,length);
1352 sprintf(meddesc,"SH");
1356 sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart,
1357 (Int_t)(medval*1000.),(Int_t)length);
1359 TH1F *hist = new TH1F("hist",hname,fgkBins,0.,fgkMaxBin*wc);
1360 hist->SetXTitle("#Delta E [GeV]");
1361 hist->SetYTitle("p(#Delta E)");
1362 hist->SetLineColor(4);
1364 Double_t rrrr = CalcR(wc,length);
1365 //loop on histogram channels
1366 for(Int_t bin=1; bin<=fgkBins; bin++) {
1367 Double_t xxxx = hist->GetBinCenter(bin)/wc;
1368 Double_t continuous,discrete;
1370 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1371 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1376 hist->SetBinContent(bin,continuous);
1381 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const
1383 // ipart = 1 for quark, 2 for gluon
1384 // medval a) qtransp = transport coefficient (GeV^2/fm)
1385 // b) mu = Debye mass (GeV)
1386 // length = path length in medium (fm)
1387 // Get from SW tables:
1388 // - continuous weight, as a function of dE/wc
1391 Char_t meddesc[100];
1393 wc=CalcWC(medval,length);
1394 sprintf(meddesc,"MS");
1396 wc=CalcWCbar(medval,length);
1397 sprintf(meddesc,"SH");
1401 sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart,
1402 (Int_t)(medval*1000.),(Int_t)length);
1404 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
1405 histx->SetXTitle("x = #Delta E/#omega_{c}");
1407 histx->SetYTitle("p(#Delta E/#omega_{c})");
1409 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1410 histx->SetLineColor(4);
1412 Double_t rrrr = CalcR(wc,length);
1413 //loop on histogram channels
1414 for(Int_t bin=1; bin<=fgkBins; bin++) {
1415 Double_t xxxx = histx->GetBinCenter(bin);
1416 Double_t continuous,discrete;
1418 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1419 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1424 histx->SetBinContent(bin,continuous);
1429 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t r) const
1431 // compute P(E) distribution for
1432 // given ipart = 1 for quark, 2 for gluon
1435 Char_t meddesc[100];
1437 sprintf(meddesc,"MS");
1439 sprintf(meddesc,"SH");
1443 sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,r);
1444 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
1445 histx->SetXTitle("x = #Delta E/#omega_{c}");
1447 histx->SetYTitle("p(#Delta E/#omega_{c})");
1449 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1450 histx->SetLineColor(4);
1453 Double_t continuous=0.,discrete=0.;
1454 //loop on histogram channels
1455 for(Int_t bin=1; bin<=fgkBins; bin++) {
1456 Double_t xxxx = histx->GetBinCenter(bin);
1458 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1459 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1464 histx->SetBinContent(bin,continuous);
1467 //add discrete part to distribution
1469 for(Int_t bin=2;bin<=fgkBins;bin++)
1470 histx->SetBinContent(bin,0.);
1472 Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins);
1473 histx->Fill(0.,val);
1475 Double_t hint=histx->Integral(1,fgkBins);
1476 if(hint!=0) histx->Scale(1./hint);
1481 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
1483 // compute energy loss histogram for
1484 // parton type, medium value, length and energy
1486 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1488 dummy->SetQTransport(medval);
1491 dummy->SetMu(medval);
1492 dummy->InitSingleHard();
1494 dummy->SampleEnergyLoss();
1499 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1500 sprintf(hname,"hLossQuarks");
1502 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1503 sprintf(hname,"hLossGluons");
1506 TH1F *h = new TH1F(hname,name,250,0,250);
1507 for(Int_t i=0;i<100000;i++){
1508 //if(i % 1000 == 0) cout << "." << flush;
1509 Double_t loss=dummy->GetELossRandom(ipart,l,e);
1517 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
1519 // compute energy loss histogram for
1520 // parton type, medium value,
1521 // length distribution and energy
1523 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1525 dummy->SetQTransport(medval);
1528 dummy->SetMu(medval);
1529 dummy->InitSingleHard();
1531 dummy->SampleEnergyLoss();
1536 sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1537 sprintf(hname,"hLossQuarks");
1539 sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1540 sprintf(hname,"hLossGluons");
1543 TH1F *h = new TH1F(hname,name,250,0,250);
1544 for(Int_t i=0;i<100000;i++){
1545 //if(i % 1000 == 0) cout << "." << flush;
1546 Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
1554 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t r) const
1556 // compute energy loss histogram for
1557 // parton type and given R
1559 TH1F *dummy = ComputeQWHistoX(ipart,r);
1560 if(!dummy) return 0;
1563 sprintf(hname,"hELossHistox_%d_%.2f",ipart,r);
1564 TH1F *histx = new TH1F("histxr",hname,fgkBins,0.,fgkMaxBin);
1565 for(Int_t i=0;i<100000;i++){
1566 //if(i % 1000 == 0) cout << "." << flush;
1567 Double_t loss=dummy->GetRandom();
1574 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,Double_t l) const
1576 // compute average energy loss for
1577 // parton type, medium value, length and energy
1579 TH1F *dummy = ComputeELossHisto(ipart,medval,l);
1580 if(!dummy) return 0;
1581 Double_t ret=dummy->GetMean();
1586 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEll) const
1588 // compute average energy loss for
1589 // parton type, medium value,
1590 // length distribution and energy
1592 TH1F *dummy = ComputeELossHisto(ipart,medval,hEll);
1593 if(!dummy) return 0;
1594 Double_t ret=dummy->GetMean();
1599 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t r) const
1601 // compute average energy loss over wc
1602 // for parton type and given R
1604 TH1F *dummy = ComputeELossHisto(ipart,r);
1605 if(!dummy) return 0;
1606 Double_t ret=dummy->GetMean();
1611 void AliQuenchingWeights::PlotDiscreteWeights(Double_t len,Double_t qm) const
1613 // plot discrete weights for given length
1617 c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
1619 c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
1622 TH2F *hframe = new TH2F("hdisc","",2,0,qm+.1,2,0,1.25);
1623 hframe->SetStats(0);
1625 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1627 hframe->SetXTitle("#mu [GeV]");
1628 //hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
1629 hframe->SetYTitle("p_{0} (discrete weight)");
1632 Int_t points=(Int_t)qm*4;
1633 TGraph *gq=new TGraph(points);
1636 for(Double_t q=0.05;q<=qm+.05;q+=0.25){
1638 CalcMult(1,1.0,q,len,cont,disc);
1639 gq->SetPoint(i,q,disc);i++;
1642 for(Double_t m=0.05;m<=qm+.05;m+=0.25){
1644 CalcSingleHard(1,1.0,m,len,cont, disc);
1645 gq->SetPoint(i,m,disc);i++;
1648 gq->SetMarkerStyle(20);
1649 gq->SetMarkerColor(1);
1650 gq->SetLineStyle(1);
1651 gq->SetLineColor(1);
1654 TGraph *gg=new TGraph(points);
1657 for(Double_t q=0.05;q<=qm+.05;q+=0.25){
1659 CalcMult(2,1.0,q,len,cont,disc);
1660 gg->SetPoint(i,q,disc);i++;
1663 for(Double_t m=0.05;m<=qm+.05;m+=0.25){
1665 CalcSingleHard(2,1.0,m,len,cont,disc);
1666 gg->SetPoint(i,m,disc);i++;
1669 gg->SetMarkerStyle(24);
1670 gg->SetMarkerColor(2);
1671 gg->SetLineStyle(2);
1672 gg->SetLineColor(2);
1675 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1676 l1a->SetFillStyle(0);
1677 l1a->SetBorderSize(0);
1679 sprintf(label,"L = %.1f fm",len);
1680 l1a->AddEntry(gq,label,"");
1681 l1a->AddEntry(gq,"quark projectile","l");
1682 l1a->AddEntry(gg,"gluon projectile","l");
1688 void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
1690 // plot continous weights for
1691 // given parton type and length
1698 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1700 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1702 medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1703 sprintf(name,"ccont-ms-%d",itype);
1706 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1708 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1710 medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1711 sprintf(name,"ccont-ms-%d",itype);
1714 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1716 TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1718 h1->SetTitle(title);
1720 h1->SetLineColor(1);
1722 TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1724 h2->SetLineColor(2);
1725 h2->DrawCopy("SAME");
1726 TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1728 h3->SetLineColor(3);
1729 h3->DrawCopy("SAME");
1731 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1732 l1a->SetFillStyle(0);
1733 l1a->SetBorderSize(0);
1735 sprintf(label,"L = %.1f fm",ell);
1736 l1a->AddEntry(h1,label,"");
1738 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1739 l1a->AddEntry(h1,label,"pl");
1740 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1741 l1a->AddEntry(h2,label,"pl");
1742 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1743 l1a->AddEntry(h3,label,"pl");
1745 sprintf(label,"#mu = %.1f GeV",medvals[0]);
1746 l1a->AddEntry(h1,label,"pl");
1747 sprintf(label,"#mu = %.1f GeV",medvals[1]);
1748 l1a->AddEntry(h2,label,"pl");
1749 sprintf(label,"#mu = %.1f GeV",medvals[2]);
1750 l1a->AddEntry(h3,label,"pl");
1757 void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
1759 // plot continous weights for
1760 // given parton type and medium value
1766 sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
1768 sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
1770 sprintf(name,"ccont2-ms-%d",itype);
1773 sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
1775 sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
1777 sprintf(name,"ccont2-sh-%d",itype);
1779 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1781 TH1F *h1=ComputeQWHisto(itype,medval,8);
1783 h1->SetTitle(title);
1785 h1->SetLineColor(1);
1787 TH1F *h2=ComputeQWHisto(itype,medval,5);
1789 h2->SetLineColor(2);
1790 h2->DrawCopy("SAME");
1791 TH1F *h3=ComputeQWHisto(itype,medval,2);
1793 h3->SetLineColor(3);
1794 h3->DrawCopy("SAME");
1796 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1797 l1a->SetFillStyle(0);
1798 l1a->SetBorderSize(0);
1801 sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval);
1803 sprintf(label,"#mu = %.1f GeV",medval);
1805 l1a->AddEntry(h1,label,"");
1806 l1a->AddEntry(h1,"L = 8 fm","pl");
1807 l1a->AddEntry(h2,"L = 5 fm","pl");
1808 l1a->AddEntry(h3,"L = 2 fm","pl");
1814 void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t qm,Double_t e) const
1816 // plot average energy loss for given length
1817 // and parton energy
1820 Error("PlotAvgELoss","Tables are not loaded.");
1827 sprintf(title,"Average Energy Loss for Multiple Scattering");
1828 sprintf(name,"cavgelossms");
1830 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1831 sprintf(name,"cavgelosssh");
1834 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1836 TH2F *hframe = new TH2F("avgloss","",2,0,qm+.1,2,0,100);
1837 hframe->SetStats(0);
1839 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1841 hframe->SetXTitle("#mu [GeV]");
1842 hframe->SetYTitle("<E_{loss}> [GeV]");
1845 TGraph *gq=new TGraph(20);
1847 for(Double_t v=0.05;v<=qm+.05;v+=0.25){
1848 TH1F *dummy=ComputeELossHisto(1,v,len,e);
1849 Double_t avgloss=dummy->GetMean();
1850 gq->SetPoint(i,v,avgloss);i++;
1853 gq->SetMarkerStyle(21);
1856 Int_t points=(Int_t)qm*4;
1857 TGraph *gg=new TGraph(points);
1859 for(Double_t v=0.05;v<=qm+.05;v+=0.25){
1860 TH1F *dummy=ComputeELossHisto(2,v,len,e);
1861 Double_t avgloss=dummy->GetMean();
1862 gg->SetPoint(i,v,avgloss);i++;
1865 gg->SetMarkerStyle(20);
1866 gg->SetMarkerColor(2);
1869 TGraph *gratio=new TGraph(points);
1870 for(i=0;i<points;i++){
1872 gg->GetPoint(i,x,y);
1873 gq->GetPoint(i,x2,y2);
1875 gratio->SetPoint(i,x,y/y2*10/2.25);
1876 else gratio->SetPoint(i,x,0);
1878 gratio->SetLineStyle(4);
1880 TLegend *l1a = new TLegend(0.15,0.60,0.50,0.90);
1881 l1a->SetFillStyle(0);
1882 l1a->SetBorderSize(0);
1884 sprintf(label,"L = %.1f fm",len);
1885 l1a->AddEntry(gq,label,"");
1886 l1a->AddEntry(gq,"quark projectile","pl");
1887 l1a->AddEntry(gg,"gluon projectile","pl");
1888 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1894 void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
1896 // plot average energy loss for given
1897 // length distribution and parton energy
1900 Error("PlotAvgELossVs","Tables are not loaded.");
1907 sprintf(title,"Average Energy Loss for Multiple Scattering");
1908 sprintf(name,"cavgelossms2");
1910 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1911 sprintf(name,"cavgelosssh2");
1914 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1916 TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1917 hframe->SetStats(0);
1919 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1921 hframe->SetXTitle("#mu [GeV]");
1922 hframe->SetYTitle("<E_{loss}> [GeV]");
1925 TGraph *gq=new TGraph(20);
1927 for(Double_t v=0.05;v<=5.05;v+=0.25){
1928 TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1929 Double_t avgloss=dummy->GetMean();
1930 gq->SetPoint(i,v,avgloss);i++;
1933 gq->SetMarkerStyle(20);
1936 TGraph *gg=new TGraph(20);
1938 for(Double_t v=0.05;v<=5.05;v+=0.25){
1939 TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1940 Double_t avgloss=dummy->GetMean();
1941 gg->SetPoint(i,v,avgloss);i++;
1944 gg->SetMarkerStyle(24);
1947 TGraph *gratio=new TGraph(20);
1950 gg->GetPoint(i,x,y);
1951 gq->GetPoint(i,x2,y2);
1953 gratio->SetPoint(i,x,y/y2*10/2.25);
1954 else gratio->SetPoint(i,x,0);
1956 gratio->SetLineStyle(4);
1959 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1960 l1a->SetFillStyle(0);
1961 l1a->SetBorderSize(0);
1963 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
1964 l1a->AddEntry(gq,label,"");
1965 l1a->AddEntry(gq,"quark","pl");
1966 l1a->AddEntry(gg,"gluon","pl");
1967 //l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1973 void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const
1975 // plot average energy loss versus ell
1978 Error("PlotAvgELossVsEll","Tables are not loaded.");
1986 sprintf(title,"Average Energy Loss for Multiple Scattering");
1987 sprintf(name,"cavgelosslms");
1990 sprintf(title,"Average Energy Loss for Single Hard Scattering");
1991 sprintf(name,"cavgelosslsh");
1995 TCanvas *c = new TCanvas(name,title,0,0,600,400);
1997 TH2F *hframe = new TH2F("avglossell",title,2,0,fLengthMax,2,0,250);
1998 hframe->SetStats(0);
1999 hframe->SetXTitle("length [fm]");
2000 hframe->SetYTitle("<E_{loss}> [GeV]");
2003 TGraph *gq=new TGraph((Int_t)fLengthMax*4);
2005 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
2006 TH1F *dummy=ComputeELossHisto(1,medval,v,e);
2007 Double_t avgloss=dummy->GetMean();
2008 gq->SetPoint(i,v,avgloss);i++;
2011 gq->SetMarkerStyle(20);
2014 TGraph *gg=new TGraph((Int_t)fLengthMax*4);
2016 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
2017 TH1F *dummy=ComputeELossHisto(2,medval,v,e);
2018 Double_t avgloss=dummy->GetMean();
2019 gg->SetPoint(i,v,avgloss);i++;
2022 gg->SetMarkerStyle(24);
2025 TGraph *gratio=new TGraph((Int_t)fLengthMax*4);
2026 for(i=0;i<=(Int_t)fLengthMax*4;i++){
2028 gg->GetPoint(i,x,y);
2029 gq->GetPoint(i,x2,y2);
2031 gratio->SetPoint(i,x,y/y2*100/2.25);
2032 else gratio->SetPoint(i,x,0);
2034 gratio->SetLineStyle(1);
2035 gratio->SetLineWidth(2);
2037 TLegend *l1a = new TLegend(0.15,0.65,.60,0.85);
2038 l1a->SetFillStyle(0);
2039 l1a->SetBorderSize(0);
2042 sprintf(label,"#hat{q} = %.2f GeV^{2}/fm",medval);
2044 sprintf(label,"#mu = %.2f GeV",medval);
2045 l1a->AddEntry(gq,label,"");
2046 l1a->AddEntry(gq,"quark","pl");
2047 l1a->AddEntry(gg,"gluon","pl");
2048 l1a->AddEntry(gratio,"gluon/quark/2.25*100","pl");
2051 TF1 *f=new TF1("f","100",0,fLengthMax);
2058 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
2060 // plot relative energy loss for given
2061 // length and parton energy versus pt
2064 Error("PlotAvgELossVsPt","Tables are not loaded.");
2071 sprintf(title,"Relative Energy Loss for Multiple Scattering");
2072 sprintf(name,"cavgelossvsptms");
2074 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
2075 sprintf(name,"cavgelossvsptsh");
2078 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2080 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2081 hframe->SetStats(0);
2082 hframe->SetXTitle("p_{T} [GeV]");
2083 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2086 TGraph *gq=new TGraph(40);
2088 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2089 TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
2090 Double_t avgloss=dummy->GetMean();
2091 gq->SetPoint(i,pt,avgloss/pt);i++;
2094 gq->SetMarkerStyle(20);
2097 TGraph *gg=new TGraph(40);
2099 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2100 TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
2101 Double_t avgloss=dummy->GetMean();
2102 gg->SetPoint(i,pt,avgloss/pt);i++;
2105 gg->SetMarkerStyle(24);
2108 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2109 l1a->SetFillStyle(0);
2110 l1a->SetBorderSize(0);
2112 sprintf(label,"L = %.1f fm",len);
2113 l1a->AddEntry(gq,label,"");
2114 l1a->AddEntry(gq,"quark","pl");
2115 l1a->AddEntry(gg,"gluon","pl");
2121 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
2123 // plot relative energy loss for given
2124 // length distribution and parton energy versus pt
2127 Error("PlotAvgELossVsPt","Tables are not loaded.");
2134 sprintf(title,"Relative Energy Loss for Multiple Scattering");
2135 sprintf(name,"cavgelossvsptms2");
2137 sprintf(title,"Relative Energy Loss for Single Hard Scattering");
2138 sprintf(name,"cavgelossvsptsh2");
2140 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2142 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2143 hframe->SetStats(0);
2144 hframe->SetXTitle("p_{T} [GeV]");
2145 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2148 TGraph *gq=new TGraph(40);
2150 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2151 TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
2152 Double_t avgloss=dummy->GetMean();
2153 gq->SetPoint(i,pt,avgloss/pt);i++;
2156 gq->SetMarkerStyle(20);
2159 TGraph *gg=new TGraph(40);
2161 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2162 TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
2163 Double_t avgloss=dummy->GetMean();
2164 gg->SetPoint(i,pt,avgloss/pt);i++;
2167 gg->SetMarkerStyle(24);
2170 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2171 l1a->SetFillStyle(0);
2172 l1a->SetBorderSize(0);
2174 sprintf(label,"<L> = %.2f fm",hEll->GetMean());
2175 l1a->AddEntry(gq,label,"");
2176 l1a->AddEntry(gq,"quark","pl");
2177 l1a->AddEntry(gg,"gluon","pl");
2183 Int_t AliQuenchingWeights::GetIndex(Double_t len) const
2185 //get the index according to length
2186 if(len>fLengthMax) len=fLengthMax;
2188 Int_t l=Int_t(len/0.25);
2189 if((len-l*0.25)>0.125) l++;