1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
19 // Implementation of the class to calculate the parton energy loss
20 // Based on the "BDMPS" quenching weights by C.A.Salgado and U.A.Wiedemann
22 // C.A.Salgado and U.A.Wiedemann, Phys.Rev.D68 (2003) 014008 [hep-ph/0302184]
23 // A.Dainese, Eur.Phys.J.C, in press, [nucl-ex/0312005]
26 // Origin: C. Loizides constantinos.loizides@cern.ch
27 // A. Dainese andrea.dainese@pd.infn.it
29 //=================== Added by C. Loizides 27/03/04 ===========================
31 // Added support for k-Quenching, where wc=I1*k and R=2I1^2/I0*k
32 // (see the AliFastGlauber class for definition of I0/I1)
33 //-----------------------------------------------------------------------------
35 #include <Riostream.h>
45 #include "AliQuenchingWeights.h"
47 ClassImp(AliQuenchingWeights)
49 // conversion from fm to GeV^-1: 1 fm = fmGeV GeV^-1
50 const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197;
53 const Double_t AliQuenchingWeights::fgkRMax = 1.e6;
56 const Int_t AliQuenchingWeights::fgkBins = 1300;
57 const Double_t AliQuenchingWeights::fgkMaxBin = 1.3;
59 // counter for histogram labels
60 Int_t AliQuenchingWeights::fgCounter = 0;
63 AliQuenchingWeights::AliQuenchingWeights()
65 fInstanceNumber(fgCounter++),
81 AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a)
83 fInstanceNumber(fgCounter++),
100 fMultSoft=a.GetMultSoft();;
103 fQTransport=a.GetQTransport();
104 fECMethod=(kECMethod)a.GetECMethod();
105 fLengthMax=a.GetLengthMax();
106 fInstanceNumber=fgCounter++;
108 snprintf(name,100, "hhistoqw_%d",fInstanceNumber);
109 fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin);
110 for(Int_t bin=1;bin<=fgkBins;bin++)
111 fHisto->SetBinContent(bin,0.);
113 //Missing in the class is the pathname
114 //to the tables, can be added if needed
117 AliQuenchingWeights::~AliQuenchingWeights()
123 void AliQuenchingWeights::Init()
127 fHisto = new TH1F(Form("hhistoqw_%d",fInstanceNumber),"",fgkBins,0.,fgkMaxBin);
128 for(Int_t bin=1;bin<=fgkBins;bin++)
129 fHisto->SetBinContent(bin,0.);
132 void AliQuenchingWeights::Reset()
134 //reset tables if there were used
137 for(Int_t l=0;l<4*fLengthMaxOld;l++){
138 delete fHistos[0][l];
139 delete fHistos[1][l];
146 void AliQuenchingWeights::SetECMethod(kECMethod type)
148 //set energy constraint method
151 if(fECMethod==kDefault)
152 Info("SetECMethod","Energy Constraint Method set to DEFAULT:\nIf (sampled energy loss > parton energy) then sampled energy loss = parton energy.");
153 else if(fECMethod==kReweight)
154 Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy.");
155 else Info("SetECMethod","Energy Constraint Method set to REWEIGHTCONT:\nRequire sampled energy loss <= parton energy (only implemented for FAST method.");
158 Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
160 // read in tables for multiple scattering approximation
161 // path to continuum and to discrete part
163 fTablesLoaded = kFALSE;
167 snprintf(fname,1024, "%s",gSystem->ExpandPathName(contall));
168 //PH ifstream fincont(fname);
169 fstream fincont(fname,ios::in);
170 #if defined(__HP_aCC) || defined(__DECCXX)
171 if(!fincont.rdbuf()->is_open()) return -1;
173 if(!fincont.is_open()) return -1;
177 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
178 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
179 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>fcaq[13][nn]>>
180 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>fcaq[18][nn]>>
181 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>fcaq[23][nn]>>
182 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>fcaq[28][nn]>>
183 fcaq[29][nn]>>fcaq[30][nn]>>fcaq[31][nn]>>fcaq[32][nn]>>fcaq[33][nn])
190 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
191 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
192 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>fcag[13][nn]>>
193 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>fcag[18][nn]>>
194 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>fcag[23][nn]>>
195 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>fcag[28][nn]>>
196 fcag[29][nn]>>fcag[30][nn]>>fcag[31][nn]>>fcag[32][nn]>>fcag[33][nn])
203 snprintf(fname,1024, "%s",gSystem->ExpandPathName(discall));
204 //PH ifstream findisc(fname);
205 fstream findisc(fname,ios::in);
206 #if defined(__HP_aCC) || defined(__DECCXX)
207 if(!findisc.rdbuf()->is_open()) return -1;
209 if(!findisc.is_open()) return -1;
213 while(findisc>>frrr[nn]>>fdaq[nn]) {
218 while(findisc>>frrrg[nn]>>fdag[nn]) {
223 fTablesLoaded = kTRUE;
228 C***************************************************************************
229 C Quenching Weights for Multiple Soft Scattering
234 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
236 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
239 C This package contains quenching weights for gluon radiation in the
240 C multiple soft scattering approximation.
242 C swqmult returns the quenching weight for a quark (ipart=1) or
243 C a gluon (ipart=2) traversing a medium with transport coeficient q and
244 C length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where
245 C wc=0.5*q*L^2 and w is the energy radiated. The output values are
246 C the continuous and discrete (prefactor of the delta function) parts
247 C of the quenching weights.
249 C In order to use this routine, the files cont.all and disc.all need to be
250 C in the working directory.
252 C An initialization of the tables is needed by doing call initmult before
255 C Please, send us any comment:
257 C urs.wiedemann@cern.ch
258 C carlos.salgado@cern.ch
261 C-------------------------------------------------------------------
263 SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete)
265 REAL*8 xx(400), daq(34), caq(34,261), rrr(34)
266 COMMON /dataqua/ xx, daq, caq, rrr
268 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
269 COMMON /dataglu/ xxg, dag, cag, rrrg
271 REAL*8 rrrr,xxxx, continuous, discrete
273 INTEGER nrlow, nrhigh, nxlow, nxhigh
274 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
275 REAL*8 xfraclow, xfrachigh
286 if (rrin.lt.rrr(nr)) then
298 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
299 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
300 if (rrin.gt.10000d0) then
301 rfraclow = dlog(rrhigh/rrin)/dlog(rrhigh/rrlow)
302 rfrachigh = dlog(rrin/rrlow)/dlog(rrhigh/rrlow)
305 if (ipart.eq.1.and.rrin.ge.rrr(1)) then
312 if (ipart.ne.1.and.rrin.ge.rrrg(1)) then
319 if (xxxx.ge.xx(261)) go to 245
321 nxlow = int(xxin/0.01) + 1
323 xfraclow = (xx(nxhigh)-xxin)/0.01
324 xfrachigh = (xxin - xx(nxlow))/0.01
327 clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh)
328 chigh = xfraclow*caq(nrhigh,nxlow)+xfrachigh*caq(nrhigh,nxhigh)
330 clow = xfraclow*cag(nrlow,nxlow)+xfrachigh*cag(nrlow,nxhigh)
331 chigh = xfraclow*cag(nrhigh,nxlow)+xfrachigh*cag(nrhigh,nxhigh)
334 continuous = rfraclow*clow + rfrachigh*chigh
339 discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh)
341 discrete = rfraclow*dag(nrlow) + rfrachigh*dag(nrhigh)
347 REAL*8 xxq(400), daq(34), caq(34,261), rrr(34)
348 COMMON /dataqua/ xxq, daq, caq, rrr
350 REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34)
351 COMMON /dataglu/ xxg, dag, cag, rrrg
353 OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90)
355 read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn),
356 + caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn),
357 + caq(9,nn), caq(10,nn), caq(11,nn), caq(12,nn),
359 + caq(14,nn), caq(15,nn), caq(16,nn), caq(17,nn),
361 + caq(19,nn), caq(20,nn), caq(21,nn), caq(22,nn),
363 + caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn),
365 + caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn),
366 + caq(33,nn), caq(34,nn)
369 read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn),
370 + cag(4,nn), cag(5,nn), cag(6,nn), cag(7,nn), cag(8,nn),
371 + cag(9,nn), cag(10,nn), cag(11,nn), cag(12,nn),
373 + cag(14,nn), cag(15,nn), cag(16,nn), cag(17,nn),
375 + cag(19,nn), cag(20,nn), cag(21,nn), cag(22,nn),
377 + cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn),
379 + cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn),
380 + cag(33,nn), cag(34,nn)
384 OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91)
386 read (21,*) rrr(nn), daq(nn)
389 read (21,*) rrrg(nn), dag(nn)
394 90 PRINT*, 'input - output error'
395 91 PRINT*, 'input - output error #2'
401 =======================================================================
403 Adapted to ROOT macro by A. Dainese - 13/07/2003
404 Ported to class by C. Loizides - 12/02/2004
405 New version for extended R values added - 06/03/2004
408 Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
409 Double_t &continuous,Double_t &discrete) const
411 // Calculate Multiple Scattering approx.
412 // weights for given parton type,
413 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
419 //read-in data before first call
421 Error("CalcMult","Tables are not loaded.");
425 Error("CalcMult","Tables are not loaded for Multiple Scattering.");
429 Double_t rrin = rrrr;
430 Double_t xxin = xxxx;
432 if(xxin>fxx[260]) return -1;
433 Int_t nxlow = (Int_t)(xxin/0.01) + 1;
434 Int_t nxhigh = nxlow + 1;
435 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.01;
436 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01;
439 if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD
440 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
442 Int_t nrlow=0,nrhigh=0;
443 Double_t rrhigh=0,rrlow=0;
444 for(Int_t nr=1; nr<=34; nr++) {
445 if(rrin<frrr[nr-1]) {
448 rrhigh = frrr[nr-1-1];
458 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
459 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
462 rfraclow = TMath::Log2(rrhigh/rrin)/TMath::Log2(rrhigh/rrlow);
463 rfrachigh = TMath::Log2(rrin/rrlow)/TMath::Log2(rrhigh/rrlow);
465 if((ipart==1) && (rrin>=frrr[0]))
472 if((ipart==2) && (rrin>=frrrg[0]))
480 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
482 Double_t clow=0,chigh=0;
484 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
485 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
487 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
488 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
491 continuous = rfraclow*clow + rfrachigh*chigh;
492 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
493 //rfraclow,clow,rfrachigh,chigh,continuous);
496 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
498 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
504 Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *discall)
506 // read in tables for Single Hard Approx.
507 // path to continuum and to discrete part
509 fTablesLoaded = kFALSE;
513 snprintf(fname, 1024, "%s",gSystem->ExpandPathName(contall));
514 //PH ifstream fincont(fname);
515 fstream fincont(fname,ios::in);
516 #if defined(__HP_aCC) || defined(__DECCXX)
517 if(!fincont.rdbuf()->is_open()) return -1;
519 if(!fincont.is_open()) return -1;
523 while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
524 fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
525 fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>
527 fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>
529 fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>
531 fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>
540 while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
541 fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
542 fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>
544 fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>
546 fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>
548 fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>
556 snprintf(fname, 1024, "%s",gSystem->ExpandPathName(discall));
557 //PH ifstream findisc(fname);
558 fstream findisc(fname,ios::in);
559 #if defined(__HP_aCC) || defined(__DECCXX)
560 if(!findisc.rdbuf()->is_open()) return -1;
562 if(!findisc.is_open()) return -1;
566 while(findisc>>frrr[nn]>>fdaq[nn]) {
571 while(findisc>>frrrg[nn]>>fdag[nn]) {
577 fTablesLoaded = kTRUE;
582 C***************************************************************************
583 C Quenching Weights for Single Hard Scattering
588 C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184.
590 C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002.
593 C This package contains quenching weights for gluon radiation in the
594 C single hard scattering approximation.
596 C swqlin returns the quenching weight for a quark (ipart=1) or
597 C a gluon (ipart=2) traversing a medium with Debye screening mass mu and
598 C length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where
599 C wc=0.5*mu^2*L and w is the energy radiated. The output values are
600 C the continuous and discrete (prefactor of the delta function) parts
601 C of the quenching weights.
603 C In order to use this routine, the files contlin.all and disclin.all
604 C need to be in the working directory.
606 C An initialization of the tables is needed by doing call initlin before
609 C Please, send us any comment:
611 C urs.wiedemann@cern.ch
612 C carlos.salgado@cern.ch
615 C-------------------------------------------------------------------
618 SUBROUTINE swqlin(ipart,rrrr,xxxx,continuous,discrete)
620 REAL*8 xx(400), dalq(30), calq(30,261), rrr(30)
621 COMMON /datalinqua/ xx, dalq, calq, rrr
623 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
624 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
626 REAL*8 rrrr,xxxx, continuous, discrete
628 INTEGER nrlow, nrhigh, nxlow, nxhigh
629 REAL*8 rrhigh, rrlow, rfraclow, rfrachigh
630 REAL*8 xfraclow, xfrachigh
636 nxlow = int(xxin/0.038) + 1
638 xfraclow = (xx(nxhigh)-xxin)/0.038
639 xfrachigh = (xxin - xx(nxlow))/0.038
642 if (rrin.lt.rrr(nr)) then
654 rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
655 rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
658 clow = xfraclow*calq(nrlow,nxlow)+xfrachigh*calq(nrlow,nxhigh)
659 chigh = xfraclow*calq(nrhigh,nxlow)+xfrachigh*calq(nrhigh,nxhigh)
661 clow = xfraclow*calg(nrlow,nxlow)+xfrachigh*calg(nrlow,nxhigh)
662 chigh = xfraclow*calg(nrhigh,nxlow)+xfrachigh*calg(nrhigh,nxhigh)
665 continuous = rfraclow*clow + rfrachigh*chigh
668 discrete = rfraclow*dalq(nrlow) + rfrachigh*dalq(nrhigh)
670 discrete = rfraclow*dalg(nrlow) + rfrachigh*dalg(nrhigh)
676 REAL*8 xxlq(400), dalq(30), calq(30,261), rrr(30)
677 COMMON /datalinqua/ xxlq, dalq, calq, rrr
679 REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30)
680 COMMON /datalinglu/ xxlg, dalg, calg, rrrlg
682 OPEN(UNIT=20,FILE='contlin.all',STATUS='OLD',ERR=90)
684 read (20,*) xxlq(nn), calq(1,nn), calq(2,nn), calq(3,nn),
685 + calq(4,nn), calq(5,nn), calq(6,nn), calq(7,nn), calq(8,nn),
686 + calq(9,nn), calq(10,nn), calq(11,nn), calq(12,nn),
688 + calq(14,nn), calq(15,nn), calq(16,nn), calq(17,nn),
690 + calq(19,nn), calq(20,nn), calq(21,nn), calq(22,nn),
692 + calq(24,nn), calq(25,nn), calq(26,nn), calq(27,nn),
694 + calq(29,nn), calq(30,nn)
697 read (20,*) xxlg(nn), calg(1,nn), calg(2,nn), calg(3,nn),
698 + calg(4,nn), calg(5,nn), calg(6,nn), calg(7,nn), calg(8,nn),
699 + calg(9,nn), calg(10,nn), calg(11,nn), calg(12,nn),
701 + calg(14,nn), calg(15,nn), calg(16,nn), calg(17,nn),
703 + calg(19,nn), calg(20,nn), calg(21,nn), calg(22,nn),
705 + calg(24,nn), calg(25,nn), calg(26,nn), calg(27,nn),
707 + calg(29,nn), calg(30,nn)
711 OPEN(UNIT=21,FILE='disclin.all',STATUS='OLD',ERR=91)
713 read (21,*) rrr(nn), dalq(nn)
716 read (21,*) rrrlg(nn), dalg(nn)
721 90 PRINT*, 'input - output error'
722 91 PRINT*, 'input - output error #2'
727 =======================================================================
729 Ported to class by C. Loizides - 17/02/2004
733 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx,
734 Double_t &continuous,Double_t &discrete) const
736 // calculate Single Hard approx.
737 // weights for given parton type,
738 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
740 // read-in data before first call
742 Error("CalcSingleHard","Tables are not loaded.");
746 Error("CalcSingleHard","Tables are not loaded for Single Hard Scattering.");
750 Double_t rrin = rrrr;
751 Double_t xxin = xxxx;
753 Int_t nxlow = (Int_t)(xxin/0.038) + 1;
754 Int_t nxhigh = nxlow + 1;
755 Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.038;
756 Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.038;
759 if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD
760 if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD
762 Int_t nrlow=0,nrhigh=0;
763 Double_t rrhigh=0,rrlow=0;
764 for(Int_t nr=1; nr<=30; nr++) {
765 if(rrin<frrr[nr-1]) {
768 rrhigh = frrr[nr-1-1];
778 Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
779 Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
781 //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD
783 Double_t clow=0,chigh=0;
785 clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1];
786 chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1];
788 clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1];
789 chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1];
792 continuous = rfraclow*clow + rfrachigh*chigh;
793 //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
794 // rfraclow,clow,rfrachigh,chigh,continuous);
797 discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
799 discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1];
805 Int_t AliQuenchingWeights::CalcMult(Int_t ipart,
806 Double_t w,Double_t qtransp,Double_t length,
807 Double_t &continuous,Double_t &discrete) const
809 // Calculate Multiple Scattering approx.
810 // weights for given parton type,
811 // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
813 Double_t wc=CalcWC(qtransp,length);
814 Double_t rrrr=CalcR(wc,length);
816 return CalcMult(ipart,rrrr,xxxx,continuous,discrete);
819 Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
820 Double_t w,Double_t mu,Double_t length,
821 Double_t &continuous,Double_t &discrete) const
823 // calculate Single Hard approx.
824 // weights for given parton type,
825 // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
827 Double_t wcbar=CalcWCbar(mu,length);
828 Double_t rrrr=CalcR(wcbar,length);
829 Double_t xxxx=w/wcbar;
830 return CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
833 Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const
835 //calculate r value and
836 //check if it is less then maximum
838 Double_t r = wc*l*fgkConvFmToInvGeV;
840 Warning("CalcR","Value of r = %.2f; should be less than %.2f", r, fgkRMax);
846 Double_t AliQuenchingWeights::CalcRk(Double_t k, Double_t I0, Double_t I1) const
848 //calculate R value and
849 //check if it is less then maximum
851 Double_t r = fgkRMax-1;
855 Warning("CalcRk","Value of r = %.2f; should be less than %.2f",r,fgkRMax);
861 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const
863 // return DeltaE for MS or SH scattering
864 // for given parton type, length and energy
865 // Dependant on ECM (energy constraint method)
866 // e is used to determine where to set bins to zero.
869 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
872 if((ipart<1) || (ipart>2)) {
873 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
877 Int_t l=GetIndex(length);
880 if(fECMethod==kReweight){
884 ret=fHistos[ipart-1][l-1]->GetRandom();
886 Warning("GetELossRandom",
887 "Stopping reweighting; maximum loss assigned after 1e6 trials.");
894 Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
899 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const
901 //return quenched parton energy
902 //for given parton type, length and energy
904 Double_t loss=GetELossRandom(ipart,length,e);
908 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const
910 // return DeltaE for MS or SH scattering
911 // for given parton type, length distribution and energy
912 // Dependant on ECM (energy constraint method)
913 // e is used to determine where to set bins to zero.
916 Warning("GetELossRandom","Pointer to length distribution is NULL.");
919 Double_t ell=hell->GetRandom();
920 return GetELossRandom(ipart,ell,e);
923 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e) const
925 //return quenched parton energy
926 //for given parton type, length distribution and energy
928 Double_t loss=GetELossRandom(ipart,hell,e);
932 Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
934 // return DeltaE for new dynamic version
935 // for given parton type, I0 and I1 value and energy
936 // Dependant on ECM (energy constraint method)
937 // e is used to determine where to set bins to zero.
939 // read-in data before first call
941 Fatal("GetELossRandomK","Tables are not loaded.");
944 if((ipart<1) || (ipart>2)) {
945 Fatal("GetELossRandomK","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
949 Double_t r=CalcRk(I0,I1);
951 Fatal("GetELossRandomK","R should not be negative");
954 Double_t wc=CalcWCk(I1);
956 Fatal("GetELossRandomK","wc should be greater than zero");
959 if(SampleEnergyLoss(ipart,r)!=0){
960 Fatal("GetELossRandomK","Could not sample energy loss");
964 if(fECMethod==kReweight){
968 ret=fHisto->GetRandom();
970 Warning("GetELossRandomK",
971 "Stopping reweighting; maximum loss assigned after 1e6 trials.");
979 Double_t ret=fHisto->GetRandom()*wc;
984 Double_t AliQuenchingWeights::CalcQuenchedEnergyK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
986 //return quenched parton energy
987 //for given parton type, I0 and I1 value and energy
989 Double_t loss=GetELossRandomK(ipart,I0,I1,e);
993 Double_t AliQuenchingWeights::GetELossRandomKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
995 // return DeltaE for new dynamic version
996 // for given parton type, I0 and I1 value and energy
997 // Dependant on ECM (energy constraint method)
998 // e is used to determine where to set bins to zero.
999 // method is optimized and should only be used if
1000 // all parameters are well within the bounds.
1001 // read-in data tables before first call
1003 Double_t r=CalcRk(I0,I1);
1008 Double_t wc=CalcWCk(I1);
1013 return GetELossRandomKFastR(ipart,r,wc,e);
1016 Double_t AliQuenchingWeights::GetELossRandomKFastR(Int_t ipart, Double_t r, Double_t wc, Double_t e)
1018 // return DeltaE for new dynamic version
1019 // for given parton type, R and wc value and energy
1020 // Dependant on ECM (energy constraint method)
1021 // e is used to determine where to set bins to zero.
1022 // method is optimized and should only be used if
1023 // all parameters are well within the bounds.
1024 // read-in data tables before first call
1030 Double_t discrete=0.;
1031 Double_t continuous=0.;
1033 Double_t xxxx = fHisto->GetBinCenter(bin);
1035 CalcMult(ipart,r,xxxx,continuous,discrete);
1037 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1040 return 0.; //no energy loss
1042 if (!fHisto) Init();
1044 fHisto->SetBinContent(bin,continuous);
1045 Int_t kbinmax=fHisto->FindBin(e/wc);
1046 if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
1047 if(kbinmax==1) return e; //maximum energy loss
1050 for(bin=2; bin<=kbinmax; bin++) {
1051 xxxx = fHisto->GetBinCenter(bin);
1052 CalcMult(ipart,r,xxxx,continuous,discrete);
1053 fHisto->SetBinContent(bin,continuous);
1056 for(bin=2; bin<=kbinmax; bin++) {
1057 xxxx = fHisto->GetBinCenter(bin);
1058 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1059 fHisto->SetBinContent(bin,continuous);
1063 if(fECMethod==kReweight){
1064 fHisto->SetBinContent(kbinmax+1,0);
1065 fHisto->Fill(0.,discrete*fgkBins/fgkMaxBin);
1066 } else if (fECMethod==kReweightCont) {
1067 fHisto->SetBinContent(kbinmax+1,0);
1068 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1069 fHisto->Scale(1./kdelta*(1-discrete));
1070 fHisto->Fill(0.,discrete);
1072 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1073 Double_t val=discrete*fgkBins/fgkMaxBin;
1074 fHisto->Fill(0.,val);
1075 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1077 for(bin=kbinmax+2; bin<=fgkBins; bin++) {
1078 fHisto->SetBinContent(bin,0);
1080 //cout << kbinmax << " " << discrete << " " << fHisto->Integral() << endl;
1081 Double_t ret=fHisto->GetRandom()*wc;
1086 Double_t AliQuenchingWeights::CalcQuenchedEnergyKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
1088 //return quenched parton energy (fast method)
1089 //for given parton type, I0 and I1 value and energy
1091 Double_t loss=GetELossRandomKFast(ipart,I0,I1,e);
1095 Double_t AliQuenchingWeights::GetDiscreteWeight(Int_t ipart, Double_t I0, Double_t I1)
1097 // return discrete weight
1099 Double_t r=CalcRk(I0,I1);
1103 return GetDiscreteWeightR(ipart,r);
1106 Double_t AliQuenchingWeights::GetDiscreteWeightR(Int_t ipart, Double_t r)
1108 // return discrete weight
1114 Double_t discrete=0.;
1115 Double_t continuous=0.;
1117 Double_t xxxx = fHisto->GetBinCenter(bin);
1119 CalcMult(ipart,r,xxxx,continuous,discrete);
1121 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1125 void AliQuenchingWeights::GetZeroLossProb(Double_t &p,Double_t &prw,Double_t &prwcont,
1126 Int_t ipart,Double_t I0,Double_t I1,Double_t e)
1128 //calculate the probabilty that there is no energy
1129 //loss for different ways of energy constraint
1130 p=1.;prw=1.;prwcont=1.;
1131 Double_t r=CalcRk(I0,I1);
1135 Double_t wc=CalcWCk(I1);
1139 GetZeroLossProbR(p,prw,prwcont,ipart,r,wc,e);
1142 void AliQuenchingWeights::GetZeroLossProbR(Double_t &p,Double_t &prw,Double_t &prwcont,
1143 Int_t ipart, Double_t r,Double_t wc,Double_t e)
1145 //calculate the probabilty that there is no energy
1146 //loss for different ways of energy constraint
1151 Double_t discrete=0.;
1152 Double_t continuous=0.;
1153 if (!fHisto) Init();
1154 Int_t kbinmax=fHisto->FindBin(e/wc);
1155 if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
1157 for(Int_t bin=1; bin<=kbinmax; bin++) {
1158 Double_t xxxx = fHisto->GetBinCenter(bin);
1159 CalcMult(ipart,r,xxxx,continuous,discrete);
1160 fHisto->SetBinContent(bin,continuous);
1163 for(Int_t bin=1; bin<=kbinmax; bin++) {
1164 Double_t xxxx = fHisto->GetBinCenter(bin);
1165 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1166 fHisto->SetBinContent(bin,continuous);
1170 //non-reweighted P(Delta E = 0)
1171 const Double_t kdelta=fHisto->Integral(1,kbinmax);
1172 Double_t val=discrete*fgkBins/fgkMaxBin;
1173 fHisto->Fill(0.,val);
1174 fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
1175 Double_t hint=fHisto->Integral(1,kbinmax+1);
1176 p=fHisto->GetBinContent(1)/hint;
1179 hint=fHisto->Integral(1,kbinmax);
1180 prw=fHisto->GetBinContent(1)/hint;
1182 Double_t xxxx = fHisto->GetBinCenter(1);
1183 CalcMult(ipart,r,xxxx,continuous,discrete);
1184 fHisto->SetBinContent(1,continuous);
1185 hint=fHisto->Integral(1,kbinmax);
1186 fHisto->Scale(1./hint*(1-discrete));
1187 fHisto->Fill(0.,discrete);
1188 prwcont=fHisto->GetBinContent(1);
1191 Int_t AliQuenchingWeights::SampleEnergyLoss()
1193 // Has to be called to fill the histograms
1195 // For stored values fQTransport loop over
1196 // particle type and length = 1 to fMaxLength (fm)
1197 // to fill energy loss histos
1199 // Take histogram of continuous weights
1200 // Take discrete_weight
1201 // If discrete_weight > 1, put all channels to 0, except channel 1
1202 // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral
1204 // read-in data before first call
1206 Error("SampleEnergyLoss","Tables are not loaded.");
1211 Int_t lmax=CalcLengthMax(fQTransport);
1212 if(fLengthMax>lmax){
1213 Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
1217 Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested.");
1221 fHistos=new TH1F**[2];
1222 fHistos[0]=new TH1F*[4*fLengthMax];
1223 fHistos[1]=new TH1F*[4*fLengthMax];
1224 fLengthMaxOld=fLengthMax; //remember old value in case
1225 //user wants to reset
1228 Char_t meddesc[100];
1230 medvalue=(Int_t)(fQTransport*1000.);
1231 snprintf(meddesc, 100, "MS");
1233 medvalue=(Int_t)(fMu*1000.);
1234 snprintf(meddesc, 100, "SH");
1237 for(Int_t ipart=1;ipart<=2;ipart++){
1238 for(Int_t l=1;l<=4*fLengthMax;l++){
1240 snprintf(hname, 100, "hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
1242 Double_t wc = CalcWC(len);
1243 fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgkBins,0.,fgkMaxBin*wc);
1244 fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
1245 fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
1246 fHistos[ipart-1][l-1]->SetLineColor(4);
1248 Double_t rrrr = CalcR(wc,len);
1249 Double_t discrete=0.;
1250 // loop on histogram channels
1251 for(Int_t bin=1; bin<=fgkBins; bin++) {
1252 Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
1253 Double_t continuous;
1255 CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1257 CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1258 fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
1260 // add discrete part to distribution
1262 for(Int_t bin=2;bin<=fgkBins;bin++)
1263 fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
1265 Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1266 fHistos[ipart-1][l-1]->Fill(0.,val);
1268 Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins);
1269 fHistos[ipart-1][l-1]->Scale(1./hint);
1275 Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t r)
1277 // Sample energy loss directly for one particle type
1278 // choses R (safe it and keep it until next call of function)
1280 // read-in data before first call
1282 Error("SampleEnergyLoss","Tables are not loaded.");
1286 Double_t discrete=0.;
1287 Double_t continuous=0;;
1289 if (!fHisto) Init();
1290 Double_t xxxx = fHisto->GetBinCenter(bin);
1292 CalcMult(ipart,r,xxxx,continuous,discrete);
1294 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1297 fHisto->SetBinContent(1,1.);
1298 for(bin=2;bin<=fgkBins;bin++)
1299 fHisto->SetBinContent(bin,0.);
1303 fHisto->SetBinContent(bin,continuous);
1304 for(bin=2; bin<=fgkBins; bin++) {
1305 xxxx = fHisto->GetBinCenter(bin);
1307 CalcMult(ipart,r,xxxx,continuous,discrete);
1309 CalcSingleHard(ipart,r,xxxx,continuous,discrete);
1310 fHisto->SetBinContent(bin,continuous);
1313 Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
1314 fHisto->Fill(0.,val);
1315 Double_t hint=fHisto->Integral(1,fgkBins);
1317 fHisto->Scale(1./hint);
1319 //cout << discrete << " " << hint << " " << continuous << endl;
1325 const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) const
1327 //return quenching histograms
1328 //for ipart and length
1331 Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
1334 if((ipart<1) || (ipart>2)) {
1335 Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
1339 Int_t l=GetIndex(length);
1341 return fHistos[ipart-1][l-1];
1344 TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const
1346 // ipart = 1 for quark, 2 for gluon
1347 // medval a) qtransp = transport coefficient (GeV^2/fm)
1348 // b) mu = Debye mass (GeV)
1349 // length = path length in medium (fm)
1350 // Get from SW tables:
1351 // - continuous weight, as a function of dE/wc
1354 Char_t meddesc[100];
1356 wc=CalcWC(medval,length);
1357 snprintf(meddesc, 100, "MS");
1359 wc=CalcWCbar(medval,length);
1360 snprintf(meddesc, 100, "SH");
1364 snprintf(hname, 100, "hContQWHisto_%s_%d_%d_%d",meddesc,ipart,
1365 (Int_t)(medval*1000.),(Int_t)length);
1367 TH1F *hist = new TH1F("hist",hname,fgkBins,0.,fgkMaxBin*wc);
1368 hist->SetXTitle("#Delta E [GeV]");
1369 hist->SetYTitle("p(#Delta E)");
1370 hist->SetLineColor(4);
1372 Double_t rrrr = CalcR(wc,length);
1373 //loop on histogram channels
1374 for(Int_t bin=1; bin<=fgkBins; bin++) {
1375 Double_t xxxx = hist->GetBinCenter(bin)/wc;
1376 Double_t continuous,discrete;
1378 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1379 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1384 hist->SetBinContent(bin,continuous);
1389 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const
1391 // ipart = 1 for quark, 2 for gluon
1392 // medval a) qtransp = transport coefficient (GeV^2/fm)
1393 // b) mu = Debye mass (GeV)
1394 // length = path length in medium (fm)
1395 // Get from SW tables:
1396 // - continuous weight, as a function of dE/wc
1399 Char_t meddesc[100];
1401 wc=CalcWC(medval,length);
1402 snprintf(meddesc, 100, "MS");
1404 wc=CalcWCbar(medval,length);
1405 snprintf(meddesc, 100, "SH");
1409 snprintf(hname, 100, "hContQWHistox_%s_%d_%d_%d",meddesc,ipart,
1410 (Int_t)(medval*1000.),(Int_t)length);
1412 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
1413 histx->SetXTitle("x = #Delta E/#omega_{c}");
1415 histx->SetYTitle("p(#Delta E/#omega_{c})");
1417 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1418 histx->SetLineColor(4);
1420 Double_t rrrr = CalcR(wc,length);
1421 //loop on histogram channels
1422 for(Int_t bin=1; bin<=fgkBins; bin++) {
1423 Double_t xxxx = histx->GetBinCenter(bin);
1424 Double_t continuous,discrete;
1426 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1427 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1432 histx->SetBinContent(bin,continuous);
1437 TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t r) const
1439 // compute P(E) distribution for
1440 // given ipart = 1 for quark, 2 for gluon
1443 Char_t meddesc[100];
1445 snprintf(meddesc, 100, "MS");
1447 snprintf(meddesc, 100, "SH");
1451 snprintf(hname, 100, "hQWHistox_%s_%d_%.2f",meddesc,ipart,r);
1452 TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
1453 histx->SetXTitle("x = #Delta E/#omega_{c}");
1455 histx->SetYTitle("p(#Delta E/#omega_{c})");
1457 histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
1458 histx->SetLineColor(4);
1461 Double_t continuous=0.,discrete=0.;
1462 //loop on histogram channels
1463 for(Int_t bin=1; bin<=fgkBins; bin++) {
1464 Double_t xxxx = histx->GetBinCenter(bin);
1466 if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
1467 else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
1472 histx->SetBinContent(bin,continuous);
1475 //add discrete part to distribution
1477 for(Int_t bin=2;bin<=fgkBins;bin++)
1478 histx->SetBinContent(bin,0.);
1480 Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins);
1481 histx->Fill(0.,val);
1483 Double_t hint=histx->Integral(1,fgkBins);
1484 if(hint!=0) histx->Scale(1./hint);
1489 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
1491 // compute energy loss histogram for
1492 // parton type, medium value, length and energy
1494 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1496 dummy->SetQTransport(medval);
1499 dummy->SetMu(medval);
1500 dummy->InitSingleHard();
1502 dummy->SampleEnergyLoss();
1507 snprintf(name, 100, "Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1508 snprintf(hname,100, "hLossQuarks");
1510 snprintf(name, 100, "Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1511 snprintf(hname, 100, "hLossGluons");
1514 TH1F *h = new TH1F(hname,name,250,0,250);
1515 for(Int_t i=0;i<100000;i++){
1516 //if(i % 1000 == 0) cout << "." << flush;
1517 Double_t loss=dummy->GetELossRandom(ipart,l,e);
1525 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
1527 // compute energy loss histogram for
1528 // parton type, medium value,
1529 // length distribution and energy
1531 AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
1533 dummy->SetQTransport(medval);
1536 dummy->SetMu(medval);
1537 dummy->InitSingleHard();
1539 dummy->SampleEnergyLoss();
1544 snprintf(name, 100, "Energy Loss Distribution - Quarks;E_{loss} (GeV);#");
1545 snprintf(hname,100, "hLossQuarks");
1547 snprintf(name,100, "Energy Loss Distribution - Gluons;E_{loss} (GeV);#");
1548 snprintf(hname, 100, "hLossGluons");
1551 TH1F *h = new TH1F(hname,name,250,0,250);
1552 for(Int_t i=0;i<100000;i++){
1553 //if(i % 1000 == 0) cout << "." << flush;
1554 Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
1562 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t r) const
1564 // compute energy loss histogram for
1565 // parton type and given R
1567 TH1F *dummy = ComputeQWHistoX(ipart,r);
1568 if(!dummy) return 0;
1571 snprintf(hname, 100, "hELossHistox_%d_%.2f",ipart,r);
1572 TH1F *histx = new TH1F("histxr",hname,fgkBins,0.,fgkMaxBin);
1573 for(Int_t i=0;i<100000;i++){
1574 //if(i % 1000 == 0) cout << "." << flush;
1575 Double_t loss=dummy->GetRandom();
1582 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,Double_t l) const
1584 // compute average energy loss for
1585 // parton type, medium value, length and energy
1587 TH1F *dummy = ComputeELossHisto(ipart,medval,l);
1588 if(!dummy) return 0;
1589 Double_t ret=dummy->GetMean();
1594 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEll) const
1596 // compute average energy loss for
1597 // parton type, medium value,
1598 // length distribution and energy
1600 TH1F *dummy = ComputeELossHisto(ipart,medval,hEll);
1601 if(!dummy) return 0;
1602 Double_t ret=dummy->GetMean();
1607 Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t r) const
1609 // compute average energy loss over wc
1610 // for parton type and given R
1612 TH1F *dummy = ComputeELossHisto(ipart,r);
1613 if(!dummy) return 0;
1614 Double_t ret=dummy->GetMean();
1619 void AliQuenchingWeights::PlotDiscreteWeights(Double_t len,Double_t qm) const
1621 // plot discrete weights for given length
1625 c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
1627 c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
1630 TH2F *hframe = new TH2F("hdisc","",2,0,qm+.1,2,0,1.25);
1631 hframe->SetStats(0);
1633 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1635 hframe->SetXTitle("#mu [GeV]");
1636 //hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
1637 hframe->SetYTitle("p_{0} (discrete weight)");
1640 Int_t points=(Int_t)qm*4;
1641 TGraph *gq=new TGraph(points);
1644 for(Double_t q=0.05;q<=qm+.05;q+=0.25){
1646 CalcMult(1,1.0,q,len,cont,disc);
1647 gq->SetPoint(i,q,disc);i++;
1650 for(Double_t m=0.05;m<=qm+.05;m+=0.25){
1652 CalcSingleHard(1,1.0,m,len,cont, disc);
1653 gq->SetPoint(i,m,disc);i++;
1656 gq->SetMarkerStyle(20);
1657 gq->SetMarkerColor(1);
1658 gq->SetLineStyle(1);
1659 gq->SetLineColor(1);
1662 TGraph *gg=new TGraph(points);
1665 for(Double_t q=0.05;q<=qm+.05;q+=0.25){
1667 CalcMult(2,1.0,q,len,cont,disc);
1668 gg->SetPoint(i,q,disc);i++;
1671 for(Double_t m=0.05;m<=qm+.05;m+=0.25){
1673 CalcSingleHard(2,1.0,m,len,cont,disc);
1674 gg->SetPoint(i,m,disc);i++;
1677 gg->SetMarkerStyle(24);
1678 gg->SetMarkerColor(2);
1679 gg->SetLineStyle(2);
1680 gg->SetLineColor(2);
1683 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1684 l1a->SetFillStyle(0);
1685 l1a->SetBorderSize(0);
1687 snprintf(label, 100, "L = %.1f fm",len);
1688 l1a->AddEntry(gq,label,"");
1689 l1a->AddEntry(gq,"quark projectile","l");
1690 l1a->AddEntry(gg,"gluon projectile","l");
1696 void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
1698 // plot continous weights for
1699 // given parton type and length
1706 snprintf(title, 1024, "Cont. Weight for Multiple Scattering - Quarks");
1708 snprintf(title, 1024, "Cont. Weight for Multiple Scattering - Gluons");
1710 medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
1711 snprintf(name, 1024, "ccont-ms-%d",itype);
1714 snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Quarks");
1716 snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Gluons");
1718 medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
1719 snprintf(name, 1024, "ccont-ms-%d",itype);
1722 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1724 TH1F *h1=ComputeQWHisto(itype,medvals[0],ell);
1726 h1->SetTitle(title);
1728 h1->SetLineColor(1);
1730 TH1F *h2=ComputeQWHisto(itype,medvals[1],ell);
1732 h2->SetLineColor(2);
1733 h2->DrawCopy("SAME");
1734 TH1F *h3=ComputeQWHisto(itype,medvals[2],ell);
1736 h3->SetLineColor(3);
1737 h3->DrawCopy("SAME");
1739 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1740 l1a->SetFillStyle(0);
1741 l1a->SetBorderSize(0);
1743 snprintf(label, 100, "L = %.1f fm",ell);
1744 l1a->AddEntry(h1,label,"");
1746 snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
1747 l1a->AddEntry(h1,label,"pl");
1748 snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
1749 l1a->AddEntry(h2,label,"pl");
1750 snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
1751 l1a->AddEntry(h3, label,"pl");
1753 snprintf(label, 100, "#mu = %.1f GeV",medvals[0]);
1754 l1a->AddEntry(h1,label,"pl");
1755 snprintf(label, 100, "#mu = %.1f GeV",medvals[1]);
1756 l1a->AddEntry(h2,label,"pl");
1757 snprintf(label, 100, "#mu = %.1f GeV",medvals[2]);
1758 l1a->AddEntry(h3,label,"pl");
1765 void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
1767 // plot continous weights for
1768 // given parton type and medium value
1774 snprintf(title,1024, "Cont. Weight for Multiple Scattering - Quarks");
1776 snprintf(title,1024, "Cont. Weight for Multiple Scattering - Gluons");
1778 snprintf(name,1024, "ccont2-ms-%d",itype);
1781 snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Quarks");
1783 snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Gluons");
1785 snprintf(name, 1024, "ccont2-sh-%d",itype);
1787 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1789 TH1F *h1=ComputeQWHisto(itype,medval,8);
1791 h1->SetTitle(title);
1793 h1->SetLineColor(1);
1795 TH1F *h2=ComputeQWHisto(itype,medval,5);
1797 h2->SetLineColor(2);
1798 h2->DrawCopy("SAME");
1799 TH1F *h3=ComputeQWHisto(itype,medval,2);
1801 h3->SetLineColor(3);
1802 h3->DrawCopy("SAME");
1804 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1805 l1a->SetFillStyle(0);
1806 l1a->SetBorderSize(0);
1809 snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medval);
1811 snprintf(label, 100, "#mu = %.1f GeV",medval);
1813 l1a->AddEntry(h1,label,"");
1814 l1a->AddEntry(h1,"L = 8 fm","pl");
1815 l1a->AddEntry(h2,"L = 5 fm","pl");
1816 l1a->AddEntry(h3,"L = 2 fm","pl");
1822 void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t qm,Double_t e) const
1824 // plot average energy loss for given length
1825 // and parton energy
1828 Error("PlotAvgELoss","Tables are not loaded.");
1835 snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
1836 snprintf(name, 1024, "cavgelossms");
1838 snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering");
1839 snprintf(name, 1024, "cavgelosssh");
1842 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1844 TH2F *hframe = new TH2F("avgloss","",2,0,qm+.1,2,0,100);
1845 hframe->SetStats(0);
1847 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1849 hframe->SetXTitle("#mu [GeV]");
1850 hframe->SetYTitle("<E_{loss}> [GeV]");
1853 TGraph *gq=new TGraph(20);
1855 for(Double_t v=0.05;v<=qm+.05;v+=0.25){
1856 TH1F *dummy=ComputeELossHisto(1,v,len,e);
1857 Double_t avgloss=dummy->GetMean();
1858 gq->SetPoint(i,v,avgloss);i++;
1861 gq->SetMarkerStyle(21);
1864 Int_t points=(Int_t)qm*4;
1865 TGraph *gg=new TGraph(points);
1867 for(Double_t v=0.05;v<=qm+.05;v+=0.25){
1868 TH1F *dummy=ComputeELossHisto(2,v,len,e);
1869 Double_t avgloss=dummy->GetMean();
1870 gg->SetPoint(i,v,avgloss);i++;
1873 gg->SetMarkerStyle(20);
1874 gg->SetMarkerColor(2);
1877 TGraph *gratio=new TGraph(points);
1878 for(i=0;i<points;i++){
1880 gg->GetPoint(i,x,y);
1881 gq->GetPoint(i,x2,y2);
1883 gratio->SetPoint(i,x,y/y2*10/2.25);
1884 else gratio->SetPoint(i,x,0);
1886 gratio->SetLineStyle(4);
1888 TLegend *l1a = new TLegend(0.15,0.60,0.50,0.90);
1889 l1a->SetFillStyle(0);
1890 l1a->SetBorderSize(0);
1892 snprintf(label, 100, "L = %.1f fm",len);
1893 l1a->AddEntry(gq,label,"");
1894 l1a->AddEntry(gq,"quark projectile","pl");
1895 l1a->AddEntry(gg,"gluon projectile","pl");
1896 l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1902 void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
1904 // plot average energy loss for given
1905 // length distribution and parton energy
1908 Error("PlotAvgELossVs","Tables are not loaded.");
1915 snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
1916 snprintf(name, 1024, "cavgelossms2");
1918 snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering");
1919 snprintf(name, 1024, "cavgelosssh2");
1922 TCanvas *c = new TCanvas(name,title,0,0,500,400);
1924 TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100);
1925 hframe->SetStats(0);
1927 hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
1929 hframe->SetXTitle("#mu [GeV]");
1930 hframe->SetYTitle("<E_{loss}> [GeV]");
1933 TGraph *gq=new TGraph(20);
1935 for(Double_t v=0.05;v<=5.05;v+=0.25){
1936 TH1F *dummy=ComputeELossHisto(1,v,hEll,e);
1937 Double_t avgloss=dummy->GetMean();
1938 gq->SetPoint(i,v,avgloss);i++;
1941 gq->SetMarkerStyle(20);
1944 TGraph *gg=new TGraph(20);
1946 for(Double_t v=0.05;v<=5.05;v+=0.25){
1947 TH1F *dummy=ComputeELossHisto(2,v,hEll,e);
1948 Double_t avgloss=dummy->GetMean();
1949 gg->SetPoint(i,v,avgloss);i++;
1952 gg->SetMarkerStyle(24);
1955 TGraph *gratio=new TGraph(20);
1958 gg->GetPoint(i,x,y);
1959 gq->GetPoint(i,x2,y2);
1961 gratio->SetPoint(i,x,y/y2*10/2.25);
1962 else gratio->SetPoint(i,x,0);
1964 gratio->SetLineStyle(4);
1967 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
1968 l1a->SetFillStyle(0);
1969 l1a->SetBorderSize(0);
1971 snprintf(label, 100, "<L> = %.2f fm",hEll->GetMean());
1972 l1a->AddEntry(gq,label,"");
1973 l1a->AddEntry(gq,"quark","pl");
1974 l1a->AddEntry(gg,"gluon","pl");
1975 //l1a->AddEntry(gratio,"gluon/quark/2.25*10","pl");
1981 void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const
1983 // plot average energy loss versus ell
1986 Error("PlotAvgELossVsEll","Tables are not loaded.");
1994 snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
1995 snprintf(name, 1024, "cavgelosslms");
1998 snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering");
1999 snprintf(name, 1024, "cavgelosslsh");
2003 TCanvas *c = new TCanvas(name,title,0,0,600,400);
2005 TH2F *hframe = new TH2F("avglossell",title,2,0,fLengthMax,2,0,250);
2006 hframe->SetStats(0);
2007 hframe->SetXTitle("length [fm]");
2008 hframe->SetYTitle("<E_{loss}> [GeV]");
2011 TGraph *gq=new TGraph((Int_t)fLengthMax*4);
2013 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
2014 TH1F *dummy=ComputeELossHisto(1,medval,v,e);
2015 Double_t avgloss=dummy->GetMean();
2016 gq->SetPoint(i,v,avgloss);i++;
2019 gq->SetMarkerStyle(20);
2022 TGraph *gg=new TGraph((Int_t)fLengthMax*4);
2024 for(Double_t v=0.25;v<=fLengthMax;v+=0.25){
2025 TH1F *dummy=ComputeELossHisto(2,medval,v,e);
2026 Double_t avgloss=dummy->GetMean();
2027 gg->SetPoint(i,v,avgloss);i++;
2030 gg->SetMarkerStyle(24);
2033 TGraph *gratio=new TGraph((Int_t)fLengthMax*4);
2034 for(i=0;i<=(Int_t)fLengthMax*4;i++){
2036 gg->GetPoint(i,x,y);
2037 gq->GetPoint(i,x2,y2);
2039 gratio->SetPoint(i,x,y/y2*100/2.25);
2040 else gratio->SetPoint(i,x,0);
2042 gratio->SetLineStyle(1);
2043 gratio->SetLineWidth(2);
2045 TLegend *l1a = new TLegend(0.15,0.65,.60,0.85);
2046 l1a->SetFillStyle(0);
2047 l1a->SetBorderSize(0);
2050 snprintf(label, 100, "#hat{q} = %.2f GeV^{2}/fm",medval);
2052 snprintf(label, 100, "#mu = %.2f GeV",medval);
2053 l1a->AddEntry(gq,label,"");
2054 l1a->AddEntry(gq,"quark","pl");
2055 l1a->AddEntry(gg,"gluon","pl");
2056 l1a->AddEntry(gratio,"gluon/quark/2.25*100","pl");
2059 TF1 *f=new TF1("f","100",0,fLengthMax);
2066 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
2068 // plot relative energy loss for given
2069 // length and parton energy versus pt
2072 Error("PlotAvgELossVsPt","Tables are not loaded.");
2079 snprintf(title, 1024, "Relative Energy Loss for Multiple Scattering");
2080 snprintf(name, 1024, "cavgelossvsptms");
2082 snprintf(title, 1024, "Relative Energy Loss for Single Hard Scattering");
2083 snprintf(name, 1024, "cavgelossvsptsh");
2086 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2088 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2089 hframe->SetStats(0);
2090 hframe->SetXTitle("p_{T} [GeV]");
2091 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2094 TGraph *gq=new TGraph(40);
2096 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2097 TH1F *dummy=ComputeELossHisto(1,medval,len,pt);
2098 Double_t avgloss=dummy->GetMean();
2099 gq->SetPoint(i,pt,avgloss/pt);i++;
2102 gq->SetMarkerStyle(20);
2105 TGraph *gg=new TGraph(40);
2107 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2108 TH1F *dummy=ComputeELossHisto(2,medval,len,pt);
2109 Double_t avgloss=dummy->GetMean();
2110 gg->SetPoint(i,pt,avgloss/pt);i++;
2113 gg->SetMarkerStyle(24);
2116 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2117 l1a->SetFillStyle(0);
2118 l1a->SetBorderSize(0);
2120 snprintf(label, 100, "L = %.1f fm",len);
2121 l1a->AddEntry(gq,label,"");
2122 l1a->AddEntry(gq,"quark","pl");
2123 l1a->AddEntry(gg,"gluon","pl");
2129 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
2131 // plot relative energy loss for given
2132 // length distribution and parton energy versus pt
2135 Error("PlotAvgELossVsPt","Tables are not loaded.");
2142 snprintf(title, 1024, "Relative Energy Loss for Multiple Scattering");
2143 snprintf(name, 1024, "cavgelossvsptms2");
2145 snprintf(title, 1024, "Relative Energy Loss for Single Hard Scattering");
2146 snprintf(name, 1024, "cavgelossvsptsh2");
2148 TCanvas *c = new TCanvas(name,title,0,0,500,400);
2150 TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1);
2151 hframe->SetStats(0);
2152 hframe->SetXTitle("p_{T} [GeV]");
2153 hframe->SetYTitle("<E_{loss}>/p_{T} [GeV]");
2156 TGraph *gq=new TGraph(40);
2158 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2159 TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt);
2160 Double_t avgloss=dummy->GetMean();
2161 gq->SetPoint(i,pt,avgloss/pt);i++;
2164 gq->SetMarkerStyle(20);
2167 TGraph *gg=new TGraph(40);
2169 for(Double_t pt=2.5;pt<=100.05;pt+=2.5){
2170 TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt);
2171 Double_t avgloss=dummy->GetMean();
2172 gg->SetPoint(i,pt,avgloss/pt);i++;
2175 gg->SetMarkerStyle(24);
2178 TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
2179 l1a->SetFillStyle(0);
2180 l1a->SetBorderSize(0);
2182 snprintf(label, 100, "<L> = %.2f fm",hEll->GetMean());
2183 l1a->AddEntry(gq,label,"");
2184 l1a->AddEntry(gq,"quark","pl");
2185 l1a->AddEntry(gg,"gluon","pl");
2191 Int_t AliQuenchingWeights::GetIndex(Double_t len) const
2193 //get the index according to length
2194 if(len>fLengthMax) len=fLengthMax;
2196 Int_t l=Int_t(len/0.25);
2197 if((len-l*0.25)>0.125) l++;