X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=FASTSIM%2FAliQuenchingWeights.cxx;h=c3c1057fce5535c28a81cfa0aab70d5d442a12ee;hb=fcff8ce7274128a4975f5d9001ca3fad89abc0fc;hp=0a582fc26d39ede17d22abeb7d58297d79e6be27;hpb=e99e3ed51a94252048f99ad6659cdcc4c45ad72c;p=u%2Fmrichter%2FAliRoot.git diff --git a/FASTSIM/AliQuenchingWeights.cxx b/FASTSIM/AliQuenchingWeights.cxx index 0a582fc26d3..c3c1057fce5 100644 --- a/FASTSIM/AliQuenchingWeights.cxx +++ b/FASTSIM/AliQuenchingWeights.cxx @@ -15,65 +15,103 @@ /* $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 + +// 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] // -// 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 #include "AliQuenchingWeights.h" +using std::fstream; +using std::ios; ClassImp(AliQuenchingWeights) // conversion from fm to GeV^-1: 1 fm = fmGeV GeV^-1 -const Double_t AliQuenchingWeights::gkConvFmToInvGeV = 1./0.197; +const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197; // maximum value of R -const Double_t AliQuenchingWeights::gkRMax = 40000.; +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::gCounter = 0; +Int_t AliQuenchingWeights::fgCounter = 0; + AliQuenchingWeights::AliQuenchingWeights() - : TObject() + : TObject(), + fInstanceNumber(fgCounter++), + fMultSoft(kTRUE), + fECMethod(kDefault), + fQTransport(1.), + fMu(1.), + fK(4.e5), + fLengthMax(20), + fLengthMaxOld(0), + fHistos(0), + fHisto(0), + fTablesLoaded(kFALSE) { - fTablesLoaded=kFALSE; - fMultSoft=kTRUE; - fHistos=0; - SetMu(); - SetQTransport(); - SetECMethod(); - SetLengthMax(); - fLengthMaxOld=0; - fInstanceNumber=gCounter++; + //default constructor + } AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a) - : TObject() + : 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=gCounter++; + fInstanceNumber=fgCounter++; + Char_t name[100]; + snprintf(name,100, "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 } @@ -81,12 +119,24 @@ AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a) AliQuenchingWeights::~AliQuenchingWeights() { Reset(); + delete fHisto; +} + +void AliQuenchingWeights::Init() +{ +// Initialization + if (fHisto) return; + fHisto = new TH1F(Form("hhistoqw_%d",fInstanceNumber),"",fgkBins,0.,fgkMaxBin); + for(Int_t bin=1;bin<=fgkBins;bin++) + fHisto->SetBinContent(bin,0.); } 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 + 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) @@ -113,22 +166,23 @@ Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall) fMultSoft=kTRUE; Char_t fname[1024]; - sprintf(fname,"%s",gSystem->ExpandPathName(contall)); - ifstream fincont(fname); + snprintf(fname,1024, "%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[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; @@ -137,33 +191,35 @@ Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall) 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[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)); - ifstream findisc(fname); + snprintf(fname,1024, "%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; + if(nn==34) break; } nn=0; //gluons while(findisc>>frrrg[nn]>>fdag[nn]) { nn++; - if(nn==30) break; + if(nn==34) break; } findisc.close(); fTablesLoaded = kTRUE; @@ -200,18 +256,18 @@ C using swqmult. C C Please, send us any comment: C -C urs.wiedemann@cern.ch -C carlos.salgado@cern.ch -C -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(30), caq(30,261), rrr(30) + REAL*8 xx(400), daq(34), caq(34,261), rrr(34) COMMON /dataqua/ xx, daq, caq, rrr * - REAL*8 xxg(400), dag(30), cag(30,261), rrrg(30) + REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34) COMMON /dataglu/ xxg, dag, cag, rrrg REAL*8 rrrr,xxxx, continuous, discrete @@ -221,15 +277,14 @@ C------------------------------------------------------------------- REAL*8 xfraclow, xfrachigh REAL*8 clow, chigh * + + continuous=0.d0 + discrete=0.d0 + rrin = rrrr xxin = xxxx * - nxlow = int(xxin/0.005) + 1 - nxhigh = nxlow + 1 - xfraclow = (xx(nxhigh)-xxin)/0.005 - xfrachigh = (xxin - xx(nxlow))/0.005 -* - do 666, nr=1,30 + do 666, nr=1,34 if (rrin.lt.rrr(nr)) then rrhigh = rrr(nr) else @@ -244,6 +299,31 @@ C------------------------------------------------------------------- * 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) @@ -255,6 +335,8 @@ C------------------------------------------------------------------- continuous = rfraclow*clow + rfrachigh*chigh +245 continue + if(ipart.eq.1) then discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh) else @@ -264,13 +346,13 @@ C------------------------------------------------------------------- END subroutine initmult - REAL*8 xxq(400), daq(30), caq(30,261), rrr(30) + REAL*8 xxq(400), daq(34), caq(34,261), rrr(34) COMMON /dataqua/ xxq, daq, caq, rrr * - REAL*8 xxg(400), dag(30), cag(30,261), rrrg(30) + REAL*8 xxg(400), dag(34), cag(34,261), rrrg(34) COMMON /dataglu/ xxg, dag, cag, rrrg * - OPEN(UNIT=20,FILE='cont.all',STATUS='OLD',ERR=90) + 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), @@ -282,7 +364,8 @@ C------------------------------------------------------------------- + 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(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), @@ -295,15 +378,16 @@ C------------------------------------------------------------------- + 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(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='disc.all',STATUS='OLD',ERR=91) - do 112 nn=1,30 + 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,30 + do 113 nn=1,34 read (21,*) rrrg(nn), dag(nn) 113 continue close(21) @@ -315,10 +399,12 @@ C------------------------------------------------------------------- end + ======================================================================= Adapted to ROOT macro by A. Dainese - 13/07/2003 - ported to class by C. Loizides - 12/02/2004 + 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, @@ -328,7 +414,11 @@ Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx, // weights for given parton type, // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2 - // read-in data before first call + //set result to zero + continuous=0.; + discrete=0.; + + //read-in data before first call if(!fTablesLoaded){ Error("CalcMult","Tables are not loaded."); return -1; @@ -341,18 +431,19 @@ Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx, 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[29]) rrin = 1.05*frrr[29]; // AD + 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<=30; nr++) { + 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; @@ -382,7 +492,7 @@ Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx, continuous = rfraclow*clow + rfrachigh*chigh; //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n", - // rfraclow,clow,rfrachigh,chigh,continuous); + //rfraclow,clow,rfrachigh,chigh,continuous); if(ipart==1) { discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1]; @@ -402,9 +512,14 @@ Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *di fMultSoft=kFALSE; Char_t fname[1024]; - sprintf(fname,"%s",gSystem->ExpandPathName(contall)); - ifstream fincont(fname); + snprintf(fname, 1024, "%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]>> @@ -440,9 +555,14 @@ Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *di } fincont.close(); - sprintf(fname,"%s",gSystem->ExpandPathName(discall)); - ifstream findisc(fname); + snprintf(fname, 1024, "%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]) { @@ -621,11 +741,11 @@ Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xx // read-in data before first call if(!fTablesLoaded){ - Error("CalcMult","Tables are not loaded."); + Error("CalcSingleHard","Tables are not loaded."); return -1; } - if(!fMultSoft){ - Error("CalcMult","Tables are not loaded for Single Hard Scattering."); + if(fMultSoft){ + Error("CalcSingleHard","Tables are not loaded for Single Hard Scattering."); return -1; } @@ -688,6 +808,10 @@ Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t w,Double_t qtransp,Double_t length, 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 + Double_t wc=CalcWC(qtransp,length); Double_t rrrr=CalcR(wc,length); Double_t xxxx=w/wc; @@ -698,6 +822,10 @@ Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t w,Double_t mu,Double_t length, 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 + Double_t wcbar=CalcWCbar(mu,length); Double_t rrrr=CalcR(wcbar,length); Double_t xxxx=w/wcbar; @@ -706,12 +834,30 @@ Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const { - Double_t R = wc*l*gkConvFmToInvGeV; - if(R>gkRMax) { - Warning("CalcR","Value of R = %.2f; should be less than %.2f",R,gkRMax); - return -R; + //calculate r value and + //check if it is less then maximum + + Double_t r = wc*l*fgkConvFmToInvGeV; + if(r >= 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; + return r; } Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const @@ -727,34 +873,29 @@ Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Doubl } if((ipart<1) || (ipart>2)) { Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart); - return -1000; + return -1000.; } - Int_t l=(Int_t)length; - if((length-(Double_t)l)>0.5) l++; + Int_t l=GetIndex(length); 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.); + Double_t ret = 2.*e; + Int_t ws=0; + while(ret>e){ + ret=fHistos[ipart-1][l-1]->GetRandom(); + if(++ws==1e6){ + Warning("GetELossRandom", + "Stopping reweighting; maximum loss assigned after 1e6 trials."); + return e; + } } - 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; } + //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 @@ -768,6 +909,11 @@ Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, D 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.; @@ -785,6 +931,265 @@ Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double 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", + "Stopping 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 + } + if (!fHisto) Init(); + + 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.; + if (!fHisto) Init(); + 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 @@ -800,14 +1205,14 @@ Int_t AliQuenchingWeights::SampleEnergyLoss() // read-in data before first call if(!fTablesLoaded){ - Error("CalcMult","Tables are not loaded."); + 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,gkRMax); + Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax); fLengthMax=lmax; } } else { @@ -816,8 +1221,8 @@ Int_t AliQuenchingWeights::SampleEnergyLoss() Reset(); fHistos=new TH1F**[2]; - fHistos[0]=new TH1F*[fLengthMax]; - fHistos[1]=new TH1F*[fLengthMax]; + fHistos[0]=new TH1F*[4*fLengthMax]; + fHistos[1]=new TH1F*[4*fLengthMax]; fLengthMaxOld=fLengthMax; //remember old value in case //user wants to reset @@ -825,49 +1230,105 @@ Int_t AliQuenchingWeights::SampleEnergyLoss() Char_t meddesc[100]; if(fMultSoft) { medvalue=(Int_t)(fQTransport*1000.); - sprintf(meddesc,"MS"); + snprintf(meddesc, 100, "MS"); } else { medvalue=(Int_t)(fMu*1000.); - sprintf(meddesc,"SH"); + snprintf(meddesc, 100, "SH"); } for(Int_t ipart=1;ipart<=2;ipart++){ - for(Int_t l=1;l<=fLengthMax;l++){ - + 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 wc = CalcWC(l); - fHistos[ipart-1][l-1] = new TH1F(hname,hname,1100,0.,1.1*wc); + snprintf(hname, 100, "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,l); + Double_t rrrr = CalcR(wc,len); Double_t discrete=0.; // loop on histogram channels - for(Int_t bin=1; bin<=1100; bin++) { + for(Int_t bin=1; bin<=fgkBins; bin++) { Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc; Double_t continuous; - CalcMult(ipart,rrrr,xxxx,continuous,discrete); + 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<=1100;bin++) + 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,1100); + 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,1100); + Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins); fHistos[ipart-1][l-1]->Scale(1./hint); } } return 0; } -const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Int_t l) const +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; + if (!fHisto) Init(); + 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; @@ -877,9 +1338,8 @@ const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Int_t l) const return 0; } + Int_t l=GetIndex(length); if(l<=0) return 0; - if(l>fLengthMax) l=fLengthMax; - return fHistos[ipart-1][l-1]; } @@ -890,34 +1350,30 @@ TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t l // 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) + // - 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"); + snprintf(meddesc, 100, "MS"); } else { wc=CalcWCbar(medval,length); - sprintf(meddesc,"SH"); + snprintf(meddesc, 100, "SH"); } Char_t hname[100]; - sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart, + snprintf(hname, 100, "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); + 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<=1100; bin++) { + //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; @@ -939,27 +1395,23 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t // 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) + // - 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"); + snprintf(meddesc, 100, "MS"); } else { wc=CalcWCbar(medval,length); - sprintf(meddesc,"SH"); + snprintf(meddesc, 100, "SH"); } Char_t hname[100]; - sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart, + snprintf(hname, 100, "hContQWHistox_%s_%d_%d_%d",meddesc,ipart, (Int_t)(medval*1000.),(Int_t)length); - TH1F *histx = new TH1F("histx",hname,1100,0.,1.1); + 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})"); @@ -968,8 +1420,8 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t histx->SetLineColor(4); Double_t rrrr = CalcR(wc,length); - // loop on histogram channels - for(Int_t bin=1; bin<=1100; bin++) { + //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; @@ -984,8 +1436,63 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t 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) { + snprintf(meddesc, 100, "MS"); + } else { + snprintf(meddesc, 100, "SH"); + } + + Char_t hname[100]; + snprintf(hname, 100, "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); @@ -999,11 +1506,11 @@ TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_ Char_t name[100]; Char_t hname[100]; if(ipart==1){ - sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); - sprintf(hname,"hLossQuarks"); + snprintf(name, 100, "Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); + snprintf(hname,100, "hLossQuarks"); } else { - sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); - sprintf(hname,"hLossGluons"); + snprintf(name, 100, "Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); + snprintf(hname, 100, "hLossGluons"); } TH1F *h = new TH1F(hname,name,250,0,250); @@ -1012,15 +1519,17 @@ TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_ 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); @@ -1034,11 +1543,11 @@ TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *h Char_t name[100]; Char_t hname[100]; if(ipart==1){ - sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); - sprintf(hname,"hLossQuarks"); + snprintf(name, 100, "Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); + snprintf(hname,100, "hLossQuarks"); } else { - sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); - sprintf(hname,"hLossGluons"); + snprintf(name,100, "Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); + snprintf(hname, 100, "hLossGluons"); } TH1F *h = new TH1F(hname,name,250,0,250); @@ -1047,15 +1556,72 @@ TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *h 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 +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]; + snprintf(hname, 100, "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); @@ -1063,85 +1629,96 @@ void AliQuenchingWeights::PlotDiscreteWeights(Int_t len) const 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); + 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("Probability #Delta E = 0 , p_{0}"); + hframe->SetYTitle("p_{0} (discrete weight)"); hframe->Draw(); - TGraph *gq=new TGraph(20); + 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<=5.05;q+=0.25){ + for(Double_t q=0.05;q<=qm+.05;q+=0.25){ Double_t disc,cont; - CalcMult(1,1.0, q, len, cont, disc); + 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){ + for(Double_t m=0.05;m<=qm+.05;m+=0.25){ Double_t disc,cont; - CalcSingleHard(1,1.0, m, len, cont, disc); + CalcSingleHard(1,1.0,m,len,cont, disc); gq->SetPoint(i,m,disc);i++; } } gq->SetMarkerStyle(20); - gq->Draw("pl"); + gq->SetMarkerColor(1); + gq->SetLineStyle(1); + gq->SetLineColor(1); + gq->Draw("l"); - TGraph *gg=new TGraph(20); + TGraph *gg=new TGraph(points); i=0; if(fMultSoft){ - for(Double_t q=0.05;q<=5.05;q+=0.25){ + for(Double_t q=0.05;q<=qm+.05;q+=0.25){ Double_t disc,cont; - CalcMult(2,1.0, q, 5., cont, disc); + CalcMult(2,1.0,q,len,cont,disc); gg->SetPoint(i,q,disc);i++; } } else { - for(Double_t m=0.05;m<=5.05;m+=0.25){ + for(Double_t m=0.05;m<=qm+.05;m+=0.25){ Double_t disc,cont; - CalcSingleHard(2,1.0, m, 5., cont, disc); + CalcSingleHard(2,1.0,m,len,cont,disc); gg->SetPoint(i,m,disc);i++; } } gg->SetMarkerStyle(24); - gg->Draw("pl"); + 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 = %d fm",len); + snprintf(label, 100, "L = %.1f fm",len); l1a->AddEntry(gq,label,""); - l1a->AddEntry(gq,"quark","pl"); - l1a->AddEntry(gg,"gluon","pl"); + l1a->AddEntry(gq,"quark projectile","l"); + l1a->AddEntry(gg,"gluon projectile","l"); l1a->Draw(); c->Update(); } -void AliQuenchingWeights::PlotContWeights(Int_t itype,Int_t ell) const +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"); + snprintf(title, 1024, "Cont. Weight for Multiple Scattering - Quarks"); else if(itype==2) - sprintf(title,"Cont. Weight for Multiple Scattering - Gluons"); + snprintf(title, 1024, "Cont. Weight for Multiple Scattering - Gluons"); else return; medvals[0]=4;medvals[1]=1;medvals[2]=0.5; - sprintf(name,"ccont-ms-%d",itype); + snprintf(name, 1024, "ccont-ms-%d",itype); } else { if(itype==1) - sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks"); + snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Quarks"); else if(itype==2) - sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons"); + snprintf(title, 1024, "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); + snprintf(name, 1024, "ccont-ms-%d",itype); } TCanvas *c = new TCanvas(name,title,0,0,500,400); @@ -1165,21 +1742,21 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Int_t ell) const l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; - sprintf(label,"L = %d fm",ell); + snprintf(label, 100, "L = %.1f fm",ell); l1a->AddEntry(h1,label,""); if(fMultSoft) { - sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]); + snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[0]); l1a->AddEntry(h1,label,"pl"); - sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]); + snprintf(label, 100, "#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"); + snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[2]); + l1a->AddEntry(h3, label,"pl"); } else { - sprintf(label,"#mu = %.1f GeV",medvals[0]); + snprintf(label, 100, "#mu = %.1f GeV",medvals[0]); l1a->AddEntry(h1,label,"pl"); - sprintf(label,"#mu = %.1f GeV",medvals[1]); + snprintf(label, 100, "#mu = %.1f GeV",medvals[1]); l1a->AddEntry(h2,label,"pl"); - sprintf(label,"#mu = %.1f GeV",medvals[2]); + snprintf(label, 100, "#mu = %.1f GeV",medvals[2]); l1a->AddEntry(h3,label,"pl"); } l1a->Draw(); @@ -1187,24 +1764,27 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Int_t ell) const c->Update(); } -void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const +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"); + snprintf(title,1024, "Cont. Weight for Multiple Scattering - Quarks"); else if(itype==2) - sprintf(title,"Cont. Weight for Multiple Scattering - Gluons"); + snprintf(title,1024, "Cont. Weight for Multiple Scattering - Gluons"); else return; - sprintf(name,"ccont2-ms-%d",itype); + snprintf(name,1024, "ccont2-ms-%d",itype); } else { if(itype==1) - sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks"); + snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Quarks"); else if(itype==2) - sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons"); + snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Gluons"); else return; - sprintf(name,"ccont2-sh-%d",itype); + snprintf(name, 1024, "ccont2-sh-%d",itype); } TCanvas *c = new TCanvas(name,title,0,0,500,400); c->cd(); @@ -1228,9 +1808,9 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const l1a->SetBorderSize(0); Char_t label[100]; if(fMultSoft) - sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval); + snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medval); else - sprintf(label,"#mu = %.1f GeV",medval); + snprintf(label, 100, "#mu = %.1f GeV",medval); l1a->AddEntry(h1,label,""); l1a->AddEntry(h1,"L = 8 fm","pl"); @@ -1241,26 +1821,29 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const c->Update(); } -void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e) const +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("CalcMult","Tables are not loaded."); + Error("PlotAvgELoss","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"); + snprintf(title, 1024, "Average Energy Loss for Multiple Scattering"); + snprintf(name, 1024, "cavgelossms"); } else { - sprintf(title,"Average Loss for Single Hard Scattering"); - sprintf(name,"cavgelosssh"); + snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering"); + snprintf(name, 1024, "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); + TH2F *hframe = new TH2F("avgloss","",2,0,qm+.1,2,0,100); hframe->SetStats(0); if(fMultSoft) hframe->SetXTitle("#hat{q} [GeV^{2}/fm]"); @@ -1271,34 +1854,48 @@ void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e) const TGraph *gq=new TGraph(20); Int_t i=0; - for(Double_t v=0.05;v<=5.05;v+=0.25){ + 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(20); + gq->SetMarkerStyle(21); gq->Draw("pl"); - TGraph *gg=new TGraph(20); + Int_t points=(Int_t)qm*4; + TGraph *gg=new TGraph(points); i=0; - for(Double_t v=0.05;v<=5.05;v+=0.25){ + 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(24); + gg->SetMarkerStyle(20); + gg->SetMarkerColor(2); gg->Draw("pl"); - TLegend *l1a = new TLegend(0.5,0.6,.95,0.8); + 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 = %d fm",len); + snprintf(label, 100, "L = %.1f fm",len); l1a->AddEntry(gq,label,""); - l1a->AddEntry(gq,"quark","pl"); - l1a->AddEntry(gg,"gluon","pl"); + 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(); @@ -1306,19 +1903,22 @@ void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e) const 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."); + Error("PlotAvgELossVs","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"); + snprintf(title, 1024, "Average Energy Loss for Multiple Scattering"); + snprintf(name, 1024, "cavgelossms2"); } else { - sprintf(title,"Average Loss for Single Hard Scattering"); - sprintf(name,"cavgelosssh2"); + snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering"); + snprintf(name, 1024, "cavgelosssh2"); } TCanvas *c = new TCanvas(name,title,0,0,500,400); @@ -1354,34 +1954,135 @@ void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const 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()); + snprintf(label, 100, " = %.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::PlotAvgELossVsPt(Double_t medval,Int_t len) const +void AliQuenchingWeights::PlotAvgELossVsL(Double_t e) const { + // plot average energy loss versus ell + if(!fTablesLoaded){ - Error("CalcMult","Tables are not loaded."); + Error("PlotAvgELossVsEll","Tables are not loaded."); + return; + } + + Char_t title[1024]; + Char_t name[1024]; + Float_t medval; + if(fMultSoft){ + snprintf(title, 1024, "Average Energy Loss for Multiple Scattering"); + snprintf(name, 1024, "cavgelosslms"); + medval=fQTransport; + } else { + snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering"); + snprintf(name, 1024, "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) + snprintf(label, 100, "#hat{q} = %.2f GeV^{2}/fm",medval); + else + snprintf(label, 100, "#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 Loss for Multiple Scattering"); - sprintf(name,"cavgelossvsptms"); + snprintf(title, 1024, "Relative Energy Loss for Multiple Scattering"); + snprintf(name, 1024, "cavgelossvsptms"); } else { - sprintf(title,"Relative Loss for Single Hard Scattering"); - sprintf(name,"cavgelossvsptsh"); + snprintf(title, 1024, "Relative Energy Loss for Single Hard Scattering"); + snprintf(name, 1024, "cavgelossvsptsh"); } TCanvas *c = new TCanvas(name,title,0,0,500,400); @@ -1418,7 +2119,7 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Int_t len) const l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; - sprintf(label,"L = %d fm",len); + snprintf(label, 100, "L = %.1f fm",len); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark","pl"); l1a->AddEntry(gg,"gluon","pl"); @@ -1429,19 +2130,22 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Int_t len) const 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."); + Error("PlotAvgELossVsPt","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"); + snprintf(title, 1024, "Relative Energy Loss for Multiple Scattering"); + snprintf(name, 1024, "cavgelossvsptms2"); } else { - sprintf(title,"Relative Loss for Single Hard Scattering"); - sprintf(name,"cavgelossvsptsh2"); + snprintf(title, 1024, "Relative Energy Loss for Single Hard Scattering"); + snprintf(name, 1024, "cavgelossvsptsh2"); } TCanvas *c = new TCanvas(name,title,0,0,500,400); c->cd(); @@ -1477,7 +2181,7 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; - sprintf(label," = %.2f fm",hEll->GetMean()); + snprintf(label, 100, " = %.2f fm",hEll->GetMean()); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark","pl"); l1a->AddEntry(gg,"gluon","pl"); @@ -1486,3 +2190,13 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const 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; +} +