AliQuenchingWeight first commit (K. Loizides, A. Dainese)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 27 Feb 2004 15:40:09 +0000 (15:40 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 27 Feb 2004 15:40:09 +0000 (15:40 +0000)
FASTSIM/AliQuenchingWeights.cxx [new file with mode: 0644]
FASTSIM/AliQuenchingWeights.h [new file with mode: 0644]

diff --git a/FASTSIM/AliQuenchingWeights.cxx b/FASTSIM/AliQuenchingWeights.cxx
new file mode 100644 (file)
index 0000000..a846437
--- /dev/null
@@ -0,0 +1,1488 @@
+/**************************************************************************
+ * 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();
+}
+
diff --git a/FASTSIM/AliQuenchingWeights.h b/FASTSIM/AliQuenchingWeights.h
new file mode 100644 (file)
index 0000000..f87d96e
--- /dev/null
@@ -0,0 +1,131 @@
+#ifndef ALIQUENCHINGWEIGHTS_H
+#define ALIQUENCHINGWEIGHTS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $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   constantin.loizides@cern.ch
+//                     A. Dainese    andrea.dainese@pd.infn.it            
+//----------------------------------------------------------------------------
+
+#include <TObject.h>
+class TH1F;
+
+class AliQuenchingWeights : public TObject {
+ public:
+  enum kECMethod {kDefault=0,kReweight=1};
+
+  AliQuenchingWeights();
+  AliQuenchingWeights(const AliQuenchingWeights& a);
+  virtual ~AliQuenchingWeights();
+
+  void Reset();
+  Int_t SampleEnergyLoss();
+  Double_t GetELossRandom(Int_t ipart, Double_t length, Double_t e=1.e6) const;
+  Double_t CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e)  const;
+  Double_t GetELossRandom(Int_t ipart, TH1F *hell, Double_t e=1.e6) const;
+  Double_t CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double_t e)  const;
+
+  //multiple soft scattering approximation
+  Int_t InitMult(const Char_t *contall="$(ALICE_ROOT)/FASTSIM/data/cont_mult.all",
+                 const Char_t *discall="$(ALICE_ROOT)/FASTSIM/data/disc_mult.all"); 
+
+  //single hard scattering approximation
+  Int_t InitSingleHard(const Char_t *contall="$(ALICE_ROOT)/FASTSIM/data/cont_lin.all",
+                       const Char_t *discall="$(ALICE_ROOT)/FASTSIM/data/disc_lin.all"); 
+
+  Int_t CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
+                 Double_t &continuous,Double_t &discrete) const;
+  Int_t CalcMult(Int_t ipart, 
+                Double_t w, Double_t qtransp, Double_t length,
+                 Double_t &continuous,Double_t &discrete) const;
+  Int_t CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xxxx,
+                      Double_t &continuous,Double_t &discrete) const;
+  Int_t CalcSingleHard(Int_t ipart, 
+                      Double_t w, Double_t mu, Double_t length,
+                       Double_t &continuous,Double_t &discrete) const;
+
+  Double_t CalcWC(Double_t q, Double_t l) const 
+    {return 0.5*q*l*l*gkConvFmToInvGeV;}
+
+  Double_t CalcWCbar(Double_t mu, Double_t l) const 
+    {return 0.5*mu*mu*l*gkConvFmToInvGeV;}
+
+  Double_t CalcWC(Double_t l) const 
+    {if(fMultSoft) return CalcWC(fQTransport,l);
+     else return CalcWCbar(fMu,l);}
+
+  Double_t CalcR(Double_t wc, Double_t l) const; 
+
+  Int_t CalcLengthMax(Double_t q) const
+    {Double_t l3max=gkRMax/.5/q/gkConvFmToInvGeV/gkConvFmToInvGeV;
+     return (Int_t)TMath::Power(l3max,1./3.);} 
+
+  const TH1F* GetHisto(Int_t ipart,Int_t l) const;
+
+  void SetMu(Double_t m=1.) {fMu=m;}
+  void SetQTransport(Double_t q=1.) {fQTransport=q;}
+  void SetECMethod(kECMethod type=kDefault);
+  void SetLengthMax(Int_t l=20) {fLengthMax=l;}
+
+  Float_t GetMu()           const {return fMu;}
+  Float_t GetQTransport()   const {return fQTransport;}
+  Bool_t  GetECMethod()     const {return fECMethod;}
+  Bool_t  GetTablesLoaded() const {return fTablesLoaded;}
+  Bool_t  GetMultSoft()     const {return fMultSoft;}
+  Int_t   GetLengthMax()    const {return fLengthMax;}
+
+  TH1F* ComputeQWHisto (Int_t ipart,Double_t medval,Double_t length)  const; 
+  TH1F* ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t length)  const; 
+  TH1F* ComputeELossHisto (Int_t ipart,Double_t medval,Double_t l,Double_t e=1.e6) const; 
+  TH1F* ComputeELossHisto (Int_t ipart,Double_t medval,TH1F *hEll,Double_t e=1.e6) const; 
+  
+  void PlotDiscreteWeights(Int_t len=4) const;
+  void PlotContWeights(Int_t itype,Int_t len) const;
+  void PlotContWeights(Int_t itype,Double_t medval) const;
+  void PlotAvgELoss(Int_t len ,Double_t e=1.e6) const;
+  void PlotAvgELoss(TH1F *hEll,Double_t e=1.e6) const;
+  void PlotAvgELossVsPt(Double_t medval,Int_t len) const;
+  void PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const;
+
+ protected:
+  static const Double_t gkConvFmToInvGeV; 
+  static const Double_t gkRMax; 
+  static Int_t gCounter;//static instance counter
+  Int_t fInstanceNumber;//instance number of class
+
+  Bool_t fMultSoft;     //approximation type
+  Bool_t fECMethod;     //energy constraint method
+  Double_t fQTransport; //transport coefficient
+  Double_t fMu;         //Debye screening mass
+  Int_t fLengthMax;     //maximum length
+  Int_t fLengthMaxOld;  //maximum length used for histos
+
+  //discrete and cont part of quenching for
+  //both parton type and different lengths
+  TH1F ***fHistos; //!
+
+  // data strucs for tables
+  Double_t fxx[400];
+  Double_t fxxg[400];
+  Double_t fdaq[30];
+  Double_t fdag[30];
+  Double_t fcaq[30][261];
+  Double_t fcag[30][261];  
+  Double_t frrr[30];
+  Double_t frrrg[30];
+  Bool_t fTablesLoaded;
+
+  ClassDef(AliQuenchingWeights,1)    // Base class for Quenching Weights
+};
+
+#endif