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