/* $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 <Riostream.h>
+#include <TF1.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TROOT.h>
#include <TSystem.h>
+#include <TString.h>
#include <TLegend.h>
#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
}
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<fLengthMaxOld;l++){
+ for(Int_t l=0;l<4*fLengthMaxOld;l++){
delete fHistos[0][l];
delete fHistos[1][l];
}
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
+ 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)
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;
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;
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
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
*
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)
continuous = rfraclow*clow + rfrachigh*chigh
+245 continue
+
if(ipart.eq.1) then
discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh)
else
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),
+ 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),
+ 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)
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,
// 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;
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(rrin<frrr[nr-1]) {
rrhigh = frrr[nr-1];
} else {
Double_t rfraclow = (rrhigh-rrin)/(rrhigh-rrlow);
Double_t rfrachigh = (rrin-rrlow)/(rrhigh-rrlow);
+ if(rrin>1.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;
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];
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]>>
}
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]) {
// 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;
}
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;
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;
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
}
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
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.;
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
// 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 {
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
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;
return 0;
}
+ Int_t l=GetIndex(length);
if(l<=0) return 0;
- if(l>fLengthMax) l=fLengthMax;
-
return fHistos[ipart-1][l-1];
}
// 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_<ipart>_<medval*1000>_<l>
- // 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;
// 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_<ipart>_<medval*1000>_<l>
- // 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})");
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;
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);
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);
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);
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);
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);
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);
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();
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();
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");
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]");
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;i<points;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.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();
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);
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,"<L> = %.2f fm",hEll->GetMean());
+ snprintf(label, 100, "<L> = %.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("<E_{loss}> [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);
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");
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();
l1a->SetFillStyle(0);
l1a->SetBorderSize(0);
Char_t label[100];
- sprintf(label,"<L> = %.2f fm",hEll->GetMean());
+ snprintf(label, 100, "<L> = %.2f fm",hEll->GetMean());
l1a->AddEntry(gq,label,"");
l1a->AddEntry(gq,"quark","pl");
l1a->AddEntry(gg,"gluon","pl");
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;
+}
+