/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* $Id$ */ //---------------------------------------------------------------------------- // Implementation of the class to calculate the parton energy loss // Based on the "BDMPS" quenching weights by C.A.Salgado and U.A.Wiedemann // // References: // C.A.Salgado and U.A.Wiedemann, Phys.Rev.D68 (2003) 014008 [hep-ph/0302184] // A.Dainese, Eur.Phys.J.C, in press, [nucl-ex/0312005] // // // Origin: C. Loizides constantinos.loizides@cern.ch // A. Dainese andrea.dainese@pd.infn.it // //---------------------------------------------------------------------------- #include #include #include #include #include #include #include #include #include "AliQuenchingWeights.h" ClassImp(AliQuenchingWeights) // conversion from fm to GeV^-1: 1 fm = fmGeV GeV^-1 const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197; // maximum value of R const Double_t AliQuenchingWeights::fgkRMax = 1.e6; // counter for histogram labels Int_t AliQuenchingWeights::fgCounter = 0; AliQuenchingWeights::AliQuenchingWeights() : TObject() { //default constructor fTablesLoaded=kFALSE; fMultSoft=kTRUE; fHistos=0; SetMu(); SetQTransport(); SetECMethod(); SetLengthMax(); fLengthMaxOld=0; fInstanceNumber=fgCounter++; } AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a) : TObject() { //copy constructor fTablesLoaded=kFALSE; fHistos=0; fLengthMaxOld=0; fMultSoft=a.GetMultSoft();; fMu=a.GetMu(); fQTransport=a.GetQTransport(); fECMethod=(kECMethod)a.GetECMethod(); fLengthMax=a.GetLengthMax(); fInstanceNumber=fgCounter++; //Missing in the class is the pathname //to the tables, can be added if needed } AliQuenchingWeights::~AliQuenchingWeights() { Reset(); } void AliQuenchingWeights::Reset() { //reset tables if there were used if(!fHistos) return; for(Int_t l=0;l parton energy) then sampled energy loss = parton energy."); else Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy."); } Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall) { // read in tables for multiple scattering approximation // path to continuum and to discrete part fTablesLoaded = kFALSE; fMultSoft=kTRUE; Char_t fname[1024]; sprintf(fname,"%s",gSystem->ExpandPathName(contall)); //PH ifstream fincont(fname); fstream fincont(fname,ios::in); #if defined(__HP_aCC) || defined(__DECCXX) if(!fincont.rdbuf()->is_open()) return -1; #else if(!fincont.is_open()) return -1; #endif Int_t nn=0; //quarks while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>> fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>> fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>fcaq[13][nn]>> fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>fcaq[18][nn]>> fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>fcaq[23][nn]>> fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>fcaq[28][nn]>> fcaq[29][nn]>>fcaq[30][nn]>>fcaq[31][nn]>>fcaq[32][nn]>>fcaq[33][nn]) { nn++; if(nn==261) break; } nn=0; //gluons while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>> fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>> fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>fcag[13][nn]>> fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>fcag[18][nn]>> fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>fcag[23][nn]>> fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>fcag[28][nn]>> fcag[29][nn]>>fcag[30][nn]>>fcag[31][nn]>>fcag[32][nn]>>fcag[33][nn]) { nn++; if(nn==261) break; } fincont.close(); sprintf(fname,"%s",gSystem->ExpandPathName(discall)); //PH ifstream findisc(fname); fstream findisc(fname,ios::in); #if defined(__HP_aCC) || defined(__DECCXX) if(!findisc.rdbuf()->is_open()) return -1; #else if(!findisc.is_open()) return -1; #endif nn=0; //quarks while(findisc>>frrr[nn]>>fdaq[nn]) { nn++; if(nn==34) break; } nn=0; //gluons while(findisc>>frrrg[nn]>>fdag[nn]) { nn++; if(nn==34) break; } findisc.close(); fTablesLoaded = kTRUE; return 0; } /* C*************************************************************************** C Quenching Weights for Multiple Soft Scattering C February 10, 2003 C C Refs: C C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184. C C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002. C C C This package contains quenching weights for gluon radiation in the C multiple soft scattering approximation. C C swqmult returns the quenching weight for a quark (ipart=1) or C a gluon (ipart=2) traversing a medium with transport coeficient q and C length L. The input values are rrrr=0.5*q*L^3 and xxxx=w/wc, where C wc=0.5*q*L^2 and w is the energy radiated. The output values are C the continuous and discrete (prefactor of the delta function) parts C of the quenching weights. C C In order to use this routine, the files cont.all and disc.all need to be C in the working directory. C C An initialization of the tables is needed by doing call initmult before C using swqmult. C C Please, send us any comment: C C urs.wiedemann@cern.ch C carlos.salgado@cern.ch C C C------------------------------------------------------------------- SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete) * REAL*8 xx(400), daq(34), caq(34,261), rrr(34) COMMON /dataqua/ xx, daq, caq, rrr * REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34) COMMON /dataglu/ xxg, dag, cag, rrrg REAL*8 rrrr,xxxx, continuous, discrete REAL*8 rrin, xxin INTEGER nrlow, nrhigh, nxlow, nxhigh REAL*8 rrhigh, rrlow, rfraclow, rfrachigh REAL*8 xfraclow, xfrachigh REAL*8 clow, chigh * continuous=0.d0 discrete=0.d0 rrin = rrrr xxin = xxxx * do 666, nr=1,34 if (rrin.lt.rrr(nr)) then rrhigh = rrr(nr) else rrhigh = rrr(nr-1) rrlow = rrr(nr) nrlow = nr nrhigh = nr-1 goto 665 endif 666 enddo 665 continue * rfraclow = (rrhigh-rrin)/(rrhigh-rrlow) rfrachigh = (rrin-rrlow)/(rrhigh-rrlow) if (rrin.gt.10000d0) then rfraclow = dlog(rrhigh/rrin)/dlog(rrhigh/rrlow) rfrachigh = dlog(rrin/rrlow)/dlog(rrhigh/rrlow) endif * if (ipart.eq.1.and.rrin.ge.rrr(1)) then nrlow=1 nrhigh=1 rfraclow=1 rfrachigh=0 endif if (ipart.ne.1.and.rrin.ge.rrrg(1)) then nrlow=1 nrhigh=1 rfraclow=1 rfrachigh=0 endif if (xxxx.ge.xx(261)) go to 245 nxlow = int(xxin/0.01) + 1 nxhigh = nxlow + 1 xfraclow = (xx(nxhigh)-xxin)/0.01 xfrachigh = (xxin - xx(nxlow))/0.01 * if(ipart.eq.1) then clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh) chigh = xfraclow*caq(nrhigh,nxlow)+xfrachigh*caq(nrhigh,nxhigh) else clow = xfraclow*cag(nrlow,nxlow)+xfrachigh*cag(nrlow,nxhigh) chigh = xfraclow*cag(nrhigh,nxlow)+xfrachigh*cag(nrhigh,nxhigh) endif continuous = rfraclow*clow + rfrachigh*chigh 245 continue if(ipart.eq.1) then discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh) else discrete = rfraclow*dag(nrlow) + rfrachigh*dag(nrhigh) endif * END subroutine initmult REAL*8 xxq(400), daq(34), caq(34,261), rrr(34) COMMON /dataqua/ xxq, daq, caq, rrr * REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34) COMMON /dataglu/ xxg, dag, cag, rrrg * OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90) do 110 nn=1,261 read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn), + caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn), + caq(9,nn), caq(10,nn), caq(11,nn), caq(12,nn), + caq(13,nn), + caq(14,nn), caq(15,nn), caq(16,nn), caq(17,nn), + caq(18,nn), + caq(19,nn), caq(20,nn), caq(21,nn), caq(22,nn), + caq(23,nn), + caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn), + caq(28,nn), + caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn), + caq(33,nn), caq(34,nn) 110 continue do 111 nn=1,261 read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn), + cag(4,nn), cag(5,nn), cag(6,nn), cag(7,nn), cag(8,nn), + cag(9,nn), cag(10,nn), cag(11,nn), cag(12,nn), + cag(13,nn), + cag(14,nn), cag(15,nn), cag(16,nn), cag(17,nn), + cag(18,nn), + cag(19,nn), cag(20,nn), cag(21,nn), cag(22,nn), + cag(23,nn), + cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn), + cag(28,nn), + cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn), + cag(33,nn), cag(34,nn) 111 continue close(20) * OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91) do 112 nn=1,34 read (21,*) rrr(nn), daq(nn) 112 continue do 113 nn=1,34 read (21,*) rrrg(nn), dag(nn) 113 continue close(21) * goto 888 90 PRINT*, 'input - output error' 91 PRINT*, 'input - output error #2' 888 continue end ======================================================================= Adapted to ROOT macro by A. Dainese - 13/07/2003 Ported to class by C. Loizides - 12/02/2004 New version for extended R values added - 06/03/2004 */ Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx, Double_t &continuous,Double_t &discrete) const { // Calculate Multiple Scattering approx. // weights for given parton type, // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2 //set result to zero continuous=0.; discrete=0.; // read-in data before first call if(!fTablesLoaded){ Error("CalcMult","Tables are not loaded."); return -1; } if(!fMultSoft){ Error("CalcMult","Tables are not loaded for Multiple Scattering."); return -1; } Double_t rrin = rrrr; Double_t xxin = xxxx; if(xxin>fxx[260]) return -1; Int_t nxlow = (Int_t)(xxin/0.01) + 1; Int_t nxhigh = nxlow + 1; Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.01; Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01; //why this? if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD Int_t nrlow=0,nrhigh=0; Double_t rrhigh=0,rrlow=0; for(Int_t nr=1; nr<=34; nr++) { if(rrin1.e4){ rfraclow = TMath::Log2(rrhigh/rrin)/TMath::Log2(rrhigh/rrlow); rfrachigh = TMath::Log2(rrin/rrlow)/TMath::Log2(rrhigh/rrlow); } if((ipart==1) && (rrin>=frrr[0])) { nrlow=1; nrhigh=1; rfraclow=1.; rfrachigh=0.; } if((ipart==2) && (rrin>=frrrg[0])) { nrlow=1; nrhigh=1; rfraclow=1.; rfrachigh=0.; } //printf("R = %f,\nRlow = %f, Rhigh = %f,\nRfraclow = %f, Rfrachigh = %f\n",rrin,rrlow,rrhigh,rfraclow,rfrachigh); // AD Double_t clow=0,chigh=0; if(ipart==1) { clow = xfraclow*fcaq[nrlow-1][nxlow-1]+xfrachigh*fcaq[nrlow-1][nxhigh-1]; chigh = xfraclow*fcaq[nrhigh-1][nxlow-1]+xfrachigh*fcaq[nrhigh-1][nxhigh-1]; } else { clow = xfraclow*fcag[nrlow-1][nxlow-1]+xfrachigh*fcag[nrlow-1][nxhigh-1]; chigh = xfraclow*fcag[nrhigh-1][nxlow-1]+xfrachigh*fcag[nrhigh-1][nxhigh-1]; } continuous = rfraclow*clow + rfrachigh*chigh; //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n", // rfraclow,clow,rfrachigh,chigh,continuous); if(ipart==1) { discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1]; } else { discrete = rfraclow*fdag[nrlow-1] + rfrachigh*fdag[nrhigh-1]; } return 0; } Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *discall) { // read in tables for Single Hard Approx. // path to continuum and to discrete part fTablesLoaded = kFALSE; fMultSoft=kFALSE; Char_t fname[1024]; sprintf(fname,"%s",gSystem->ExpandPathName(contall)); //PH ifstream fincont(fname); fstream fincont(fname,ios::in); #if defined(__HP_aCC) || defined(__DECCXX) if(!fincont.rdbuf()->is_open()) return -1; #else if(!fincont.is_open()) return -1; #endif Int_t nn=0; //quarks while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>> fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>> fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>> fcaq[13][nn]>> fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>> fcaq[18][nn]>> fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>> fcaq[23][nn]>> fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>> fcaq[28][nn]>> fcaq[29][nn]) { nn++; if(nn==261) break; } nn=0; //gluons while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>> fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>> fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>> fcag[13][nn]>> fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>> fcag[18][nn]>> fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>> fcag[23][nn]>> fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>> fcag[28][nn]>> fcag[29][nn]) { nn++; if(nn==261) break; } fincont.close(); sprintf(fname,"%s",gSystem->ExpandPathName(discall)); //PH ifstream findisc(fname); fstream findisc(fname,ios::in); #if defined(__HP_aCC) || defined(__DECCXX) if(!findisc.rdbuf()->is_open()) return -1; #else if(!findisc.is_open()) return -1; #endif nn=0; //quarks while(findisc>>frrr[nn]>>fdaq[nn]) { nn++; if(nn==30) break; } nn=0; //gluons while(findisc>>frrrg[nn]>>fdag[nn]) { nn++; if(nn==30) break; } findisc.close(); fTablesLoaded = kTRUE; return 0; } /* C*************************************************************************** C Quenching Weights for Single Hard Scattering C February 20, 2003 C C Refs: C C Carlos A. Salgado and Urs A. Wiedemann, hep-ph/0302184. C C Carlos A. Salgado and Urs A. Wiedemann Phys.Rev.Lett.89:092303,2002. C C C This package contains quenching weights for gluon radiation in the C single hard scattering approximation. C C swqlin returns the quenching weight for a quark (ipart=1) or C a gluon (ipart=2) traversing a medium with Debye screening mass mu and C length L. The input values are rrrr=0.5*mu^2*L^2 and xxxx=w/wc, where C wc=0.5*mu^2*L and w is the energy radiated. The output values are C the continuous and discrete (prefactor of the delta function) parts C of the quenching weights. C C In order to use this routine, the files contlin.all and disclin.all C need to be in the working directory. C C An initialization of the tables is needed by doing call initlin before C using swqlin. C C Please, send us any comment: C C urs.wiedemann@cern.ch C carlos.salgado@cern.ch C C C------------------------------------------------------------------- SUBROUTINE swqlin(ipart,rrrr,xxxx,continuous,discrete) * REAL*8 xx(400), dalq(30), calq(30,261), rrr(30) COMMON /datalinqua/ xx, dalq, calq, rrr * REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30) COMMON /datalinglu/ xxlg, dalg, calg, rrrlg REAL*8 rrrr,xxxx, continuous, discrete REAL*8 rrin, xxin INTEGER nrlow, nrhigh, nxlow, nxhigh REAL*8 rrhigh, rrlow, rfraclow, rfrachigh REAL*8 xfraclow, xfrachigh REAL*8 clow, chigh * rrin = rrrr xxin = xxxx * nxlow = int(xxin/0.038) + 1 nxhigh = nxlow + 1 xfraclow = (xx(nxhigh)-xxin)/0.038 xfrachigh = (xxin - xx(nxlow))/0.038 * do 666, nr=1,30 if (rrin.lt.rrr(nr)) then rrhigh = rrr(nr) else rrhigh = rrr(nr-1) rrlow = rrr(nr) nrlow = nr nrhigh = nr-1 goto 665 endif 666 enddo 665 continue * rfraclow = (rrhigh-rrin)/(rrhigh-rrlow) rfrachigh = (rrin-rrlow)/(rrhigh-rrlow) * if(ipart.eq.1) then clow = xfraclow*calq(nrlow,nxlow)+xfrachigh*calq(nrlow,nxhigh) chigh = xfraclow*calq(nrhigh,nxlow)+xfrachigh*calq(nrhigh,nxhigh) else clow = xfraclow*calg(nrlow,nxlow)+xfrachigh*calg(nrlow,nxhigh) chigh = xfraclow*calg(nrhigh,nxlow)+xfrachigh*calg(nrhigh,nxhigh) endif continuous = rfraclow*clow + rfrachigh*chigh if(ipart.eq.1) then discrete = rfraclow*dalq(nrlow) + rfrachigh*dalq(nrhigh) else discrete = rfraclow*dalg(nrlow) + rfrachigh*dalg(nrhigh) endif * END subroutine initlin REAL*8 xxlq(400), dalq(30), calq(30,261), rrr(30) COMMON /datalinqua/ xxlq, dalq, calq, rrr * REAL*8 xxlg(400), dalg(30), calg(30,261), rrrlg(30) COMMON /datalinglu/ xxlg, dalg, calg, rrrlg * OPEN(UNIT=20,FILE='contlin.all',STATUS='OLD',ERR=90) do 110 nn=1,261 read (20,*) xxlq(nn), calq(1,nn), calq(2,nn), calq(3,nn), + calq(4,nn), calq(5,nn), calq(6,nn), calq(7,nn), calq(8,nn), + calq(9,nn), calq(10,nn), calq(11,nn), calq(12,nn), + calq(13,nn), + calq(14,nn), calq(15,nn), calq(16,nn), calq(17,nn), + calq(18,nn), + calq(19,nn), calq(20,nn), calq(21,nn), calq(22,nn), + calq(23,nn), + calq(24,nn), calq(25,nn), calq(26,nn), calq(27,nn), + calq(28,nn), + calq(29,nn), calq(30,nn) 110 continue do 111 nn=1,261 read (20,*) xxlg(nn), calg(1,nn), calg(2,nn), calg(3,nn), + calg(4,nn), calg(5,nn), calg(6,nn), calg(7,nn), calg(8,nn), + calg(9,nn), calg(10,nn), calg(11,nn), calg(12,nn), + calg(13,nn), + calg(14,nn), calg(15,nn), calg(16,nn), calg(17,nn), + calg(18,nn), + calg(19,nn), calg(20,nn), calg(21,nn), calg(22,nn), + calg(23,nn), + calg(24,nn), calg(25,nn), calg(26,nn), calg(27,nn), + calg(28,nn), + calg(29,nn), calg(30,nn) 111 continue close(20) * OPEN(UNIT=21,FILE='disclin.all',STATUS='OLD',ERR=91) do 112 nn=1,30 read (21,*) rrr(nn), dalq(nn) 112 continue do 113 nn=1,30 read (21,*) rrrlg(nn), dalg(nn) 113 continue close(21) * goto 888 90 PRINT*, 'input - output error' 91 PRINT*, 'input - output error #2' 888 continue end ======================================================================= Ported to class by C. Loizides - 17/02/2004 */ Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx, Double_t &continuous,Double_t &discrete) const { // calculate Single Hard approx. // weights for given parton type, // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L // read-in data before first call if(!fTablesLoaded){ Error("CalcMult","Tables are not loaded."); return -1; } if(!fMultSoft){ Error("CalcMult","Tables are not loaded for Single Hard Scattering."); return -1; } Double_t rrin = rrrr; Double_t xxin = xxxx; Int_t nxlow = (Int_t)(xxin/0.038) + 1; Int_t nxhigh = nxlow + 1; Double_t xfraclow = (fxx[nxhigh-1]-xxin)/0.038; Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.038; //why this? if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD if(rrin>=frrr[0]) rrin = 0.95*frrr[0]; // AD Int_t nrlow=0,nrhigh=0; Double_t rrhigh=0,rrlow=0; for(Int_t nr=1; nr<=30; nr++) { if(rrinfgkRMax) { Warning("CalcR","Value of R = %.2f; should be less than %.2f",R,fgkRMax); return -R; } return R; } Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const { // return DeltaE for MS or SH scattering // for given parton type, length and energy // Dependant on ECM (energy constraint method) // e is used to determine where to set bins to zero. if(!fHistos){ Fatal("GetELossRandom","Call SampleEnergyLoss method before!"); return -1000.; } if((ipart<1) || (ipart>2)) { Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart); return -1000; } Int_t l=(Int_t)length; if((length-(Double_t)l)>0.5) l++; if(l<=0) return 0.; if(l>fLengthMax) l=fLengthMax; if(fECMethod==kReweight){ TH1F *dummy=new TH1F(*fHistos[ipart-1][l-1]); dummy->SetName("dummy"); for(Int_t bin=dummy->FindBin(e)+1;bin<=1100;bin++){ dummy->SetBinContent(bin,0.); } dummy->Scale(1./dummy->Integral()); Double_t ret=dummy->GetRandom(); delete dummy; return ret; //****** !! ALTERNATIVE WAY OF DOING IT !! *** //Double_t ret = 2.*e; //while(ret>e) ret=fHistos[ipart-1][l-1]->GetRandom(); //return ret; //******************************************** } else { //kDefault Double_t ret=fHistos[ipart-1][l-1]->GetRandom(); if(ret>e) return e; return ret; } } Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const { //return quenched parton energy //for given parton type, length and energy Double_t loss=GetELossRandom(ipart,length,e); return e-loss; } Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e) const { // return DeltaE for MS or SH scattering // for given parton type, length distribution and energy // Dependant on ECM (energy constraint method) // e is used to determine where to set bins to zero. if(!hell){ Warning("GetELossRandom","Pointer to length distribution is NULL."); return 0.; } Double_t ell=hell->GetRandom(); return GetELossRandom(ipart,ell,e); } Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e) const { //return quenched parton energy //for given parton type, length distribution and energy Double_t loss=GetELossRandom(ipart,hell,e); return e-loss; } Int_t AliQuenchingWeights::SampleEnergyLoss() { // Has to be called to fill the histograms // // For stored values fQTransport loop over // particle type and length = 1 to fMaxLength (fm) // to fill energy loss histos // // Take histogram of continuous weights // Take discrete_weight // If discrete_weight > 1, put all channels to 0, except channel 1 // Fill channel 1 with discrete_weight/(1-discrete_weight)*integral // read-in data before first call if(!fTablesLoaded){ Error("CalcMult","Tables are not loaded."); return -1; } if(fMultSoft) { Int_t lmax=CalcLengthMax(fQTransport); if(fLengthMax>lmax){ Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax); fLengthMax=lmax; } } else { Warning("SampleEnergyLoss","Maximum length not checked,\nbecause SingeHard is not yet tested."); } Reset(); fHistos=new TH1F**[2]; fHistos[0]=new TH1F*[fLengthMax]; fHistos[1]=new TH1F*[fLengthMax]; fLengthMaxOld=fLengthMax; //remember old value in case //user wants to reset Int_t medvalue=0; Char_t meddesc[100]; if(fMultSoft) { medvalue=(Int_t)(fQTransport*1000.); sprintf(meddesc,"MS"); } else { medvalue=(Int_t)(fMu*1000.); sprintf(meddesc,"SH"); } for(Int_t ipart=1;ipart<=2;ipart++){ for(Int_t l=1;l<=fLengthMax;l++){ Char_t hname[100]; sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l); Double_t wc = CalcWC(l); fHistos[ipart-1][l-1] = new TH1F(hname,hname,1100,0.,1.1*wc); fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]"); fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)"); fHistos[ipart-1][l-1]->SetLineColor(4); Double_t rrrr = CalcR(wc,l); Double_t discrete=0.; // loop on histogram channels for(Int_t bin=1; bin<=1100; bin++) { Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc; Double_t continuous; CalcMult(ipart,rrrr,xxxx,continuous,discrete); fHistos[ipart-1][l-1]->SetBinContent(bin,continuous); } // add discrete part to distribution if(discrete>=1.) for(Int_t bin=2;bin<=1100;bin++) fHistos[ipart-1][l-1]->SetBinContent(bin,0.); else { Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,1100); fHistos[ipart-1][l-1]->Fill(0.,val); } Double_t hint=fHistos[ipart-1][l-1]->Integral(1,1100); fHistos[ipart-1][l-1]->Scale(1./hint); } } return 0; } const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Int_t l) const { //return quenching histograms //for ipart and length if(!fHistos){ Fatal("GetELossRandom","Call SampleEnergyLoss method before!"); return 0; } if((ipart<1) || (ipart>2)) { Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart); return 0; } if(l<=0) return 0; if(l>fLengthMax) l=fLengthMax; return fHistos[ipart-1][l-1]; } TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t length) const { // ipart = 1 for quark, 2 for gluon // medval a) qtransp = transport coefficient (GeV^2/fm) // b) mu = Debye mass (GeV) // length = path length in medium (fm) // Get from SW tables: // - discrete weight (probability to have NO energy loss) // - continuous weight, as a function of dE // compute up to max dE = 1.1*wc // (as an histogram named hContQW___ // and range [0,1.1*wc] (1100 channels) Double_t wc = 0; Char_t meddesc[100]; if(fMultSoft) { wc=CalcWC(medval,length); sprintf(meddesc,"MS"); } else { wc=CalcWCbar(medval,length); sprintf(meddesc,"SH"); } Char_t hname[100]; sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart, (Int_t)(medval*1000.),(Int_t)length); TH1F *hist = new TH1F("hist",hname,1100,0.,1.1*wc); hist->SetXTitle("#Delta E [GeV]"); hist->SetYTitle("p(#Delta E)"); hist->SetLineColor(4); Double_t rrrr = CalcR(wc,length); // loop on histogram channels for(Int_t bin=1; bin<=1100; bin++) { Double_t xxxx = hist->GetBinCenter(bin)/wc; Double_t continuous,discrete; Int_t ret=0; if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete); else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete); if(ret!=0){ delete hist; return 0; }; hist->SetBinContent(bin,continuous); } return hist; } TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length) const { // ipart = 1 for quark, 2 for gluon // medval a) qtransp = transport coefficient (GeV^2/fm) // b) mu = Debye mass (GeV) // length = path length in medium (fm) // Get from SW tables: // - discrete weight (probability to have NO energy loss) // - continuous weight, as a function of dE // compute up to max dE = 1.1*wc // (as an histogram named hContQW___ // and range [0,1.1] (1100 channels) Double_t wc = 0; Char_t meddesc[100]; if(fMultSoft) { wc=CalcWC(medval,length); sprintf(meddesc,"MS"); } else { wc=CalcWCbar(medval,length); sprintf(meddesc,"SH"); } Char_t hname[100]; sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart, (Int_t)(medval*1000.),(Int_t)length); TH1F *histx = new TH1F("histx",hname,1100,0.,1.1); histx->SetXTitle("x = #Delta E/#omega_{c}"); if(fMultSoft) histx->SetYTitle("p(#Delta E/#omega_{c})"); else histx->SetYTitle("p(#Delta E/#bar#omega_{c})"); histx->SetLineColor(4); Double_t rrrr = CalcR(wc,length); // loop on histogram channels for(Int_t bin=1; bin<=1100; bin++) { Double_t xxxx = histx->GetBinCenter(bin); Double_t continuous,discrete; Int_t ret=0; if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete); else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete); if(ret!=0){ delete histx; return 0; }; histx->SetBinContent(bin,continuous); } return histx; } TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const { //compute energy loss histogram for //parton type, medium value, length and energy AliQuenchingWeights *dummy=new AliQuenchingWeights(*this); if(fMultSoft){ dummy->SetQTransport(medval); dummy->InitMult(); } else { dummy->SetMu(medval); dummy->InitSingleHard(); } dummy->SampleEnergyLoss(); Char_t name[100]; Char_t hname[100]; if(ipart==1){ sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); sprintf(hname,"hLossQuarks"); } else { sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); sprintf(hname,"hLossGluons"); } TH1F *h = new TH1F(hname,name,250,0,250); for(Int_t i=0;i<100000;i++){ //if(i % 1000 == 0) cout << "." << flush; Double_t loss=dummy->GetELossRandom(ipart,l,e); h->Fill(loss); } h->SetStats(kTRUE); delete dummy; return h; } TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const { //compute energy loss histogram for //parton type, medium value, //length distribution and energy AliQuenchingWeights *dummy=new AliQuenchingWeights(*this); if(fMultSoft){ dummy->SetQTransport(medval); dummy->InitMult(); } else { dummy->SetMu(medval); dummy->InitSingleHard(); } dummy->SampleEnergyLoss(); Char_t name[100]; Char_t hname[100]; if(ipart==1){ sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); sprintf(hname,"hLossQuarks"); } else { sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); sprintf(hname,"hLossGluons"); } TH1F *h = new TH1F(hname,name,250,0,250); for(Int_t i=0;i<100000;i++){ //if(i % 1000 == 0) cout << "." << flush; Double_t loss=dummy->GetELossRandom(ipart,hEll,e); h->Fill(loss); } h->SetStats(kTRUE); delete dummy; return h; } void AliQuenchingWeights::PlotDiscreteWeights(Int_t len) const { //plot discrete weights for given length TCanvas *c; if(fMultSoft) c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400); else c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400); c->cd(); TH2F *hframe = new TH2F("hdisc","",2,0,5.1,2,0,1); hframe->SetStats(0); if(fMultSoft) hframe->SetXTitle("#hat{q} [GeV^{2}/fm]"); else hframe->SetXTitle("#mu [GeV]"); hframe->SetYTitle("Probability #Delta E = 0 , p_{0}"); hframe->Draw(); TGraph *gq=new TGraph(20); Int_t i=0; if(fMultSoft) { for(Double_t q=0.05;q<=5.05;q+=0.25){ Double_t disc,cont; CalcMult(1,1.0, q, len, cont, disc); gq->SetPoint(i,q,disc);i++; } } else { for(Double_t m=0.05;m<=5.05;m+=0.25){ Double_t disc,cont; CalcSingleHard(1,1.0, m, len, cont, disc); gq->SetPoint(i,m,disc);i++; } } gq->SetMarkerStyle(20); gq->Draw("pl"); TGraph *gg=new TGraph(20); i=0; if(fMultSoft){ for(Double_t q=0.05;q<=5.05;q+=0.25){ Double_t disc,cont; CalcMult(2,1.0, q, 5., cont, disc); gg->SetPoint(i,q,disc);i++; } } else { for(Double_t m=0.05;m<=5.05;m+=0.25){ Double_t disc,cont; CalcSingleHard(2,1.0, m, 5., cont, disc); gg->SetPoint(i,m,disc);i++; } } gg->SetMarkerStyle(24); gg->Draw("pl"); TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; sprintf(label,"L = %d fm",len); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark","pl"); l1a->AddEntry(gg,"gluon","pl"); l1a->Draw(); c->Update(); } void AliQuenchingWeights::PlotContWeights(Int_t itype,Int_t ell) const { //plot continous weights for //given parton type and length Float_t medvals[3]; Char_t title[1024]; Char_t name[1024]; if(fMultSoft) { if(itype==1) sprintf(title,"Cont. Weight for Multiple Scattering - Quarks"); else if(itype==2) sprintf(title,"Cont. Weight for Multiple Scattering - Gluons"); else return; medvals[0]=4;medvals[1]=1;medvals[2]=0.5; sprintf(name,"ccont-ms-%d",itype); } else { if(itype==1) sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks"); else if(itype==2) sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons"); else return; medvals[0]=2;medvals[1]=1;medvals[2]=0.5; sprintf(name,"ccont-ms-%d",itype); } TCanvas *c = new TCanvas(name,title,0,0,500,400); c->cd(); TH1F *h1=ComputeQWHisto(itype,medvals[0],ell); h1->SetName("h1"); h1->SetTitle(title); h1->SetStats(0); h1->SetLineColor(1); h1->DrawCopy(); TH1F *h2=ComputeQWHisto(itype,medvals[1],ell); h2->SetName("h2"); h2->SetLineColor(2); h2->DrawCopy("SAME"); TH1F *h3=ComputeQWHisto(itype,medvals[2],ell); h3->SetName("h3"); h3->SetLineColor(3); h3->DrawCopy("SAME"); TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; sprintf(label,"L = %d fm",ell); l1a->AddEntry(h1,label,""); if(fMultSoft) { sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]); l1a->AddEntry(h1,label,"pl"); sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]); l1a->AddEntry(h2,label,"pl"); sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]); l1a->AddEntry(h3,label,"pl"); } else { sprintf(label,"#mu = %.1f GeV",medvals[0]); l1a->AddEntry(h1,label,"pl"); sprintf(label,"#mu = %.1f GeV",medvals[1]); l1a->AddEntry(h2,label,"pl"); sprintf(label,"#mu = %.1f GeV",medvals[2]); l1a->AddEntry(h3,label,"pl"); } l1a->Draw(); c->Update(); } void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const { //plot continous weights for //given parton type and medium value Char_t title[1024]; Char_t name[1024]; if(fMultSoft) { if(itype==1) sprintf(title,"Cont. Weight for Multiple Scattering - Quarks"); else if(itype==2) sprintf(title,"Cont. Weight for Multiple Scattering - Gluons"); else return; sprintf(name,"ccont2-ms-%d",itype); } else { if(itype==1) sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks"); else if(itype==2) sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons"); else return; sprintf(name,"ccont2-sh-%d",itype); } TCanvas *c = new TCanvas(name,title,0,0,500,400); c->cd(); TH1F *h1=ComputeQWHisto(itype,medval,8); h1->SetName("h1"); h1->SetTitle(title); h1->SetStats(0); h1->SetLineColor(1); h1->DrawCopy(); TH1F *h2=ComputeQWHisto(itype,medval,5); h2->SetName("h2"); h2->SetLineColor(2); h2->DrawCopy("SAME"); TH1F *h3=ComputeQWHisto(itype,medval,2); h3->SetName("h3"); h3->SetLineColor(3); h3->DrawCopy("SAME"); TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; if(fMultSoft) sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval); else sprintf(label,"#mu = %.1f GeV",medval); l1a->AddEntry(h1,label,""); l1a->AddEntry(h1,"L = 8 fm","pl"); l1a->AddEntry(h2,"L = 5 fm","pl"); l1a->AddEntry(h3,"L = 2 fm","pl"); l1a->Draw(); c->Update(); } void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e) const { //plot average energy loss for given length //and parton energy if(!fTablesLoaded){ Error("CalcMult","Tables are not loaded."); return; } Char_t title[1024]; Char_t name[1024]; if(fMultSoft){ sprintf(title,"Average Loss for Multiple Scattering"); sprintf(name,"cavgelossms"); } else { sprintf(title,"Average Loss for Single Hard Scattering"); sprintf(name,"cavgelosssh"); } TCanvas *c = new TCanvas(name,title,0,0,500,400); c->cd(); TH2F *hframe = new TH2F("avgloss",title,2,0,5.1,2,0,100); hframe->SetStats(0); if(fMultSoft) hframe->SetXTitle("#hat{q} [GeV^{2}/fm]"); else hframe->SetXTitle("#mu [GeV]"); hframe->SetYTitle(" [GeV]"); hframe->Draw(); TGraph *gq=new TGraph(20); Int_t i=0; for(Double_t v=0.05;v<=5.05;v+=0.25){ TH1F *dummy=ComputeELossHisto(1,v,len,e); Double_t avgloss=dummy->GetMean(); gq->SetPoint(i,v,avgloss);i++; delete dummy; } gq->SetMarkerStyle(20); gq->Draw("pl"); TGraph *gg=new TGraph(20); i=0; for(Double_t v=0.05;v<=5.05;v+=0.25){ TH1F *dummy=ComputeELossHisto(2,v,len,e); Double_t avgloss=dummy->GetMean(); gg->SetPoint(i,v,avgloss);i++; delete dummy; } gg->SetMarkerStyle(24); gg->Draw("pl"); TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; sprintf(label,"L = %d fm",len); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark","pl"); l1a->AddEntry(gg,"gluon","pl"); l1a->Draw(); c->Update(); } void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const { //plot average energy loss for given //length distribution and parton energy if(!fTablesLoaded){ Error("CalcMult","Tables are not loaded."); return; } Char_t title[1024]; Char_t name[1024]; if(fMultSoft){ sprintf(title,"Average Loss for Multiple Scattering"); sprintf(name,"cavgelossms2"); } else { sprintf(title,"Average Loss for Single Hard Scattering"); sprintf(name,"cavgelosssh2"); } TCanvas *c = new TCanvas(name,title,0,0,500,400); c->cd(); TH2F *hframe = new TH2F("havgloss",title,2,0,5.1,2,0,100); hframe->SetStats(0); if(fMultSoft) hframe->SetXTitle("#hat{q} [GeV^{2}/fm]"); else hframe->SetXTitle("#mu [GeV]"); hframe->SetYTitle(" [GeV]"); hframe->Draw(); TGraph *gq=new TGraph(20); Int_t i=0; for(Double_t v=0.05;v<=5.05;v+=0.25){ TH1F *dummy=ComputeELossHisto(1,v,hEll,e); Double_t avgloss=dummy->GetMean(); gq->SetPoint(i,v,avgloss);i++; delete dummy; } gq->SetMarkerStyle(20); gq->Draw("pl"); TGraph *gg=new TGraph(20); i=0; for(Double_t v=0.05;v<=5.05;v+=0.25){ TH1F *dummy=ComputeELossHisto(2,v,hEll,e); Double_t avgloss=dummy->GetMean(); gg->SetPoint(i,v,avgloss);i++; delete dummy; } gg->SetMarkerStyle(24); gg->Draw("pl"); TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; sprintf(label," = %.2f fm",hEll->GetMean()); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark","pl"); l1a->AddEntry(gg,"gluon","pl"); l1a->Draw(); c->Update(); } void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Int_t len) const { //plot relative energy loss for given //length and parton energy versus pt if(!fTablesLoaded){ Error("CalcMult","Tables are not loaded."); return; } Char_t title[1024]; Char_t name[1024]; if(fMultSoft){ sprintf(title,"Relative Loss for Multiple Scattering"); sprintf(name,"cavgelossvsptms"); } else { sprintf(title,"Relative Loss for Single Hard Scattering"); sprintf(name,"cavgelossvsptsh"); } TCanvas *c = new TCanvas(name,title,0,0,500,400); c->cd(); TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1); hframe->SetStats(0); hframe->SetXTitle("p_{T} [GeV]"); hframe->SetYTitle("/p_{T} [GeV]"); hframe->Draw(); TGraph *gq=new TGraph(40); Int_t i=0; for(Double_t pt=2.5;pt<=100.05;pt+=2.5){ TH1F *dummy=ComputeELossHisto(1,medval,len,pt); Double_t avgloss=dummy->GetMean(); gq->SetPoint(i,pt,avgloss/pt);i++; delete dummy; } gq->SetMarkerStyle(20); gq->Draw("pl"); TGraph *gg=new TGraph(40); i=0; for(Double_t pt=2.5;pt<=100.05;pt+=2.5){ TH1F *dummy=ComputeELossHisto(2,medval,len,pt); Double_t avgloss=dummy->GetMean(); gg->SetPoint(i,pt,avgloss/pt);i++; delete dummy; } gg->SetMarkerStyle(24); gg->Draw("pl"); TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; sprintf(label,"L = %d fm",len); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark","pl"); l1a->AddEntry(gg,"gluon","pl"); l1a->Draw(); c->Update(); } void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const { //plot relative energy loss for given //length distribution and parton energy versus pt if(!fTablesLoaded){ Error("CalcMult","Tables are not loaded."); return; } Char_t title[1024]; Char_t name[1024]; if(fMultSoft){ sprintf(title,"Relative Loss for Multiple Scattering"); sprintf(name,"cavgelossvsptms2"); } else { sprintf(title,"Relative Loss for Single Hard Scattering"); sprintf(name,"cavgelossvsptsh2"); } TCanvas *c = new TCanvas(name,title,0,0,500,400); c->cd(); TH2F *hframe = new TH2F("havglossvspt",title,2,0,100,2,0,1); hframe->SetStats(0); hframe->SetXTitle("p_{T} [GeV]"); hframe->SetYTitle("/p_{T} [GeV]"); hframe->Draw(); TGraph *gq=new TGraph(40); Int_t i=0; for(Double_t pt=2.5;pt<=100.05;pt+=2.5){ TH1F *dummy=ComputeELossHisto(1,medval,hEll,pt); Double_t avgloss=dummy->GetMean(); gq->SetPoint(i,pt,avgloss/pt);i++; delete dummy; } gq->SetMarkerStyle(20); gq->Draw("pl"); TGraph *gg=new TGraph(40); i=0; for(Double_t pt=2.5;pt<=100.05;pt+=2.5){ TH1F *dummy=ComputeELossHisto(2,medval,hEll,pt); Double_t avgloss=dummy->GetMean(); gg->SetPoint(i,pt,avgloss/pt);i++; delete dummy; } gg->SetMarkerStyle(24); gg->Draw("pl"); TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; sprintf(label," = %.2f fm",hEll->GetMean()); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark","pl"); l1a->AddEntry(gg,"gluon","pl"); l1a->Draw(); c->Update(); }