/************************************************************************** * 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 // //=================== Added by C. Loizides 27/03/04 =========================== // // Added support for k-Quenching, where wc=I1*k and R=2I1^2/I0*k // (see the AliFastGlauber class for definition of I0/I1) //----------------------------------------------------------------------------- #include #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; // hist binning const Int_t AliQuenchingWeights::fgkBins = 1300; const Double_t AliQuenchingWeights::fgkMaxBin = 1.3; // counter for histogram labels Int_t AliQuenchingWeights::fgCounter = 0; AliQuenchingWeights::AliQuenchingWeights() : TObject(), fInstanceNumber(fgCounter++), fMultSoft(kTRUE), fECMethod(kDefault), fQTransport(1.), fMu(1.), fK(4.e5), fLengthMax(20), fLengthMaxOld(0), fHistos(0), fHisto(0), fTablesLoaded(kFALSE) { //default constructor Char_t name[100]; sprintf(name,"hhistoqw_%d",fInstanceNumber); fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin); for(Int_t bin=1;bin<=fgkBins;bin++) fHisto->SetBinContent(bin,0.); } AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a) : TObject(), fInstanceNumber(fgCounter++), fMultSoft(kTRUE), fECMethod(kDefault), fQTransport(1.), fMu(1.), fK(4.e5), fLengthMax(20), fLengthMaxOld(0), fHistos(0), fHisto(0), fTablesLoaded(kFALSE) { // copy constructor fTablesLoaded=kFALSE; fHistos=0; fLengthMaxOld=0; fMultSoft=a.GetMultSoft();; fMu=a.GetMu(); fK=a.GetK(); fQTransport=a.GetQTransport(); fECMethod=(kECMethod)a.GetECMethod(); fLengthMax=a.GetLengthMax(); fInstanceNumber=fgCounter++; Char_t name[100]; sprintf(name,"hhistoqw_%d",fInstanceNumber); fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin); for(Int_t bin=1;bin<=fgkBins;bin++) fHisto->SetBinContent(bin,0.); //Missing in the class is the pathname //to the tables, can be added if needed } AliQuenchingWeights::~AliQuenchingWeights() { Reset(); delete fHisto; } void AliQuenchingWeights::Reset() { //reset tables if there were used if(!fHistos) return; for(Int_t l=0;l<4*fLengthMaxOld;l++){ delete fHistos[0][l]; delete fHistos[1][l]; } delete[] fHistos; fHistos=0; fLengthMaxOld=0; } void AliQuenchingWeights::SetECMethod(kECMethod type) { //set energy constraint method fECMethod=type; if(fECMethod==kDefault) Info("SetECMethod","Energy Constraint Method set to DEFAULT:\nIf (sampled energy loss > parton energy) then sampled energy loss = parton energy."); else if(fECMethod==kReweight) Info("SetECMethod","Energy Constraint Method set to REWEIGHT:\nRequire sampled energy loss <= parton energy."); else Info("SetECMethod","Energy Constraint Method set to REWEIGHTCONT:\nRequire sampled energy loss <= parton energy (only implemented for FAST method."); } 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("CalcSingleHard","Tables are not loaded."); return -1; } if(fMultSoft){ Error("CalcSingleHard","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(rrin= fgkRMax) { Warning("CalcR","Value of r = %.2f; should be less than %.2f", r, fgkRMax); return fgkRMax-1; } return r; } Double_t AliQuenchingWeights::CalcRk(Double_t k, Double_t I0, Double_t I1) const { //calculate R value and //check if it is less then maximum Double_t r = fgkRMax-1; if(I0>0) r = 2*I1*I1/I0*k; if(r>=fgkRMax) { Warning("CalcRk","Value of r = %.2f; should be less than %.2f",r,fgkRMax); return fgkRMax-1; } 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=GetIndex(length); if(l<=0) return 0.; if(fECMethod==kReweight){ Double_t ret = 2.*e; Int_t ws=0; while(ret>e){ ret=fHistos[ipart-1][l-1]->GetRandom(); if(++ws==1e6){ Warning("GetELossRandom", "Aborted reweighting; maximum loss assigned after 1e6 trials."); return e; } } return ret; } //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; } Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t I1, Double_t e) { // return DeltaE for new dynamic version // for given parton type, I0 and I1 value and energy // Dependant on ECM (energy constraint method) // e is used to determine where to set bins to zero. // read-in data before first call if(!fTablesLoaded){ Fatal("GetELossRandomK","Tables are not loaded."); return -1000.; } if((ipart<1) || (ipart>2)) { Fatal("GetELossRandomK","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart); return -1000.; } Double_t r=CalcRk(I0,I1); if(r<0.){ Fatal("GetELossRandomK","R should not be negative"); return -1000.; } Double_t wc=CalcWCk(I1); if(wc<=0.){ Fatal("GetELossRandomK","wc should be greater than zero"); return -1000.; } if(SampleEnergyLoss(ipart,r)!=0){ Fatal("GetELossRandomK","Could not sample energy loss"); return -1000.; } if(fECMethod==kReweight){ Double_t ret = 2.*e; Int_t ws=0; while(ret>e){ ret=fHisto->GetRandom(); if(++ws==1e6){ Warning("GetELossRandomK", "Aborted reweighting; maximum loss assigned after 1e6 trials."); return e; } } return ret; } //kDefault Double_t ret=fHisto->GetRandom()*wc; if(ret>e) return e; return ret; } Double_t AliQuenchingWeights::CalcQuenchedEnergyK(Int_t ipart, Double_t I0, Double_t I1, Double_t e) { //return quenched parton energy //for given parton type, I0 and I1 value and energy Double_t loss=GetELossRandomK(ipart,I0,I1,e); return e-loss; } Double_t AliQuenchingWeights::GetELossRandomKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e) { // return DeltaE for new dynamic version // for given parton type, I0 and I1 value and energy // Dependant on ECM (energy constraint method) // e is used to determine where to set bins to zero. // method is optimized and should only be used if // all parameters are well within the bounds. // read-in data tables before first call Double_t r=CalcRk(I0,I1); if(r<=0.){ return 0.; } Double_t wc=CalcWCk(I1); if(wc<=0.){ return 0.; } return GetELossRandomKFastR(ipart,r,wc,e); } Double_t AliQuenchingWeights::GetELossRandomKFastR(Int_t ipart, Double_t r, Double_t wc, Double_t e) { // return DeltaE for new dynamic version // for given parton type, R and wc value and energy // Dependant on ECM (energy constraint method) // e is used to determine where to set bins to zero. // method is optimized and should only be used if // all parameters are well within the bounds. // read-in data tables before first call if(r>=fgkRMax) { r=fgkRMax-1; } Double_t discrete=0.; Double_t continuous=0.; Int_t bin=1; Double_t xxxx = fHisto->GetBinCenter(bin); if(fMultSoft) CalcMult(ipart,r,xxxx,continuous,discrete); else CalcSingleHard(ipart,r,xxxx,continuous,discrete); if(discrete>=1.0) { return 0.; //no energy loss } fHisto->SetBinContent(bin,continuous); Int_t kbinmax=fHisto->FindBin(e/wc); if(kbinmax>=fgkBins) kbinmax=fgkBins-1; if(kbinmax==1) return e; //maximum energy loss if(fMultSoft) { for(bin=2; bin<=kbinmax; bin++) { xxxx = fHisto->GetBinCenter(bin); CalcMult(ipart,r,xxxx,continuous,discrete); fHisto->SetBinContent(bin,continuous); } } else { for(bin=2; bin<=kbinmax; bin++) { xxxx = fHisto->GetBinCenter(bin); CalcSingleHard(ipart,r,xxxx,continuous,discrete); fHisto->SetBinContent(bin,continuous); } } if(fECMethod==kReweight){ fHisto->SetBinContent(kbinmax+1,0); fHisto->Fill(0.,discrete*fgkBins/fgkMaxBin); } else if (fECMethod==kReweightCont) { fHisto->SetBinContent(kbinmax+1,0); const Double_t kdelta=fHisto->Integral(1,kbinmax); fHisto->Scale(1./kdelta*(1-discrete)); fHisto->Fill(0.,discrete); } else { const Double_t kdelta=fHisto->Integral(1,kbinmax); Double_t val=discrete*fgkBins/fgkMaxBin; fHisto->Fill(0.,val); fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta); } for(bin=kbinmax+2; bin<=fgkBins; bin++) { fHisto->SetBinContent(bin,0); } //cout << kbinmax << " " << discrete << " " << fHisto->Integral() << endl; Double_t ret=fHisto->GetRandom()*wc; if(ret>e) return e; return ret; } Double_t AliQuenchingWeights::CalcQuenchedEnergyKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e) { //return quenched parton energy (fast method) //for given parton type, I0 and I1 value and energy Double_t loss=GetELossRandomKFast(ipart,I0,I1,e); return e-loss; } Double_t AliQuenchingWeights::GetDiscreteWeight(Int_t ipart, Double_t I0, Double_t I1) { // return discrete weight Double_t r=CalcRk(I0,I1); if(r<=0.){ return 1.; } return GetDiscreteWeightR(ipart,r); } Double_t AliQuenchingWeights::GetDiscreteWeightR(Int_t ipart, Double_t r) { // return discrete weight if(r>=fgkRMax) { r=fgkRMax-1; } Double_t discrete=0.; Double_t continuous=0.; Int_t bin=1; Double_t xxxx = fHisto->GetBinCenter(bin); if(fMultSoft) CalcMult(ipart,r,xxxx,continuous,discrete); else CalcSingleHard(ipart,r,xxxx,continuous,discrete); return discrete; } void AliQuenchingWeights::GetZeroLossProb(Double_t &p,Double_t &prw,Double_t &prwcont, Int_t ipart,Double_t I0,Double_t I1,Double_t e) { //calculate the probabilty that there is no energy //loss for different ways of energy constraint p=1.;prw=1.;prwcont=1.; Double_t r=CalcRk(I0,I1); if(r<=0.){ return; } Double_t wc=CalcWCk(I1); if(wc<=0.){ return; } GetZeroLossProbR(p,prw,prwcont,ipart,r,wc,e); } void AliQuenchingWeights::GetZeroLossProbR(Double_t &p,Double_t &prw,Double_t &prwcont, Int_t ipart, Double_t r,Double_t wc,Double_t e) { //calculate the probabilty that there is no energy //loss for different ways of energy constraint if(r>=fgkRMax) { r=fgkRMax-1; } Double_t discrete=0.; Double_t continuous=0.; Int_t kbinmax=fHisto->FindBin(e/wc); if(kbinmax>=fgkBins) kbinmax=fgkBins-1; if(fMultSoft) { for(Int_t bin=1; bin<=kbinmax; bin++) { Double_t xxxx = fHisto->GetBinCenter(bin); CalcMult(ipart,r,xxxx,continuous,discrete); fHisto->SetBinContent(bin,continuous); } } else { for(Int_t bin=1; bin<=kbinmax; bin++) { Double_t xxxx = fHisto->GetBinCenter(bin); CalcSingleHard(ipart,r,xxxx,continuous,discrete); fHisto->SetBinContent(bin,continuous); } } //non-reweighted P(Delta E = 0) const Double_t kdelta=fHisto->Integral(1,kbinmax); Double_t val=discrete*fgkBins/fgkMaxBin; fHisto->Fill(0.,val); fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta); Double_t hint=fHisto->Integral(1,kbinmax+1); p=fHisto->GetBinContent(1)/hint; // reweighted hint=fHisto->Integral(1,kbinmax); prw=fHisto->GetBinContent(1)/hint; Double_t xxxx = fHisto->GetBinCenter(1); CalcMult(ipart,r,xxxx,continuous,discrete); fHisto->SetBinContent(1,continuous); hint=fHisto->Integral(1,kbinmax); fHisto->Scale(1./hint*(1-discrete)); fHisto->Fill(0.,discrete); prwcont=fHisto->GetBinContent(1); } 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("SampleEnergyLoss","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*[4*fLengthMax]; fHistos[1]=new TH1F*[4*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<=4*fLengthMax;l++){ Char_t hname[100]; sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l); Double_t len=l/4.; Double_t wc = CalcWC(len); fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgkBins,0.,fgkMaxBin*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,len); Double_t discrete=0.; // loop on histogram channels for(Int_t bin=1; bin<=fgkBins; bin++) { Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc; Double_t continuous; if(fMultSoft) CalcMult(ipart,rrrr,xxxx,continuous,discrete); else CalcSingleHard(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<=fgkBins;bin++) fHistos[ipart-1][l-1]->SetBinContent(bin,0.); else { Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgkBins); fHistos[ipart-1][l-1]->Fill(0.,val); } Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins); fHistos[ipart-1][l-1]->Scale(1./hint); } } return 0; } Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t r) { // Sample energy loss directly for one particle type // choses R (safe it and keep it until next call of function) // read-in data before first call if(!fTablesLoaded){ Error("SampleEnergyLoss","Tables are not loaded."); return -1; } Double_t discrete=0.; Double_t continuous=0;; Int_t bin=1; Double_t xxxx = fHisto->GetBinCenter(bin); if(fMultSoft) CalcMult(ipart,r,xxxx,continuous,discrete); else CalcSingleHard(ipart,r,xxxx,continuous,discrete); if(discrete>=1.) { fHisto->SetBinContent(1,1.); for(bin=2;bin<=fgkBins;bin++) fHisto->SetBinContent(bin,0.); return 0; } fHisto->SetBinContent(bin,continuous); for(bin=2; bin<=fgkBins; bin++) { xxxx = fHisto->GetBinCenter(bin); if(fMultSoft) CalcMult(ipart,r,xxxx,continuous,discrete); else CalcSingleHard(ipart,r,xxxx,continuous,discrete); fHisto->SetBinContent(bin,continuous); } Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins); fHisto->Fill(0.,val); Double_t hint=fHisto->Integral(1,fgkBins); if(hint!=0) fHisto->Scale(1./hint); else { //cout << discrete << " " << hint << " " << continuous << endl; return -1; } return 0; } const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) 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; } Int_t l=GetIndex(length); if(l<=0) return 0; 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: // - continuous weight, as a function of dE/wc 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,fgkBins,0.,fgkMaxBin*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<=fgkBins; 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: // - continuous weight, as a function of dE/wc 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,fgkBins,0.,fgkMaxBin); 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<=fgkBins; 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::ComputeQWHistoX(Int_t ipart,Double_t r) const { // compute P(E) distribution for // given ipart = 1 for quark, 2 for gluon // and R Char_t meddesc[100]; if(fMultSoft) { sprintf(meddesc,"MS"); } else { sprintf(meddesc,"SH"); } Char_t hname[100]; sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,r); TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin); 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 = r; Double_t continuous=0.,discrete=0.; //loop on histogram channels for(Int_t bin=1; bin<=fgkBins; bin++) { Double_t xxxx = histx->GetBinCenter(bin); 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); } //add discrete part to distribution if(discrete>=1.) for(Int_t bin=2;bin<=fgkBins;bin++) histx->SetBinContent(bin,0.); else { Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins); histx->Fill(0.,val); } Double_t hint=histx->Integral(1,fgkBins); if(hint!=0) histx->Scale(1./hint); 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; } TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t r) const { // compute energy loss histogram for // parton type and given R TH1F *dummy = ComputeQWHistoX(ipart,r); if(!dummy) return 0; Char_t hname[100]; sprintf(hname,"hELossHistox_%d_%.2f",ipart,r); TH1F *histx = new TH1F("histxr",hname,fgkBins,0.,fgkMaxBin); for(Int_t i=0;i<100000;i++){ //if(i % 1000 == 0) cout << "." << flush; Double_t loss=dummy->GetRandom(); histx->Fill(loss); } delete dummy; return histx; } Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,Double_t l) const { // compute average energy loss for // parton type, medium value, length and energy TH1F *dummy = ComputeELossHisto(ipart,medval,l); if(!dummy) return 0; Double_t ret=dummy->GetMean(); delete dummy; return ret; } Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEll) const { // compute average energy loss for // parton type, medium value, // length distribution and energy TH1F *dummy = ComputeELossHisto(ipart,medval,hEll); if(!dummy) return 0; Double_t ret=dummy->GetMean(); delete dummy; return ret; } Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t r) const { // compute average energy loss over wc // for parton type and given R TH1F *dummy = ComputeELossHisto(ipart,r); if(!dummy) return 0; Double_t ret=dummy->GetMean(); delete dummy; return ret; } void AliQuenchingWeights::PlotDiscreteWeights(Double_t len,Double_t qm) 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,qm+.1,2,0,1.25); 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->SetYTitle("p_{0} (discrete weight)"); hframe->Draw(); Int_t points=(Int_t)qm*4; TGraph *gq=new TGraph(points); Int_t i=0; if(fMultSoft) { for(Double_t q=0.05;q<=qm+.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<=qm+.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->SetMarkerColor(1); gq->SetLineStyle(1); gq->SetLineColor(1); gq->Draw("l"); TGraph *gg=new TGraph(points); i=0; if(fMultSoft){ for(Double_t q=0.05;q<=qm+.05;q+=0.25){ Double_t disc,cont; CalcMult(2,1.0,q,len,cont,disc); gg->SetPoint(i,q,disc);i++; } } else { for(Double_t m=0.05;m<=qm+.05;m+=0.25){ Double_t disc,cont; CalcSingleHard(2,1.0,m,len,cont,disc); gg->SetPoint(i,m,disc);i++; } } gg->SetMarkerStyle(24); gg->SetMarkerColor(2); gg->SetLineStyle(2); gg->SetLineColor(2); gg->Draw("l"); TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; sprintf(label,"L = %.1f fm",len); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark projectile","l"); l1a->AddEntry(gg,"gluon projectile","l"); l1a->Draw(); c->Update(); } void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_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 = %.1f 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::PlotContWeightsVsL(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(Double_t len,Double_t qm,Double_t e) const { // plot average energy loss for given length // and parton energy if(!fTablesLoaded){ Error("PlotAvgELoss","Tables are not loaded."); return; } Char_t title[1024]; Char_t name[1024]; if(fMultSoft){ sprintf(title,"Average Energy Loss for Multiple Scattering"); sprintf(name,"cavgelossms"); } else { sprintf(title,"Average Energy 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","",2,0,qm+.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<=qm+.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(21); gq->Draw("pl"); Int_t points=(Int_t)qm*4; TGraph *gg=new TGraph(points); i=0; for(Double_t v=0.05;v<=qm+.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(20); gg->SetMarkerColor(2); gg->Draw("pl"); TGraph *gratio=new TGraph(points); for(i=0;iGetPoint(i,x,y); gq->GetPoint(i,x2,y2); if(y2>0) gratio->SetPoint(i,x,y/y2*10/2.25); else gratio->SetPoint(i,x,0); } gratio->SetLineStyle(4); gratio->Draw(); TLegend *l1a = new TLegend(0.15,0.60,0.50,0.90); l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; sprintf(label,"L = %.1f fm",len); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark projectile","pl"); l1a->AddEntry(gg,"gluon projectile","pl"); l1a->AddEntry(gratio,"gluon/quark/2.25*10","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("PlotAvgELossVs","Tables are not loaded."); return; } Char_t title[1024]; Char_t name[1024]; if(fMultSoft){ sprintf(title,"Average Energy Loss for Multiple Scattering"); sprintf(name,"cavgelossms2"); } else { sprintf(title,"Average Energy 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"); TGraph *gratio=new TGraph(20); for(i=0;i<20;i++){ Double_t x,y,x2,y2; gg->GetPoint(i,x,y); gq->GetPoint(i,x2,y2); if(y2>0) gratio->SetPoint(i,x,y/y2*10/2.25); else gratio->SetPoint(i,x,0); } gratio->SetLineStyle(4); //gratio->Draw(); 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->AddEntry(gratio,"gluon/quark/2.25*10","pl"); l1a->Draw(); c->Update(); } void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const { // plot average energy loss versus ell if(!fTablesLoaded){ Error("PlotAvgELossVsEll","Tables are not loaded."); return; } Char_t title[1024]; Char_t name[1024]; Float_t medval; if(fMultSoft){ sprintf(title,"Average Energy Loss for Multiple Scattering"); sprintf(name,"cavgelosslms"); medval=fQTransport; } else { sprintf(title,"Average Energy Loss for Single Hard Scattering"); sprintf(name,"cavgelosslsh"); medval=fMu; } TCanvas *c = new TCanvas(name,title,0,0,600,400); c->cd(); TH2F *hframe = new TH2F("avglossell",title,2,0,fLengthMax,2,0,250); hframe->SetStats(0); hframe->SetXTitle("length [fm]"); hframe->SetYTitle(" [GeV]"); hframe->Draw(); TGraph *gq=new TGraph((Int_t)fLengthMax*4); Int_t i=0; for(Double_t v=0.25;v<=fLengthMax;v+=0.25){ TH1F *dummy=ComputeELossHisto(1,medval,v,e); Double_t avgloss=dummy->GetMean(); gq->SetPoint(i,v,avgloss);i++; delete dummy; } gq->SetMarkerStyle(20); gq->Draw("pl"); TGraph *gg=new TGraph((Int_t)fLengthMax*4); i=0; for(Double_t v=0.25;v<=fLengthMax;v+=0.25){ TH1F *dummy=ComputeELossHisto(2,medval,v,e); Double_t avgloss=dummy->GetMean(); gg->SetPoint(i,v,avgloss);i++; delete dummy; } gg->SetMarkerStyle(24); gg->Draw("pl"); TGraph *gratio=new TGraph((Int_t)fLengthMax*4); for(i=0;i<=(Int_t)fLengthMax*4;i++){ Double_t x,y,x2,y2; gg->GetPoint(i,x,y); gq->GetPoint(i,x2,y2); if(y2>0) gratio->SetPoint(i,x,y/y2*100/2.25); else gratio->SetPoint(i,x,0); } gratio->SetLineStyle(1); gratio->SetLineWidth(2); gratio->Draw(); TLegend *l1a = new TLegend(0.15,0.65,.60,0.85); l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; if(fMultSoft) sprintf(label,"#hat{q} = %.2f GeV^{2}/fm",medval); else sprintf(label,"#mu = %.2f GeV",medval); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark","pl"); l1a->AddEntry(gg,"gluon","pl"); l1a->AddEntry(gratio,"gluon/quark/2.25*100","pl"); l1a->Draw(); TF1 *f=new TF1("f","100",0,fLengthMax); f->SetLineStyle(4); f->SetLineWidth(1); f->Draw("same"); c->Update(); } void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const { // plot relative energy loss for given // length and parton energy versus pt if(!fTablesLoaded){ Error("PlotAvgELossVsPt","Tables are not loaded."); return; } Char_t title[1024]; Char_t name[1024]; if(fMultSoft){ sprintf(title,"Relative Energy Loss for Multiple Scattering"); sprintf(name,"cavgelossvsptms"); } else { sprintf(title,"Relative Energy 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 = %.1f 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("PlotAvgELossVsPt","Tables are not loaded."); return; } Char_t title[1024]; Char_t name[1024]; if(fMultSoft){ sprintf(title,"Relative Energy Loss for Multiple Scattering"); sprintf(name,"cavgelossvsptms2"); } else { sprintf(title,"Relative Energy 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(); } Int_t AliQuenchingWeights::GetIndex(Double_t len) const { //get the index according to length if(len>fLengthMax) len=fLengthMax; Int_t l=Int_t(len/0.25); if((len-l*0.25)>0.125) l++; return l; }