Do not include from subdirectories
[u/mrichter/AliRoot.git] / FASTSIM / AliQuenchingWeights.cxx
index a846437..c3c1057 100644 (file)
 
 /* $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
 }
@@ -81,12 +119,24 @@ AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a)
 AliQuenchingWeights::~AliQuenchingWeights()
 {
   Reset();
+  delete fHisto;
+}
+
+void AliQuenchingWeights::Init()
+{
+//    Initialization
+    if (fHisto) return;
+    fHisto = new TH1F(Form("hhistoqw_%d",fInstanceNumber),"",fgkBins,0.,fgkMaxBin);
+    for(Int_t bin=1;bin<=fgkBins;bin++) 
+       fHisto->SetBinContent(bin,0.);
 }
 
 void AliQuenchingWeights::Reset()
 {
+  //reset tables if there were used
+
   if(!fHistos) return;
-  for(Int_t l=0;l<fLengthMaxOld;l++){
+  for(Int_t l=0;l<4*fLengthMaxOld;l++){
     delete fHistos[0][l];
     delete fHistos[1][l];
   }
@@ -97,11 +147,14 @@ void AliQuenchingWeights::Reset()
 
 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) 
@@ -113,22 +166,23 @@ Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
   fMultSoft=kTRUE;
   
   Char_t fname[1024];
-  sprintf(fname,"%s",gSystem->ExpandPathName(contall));
-  ifstream fincont(fname);
+  snprintf(fname,1024, "%s",gSystem->ExpandPathName(contall));
+  //PH  ifstream fincont(fname);
+  fstream fincont(fname,ios::in);
+#if defined(__HP_aCC) || defined(__DECCXX)
+  if(!fincont.rdbuf()->is_open()) return -1;
+#else
   if(!fincont.is_open()) return -1;
+#endif
 
   Int_t nn=0; //quarks
   while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
        fcaq[4][nn]>>fcaq[5][nn]>>fcaq[6][nn]>>fcaq[7][nn]>>fcaq[8][nn]>>
-       fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>
-       fcaq[13][nn]>>
-       fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>
-       fcaq[18][nn]>>
-       fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>
-       fcaq[23][nn]>>
-       fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>
-       fcaq[28][nn]>>
-       fcaq[29][nn]) 
+       fcaq[9][nn]>>fcaq[10][nn]>>fcaq[11][nn]>>fcaq[12][nn]>>fcaq[13][nn]>>
+       fcaq[14][nn]>>fcaq[15][nn]>>fcaq[16][nn]>>fcaq[17][nn]>>fcaq[18][nn]>>
+       fcaq[19][nn]>>fcaq[20][nn]>>fcaq[21][nn]>>fcaq[22][nn]>>fcaq[23][nn]>>
+       fcaq[24][nn]>>fcaq[25][nn]>>fcaq[26][nn]>>fcaq[27][nn]>>fcaq[28][nn]>>
+       fcaq[29][nn]>>fcaq[30][nn]>>fcaq[31][nn]>>fcaq[32][nn]>>fcaq[33][nn]) 
     {
       nn++;
       if(nn==261) break;
@@ -137,33 +191,35 @@ Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
   nn=0;       //gluons
   while(fincont>>fxxg[nn]>>fcag[0][nn]>>fcag[1][nn]>>fcag[2][nn]>>fcag[3][nn]>>
        fcag[4][nn]>>fcag[5][nn]>>fcag[6][nn]>>fcag[7][nn]>>fcag[8][nn]>>
-       fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>
-       fcag[13][nn]>>
-       fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>
-       fcag[18][nn]>>
-       fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>
-       fcag[23][nn]>>
-       fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>
-       fcag[28][nn]>>
-       fcag[29][nn]) {
+       fcag[9][nn]>>fcag[10][nn]>>fcag[11][nn]>>fcag[12][nn]>>fcag[13][nn]>>
+        fcag[14][nn]>>fcag[15][nn]>>fcag[16][nn]>>fcag[17][nn]>>fcag[18][nn]>>
+       fcag[19][nn]>>fcag[20][nn]>>fcag[21][nn]>>fcag[22][nn]>>fcag[23][nn]>>
+       fcag[24][nn]>>fcag[25][nn]>>fcag[26][nn]>>fcag[27][nn]>>fcag[28][nn]>>
+       fcag[29][nn]>>fcag[30][nn]>>fcag[31][nn]>>fcag[32][nn]>>fcag[33][nn]) 
+  {
     nn++;
     if(nn==261) break;
   }
   fincont.close();
 
-  sprintf(fname,"%s",gSystem->ExpandPathName(discall));
-  ifstream findisc(fname); 
+  snprintf(fname,1024, "%s",gSystem->ExpandPathName(discall));
+  //PH  ifstream findisc(fname); 
+  fstream findisc(fname,ios::in); 
+#if defined(__HP_aCC) || defined(__DECCXX)
+  if(!findisc.rdbuf()->is_open()) return -1;
+#else
   if(!findisc.is_open()) return -1;
+#endif
 
   nn=0; //quarks
   while(findisc>>frrr[nn]>>fdaq[nn]) {
     nn++;
-    if(nn==30) break;
+    if(nn==34) break;
   }
   nn=0; //gluons
   while(findisc>>frrrg[nn]>>fdag[nn]) {
     nn++;
-    if(nn==30) break;
+    if(nn==34) break;
   }
   findisc.close();
   fTablesLoaded = kTRUE;
@@ -200,18 +256,18 @@ C   using swqmult.
 C
 C   Please, send us any comment:
 C
-C      urs.wiedemann@cern.ch 
-C      carlos.salgado@cern.ch 
-C 
-C 
+C      urs.wiedemann@cern.ch
+C      carlos.salgado@cern.ch
+C
+C
 C-------------------------------------------------------------------
 
       SUBROUTINE swqmult(ipart,rrrr,xxxx,continuous,discrete)
 *
-      REAL*8           xx(400), daq(30), caq(30,261), rrr(30)
+      REAL*8           xx(400), daq(34), caq(34,261), rrr(34)
       COMMON /dataqua/    xx, daq, caq, rrr
 *
-      REAL*8           xxg(400), dag(30), cag(30,261), rrrg(30)
+      REAL*8           xxg(400), dag(34), cag(34,261), rrrg(34)
       COMMON /dataglu/    xxg, dag, cag, rrrg
 
       REAL*8           rrrr,xxxx, continuous, discrete
@@ -221,15 +277,14 @@ C-------------------------------------------------------------------
       REAL*8           xfraclow, xfrachigh
       REAL*8           clow, chigh
 *
+
+      continuous=0.d0
+      discrete=0.d0
+
       rrin = rrrr
       xxin = xxxx
 *
-      nxlow = int(xxin/0.005) + 1
-      nxhigh = nxlow + 1
-      xfraclow = (xx(nxhigh)-xxin)/0.005
-      xfrachigh = (xxin - xx(nxlow))/0.005
-*
-      do 666, nr=1,30
+      do 666, nr=1,34
          if (rrin.lt.rrr(nr)) then
             rrhigh = rrr(nr)
          else
@@ -244,6 +299,31 @@ C-------------------------------------------------------------------
 *
       rfraclow = (rrhigh-rrin)/(rrhigh-rrlow)
       rfrachigh = (rrin-rrlow)/(rrhigh-rrlow)
+      if (rrin.gt.10000d0) then
+         rfraclow = dlog(rrhigh/rrin)/dlog(rrhigh/rrlow)
+         rfrachigh = dlog(rrin/rrlow)/dlog(rrhigh/rrlow)
+      endif
+*
+      if (ipart.eq.1.and.rrin.ge.rrr(1)) then
+         nrlow=1
+         nrhigh=1
+         rfraclow=1
+         rfrachigh=0
+      endif
+
+      if (ipart.ne.1.and.rrin.ge.rrrg(1)) then
+         nrlow=1
+         nrhigh=1
+         rfraclow=1
+         rfrachigh=0
+      endif
+
+      if (xxxx.ge.xx(261)) go to 245
+
+      nxlow = int(xxin/0.01) + 1
+      nxhigh = nxlow + 1
+      xfraclow = (xx(nxhigh)-xxin)/0.01
+      xfrachigh = (xxin - xx(nxlow))/0.01
 *
       if(ipart.eq.1) then
       clow = xfraclow*caq(nrlow,nxlow)+xfrachigh*caq(nrlow,nxhigh)
@@ -255,6 +335,8 @@ C-------------------------------------------------------------------
 
       continuous = rfraclow*clow + rfrachigh*chigh
 
+245   continue
+
       if(ipart.eq.1) then
       discrete = rfraclow*daq(nrlow) + rfrachigh*daq(nrhigh)
       else
@@ -264,13 +346,13 @@ C-------------------------------------------------------------------
       END
 
       subroutine initmult
-      REAL*8           xxq(400), daq(30), caq(30,261), rrr(30)
+      REAL*8           xxq(400), daq(34), caq(34,261), rrr(34)
       COMMON /dataqua/    xxq, daq, caq, rrr
 *
-      REAL*8           xxg(400), dag(30), cag(30,261), rrrg(30)
+      REAL*8           xxg(400), dag(34), cag(34,261), rrrg(34)
       COMMON /dataglu/    xxg, dag, cag, rrrg
 *
-      OPEN(UNIT=20,FILE='cont.all',STATUS='OLD',ERR=90)
+      OPEN(UNIT=20,FILE='contnew.all',STATUS='OLD',ERR=90)
       do 110 nn=1,261
       read (20,*) xxq(nn), caq(1,nn), caq(2,nn), caq(3,nn),
      +     caq(4,nn), caq(5,nn), caq(6,nn), caq(7,nn), caq(8,nn),
@@ -282,7 +364,8 @@ C-------------------------------------------------------------------
      +     caq(23,nn),
      +     caq(24,nn), caq(25,nn), caq(26,nn), caq(27,nn), 
      +     caq(28,nn),
-     +     caq(29,nn), caq(30,nn)
+     +     caq(29,nn), caq(30,nn), caq(31,nn), caq(32,nn), 
+     +     caq(33,nn), caq(34,nn)
  110     continue
       do 111 nn=1,261
       read (20,*) xxg(nn), cag(1,nn), cag(2,nn), cag(3,nn),
@@ -295,15 +378,16 @@ C-------------------------------------------------------------------
      +     cag(23,nn),
      +     cag(24,nn), cag(25,nn), cag(26,nn), cag(27,nn), 
      +     cag(28,nn),
-     +     cag(29,nn), cag(30,nn)
+     +     cag(29,nn), cag(30,nn), cag(31,nn), cag(32,nn), 
+     +     cag(33,nn), cag(34,nn)
  111     continue
       close(20)
 *
-      OPEN(UNIT=21,FILE='disc.all',STATUS='OLD',ERR=91)
-      do 112 nn=1,30
+      OPEN(UNIT=21,FILE='discnew.all',STATUS='OLD',ERR=91)
+      do 112 nn=1,34
       read (21,*) rrr(nn), daq(nn)
  112     continue
-      do 113 nn=1,30
+      do 113 nn=1,34
       read (21,*) rrrg(nn), dag(nn)
  113     continue
       close(21)
@@ -315,10 +399,12 @@ C-------------------------------------------------------------------
 
       end
 
+
 =======================================================================
 
    Adapted to ROOT macro by A. Dainese - 13/07/2003
-   ported to class by C. Loizides - 12/02/2004
+   Ported to class by C. Loizides - 12/02/2004
+   New version for extended R values added - 06/03/2004 
 */
 
 Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
@@ -328,7 +414,11 @@ Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
   // weights for given parton type, 
   // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
 
-  // read-in data before first call
+  //set result to zero
+  continuous=0.;
+  discrete=0.;
+
+  //read-in data before first call
   if(!fTablesLoaded){
     Error("CalcMult","Tables are not loaded.");
     return -1;
@@ -341,18 +431,19 @@ Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
   Double_t rrin = rrrr;
   Double_t xxin = xxxx;
 
+  if(xxin>fxx[260]) return -1;
   Int_t nxlow     = (Int_t)(xxin/0.01) + 1;
   Int_t nxhigh    = nxlow + 1;
   Double_t xfraclow  = (fxx[nxhigh-1]-xxin)/0.01;
   Double_t xfrachigh = (xxin - fxx[nxlow-1])/0.01;
 
   //why this?
-  if(rrin<=frrr[29]) rrin = 1.05*frrr[29]; // AD
+  if(rrin<=frrr[33]) rrin = 1.05*frrr[33]; // AD
   if(rrin>=frrr[0])  rrin = 0.95*frrr[0];  // AD
 
   Int_t nrlow=0,nrhigh=0;
   Double_t rrhigh=0,rrlow=0;
-  for(Int_t nr=1; nr<=30; nr++) {
+  for(Int_t nr=1; nr<=34; nr++) {
     if(rrin<frrr[nr-1]) {
       rrhigh = frrr[nr-1];
     } else {
@@ -369,6 +460,25 @@ Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
   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;
@@ -382,7 +492,7 @@ Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
 
   continuous = rfraclow*clow + rfrachigh*chigh;
   //printf("rfraclow %f, clow %f, rfrachigh %f, chigh %f,\n continuous %f\n",
-  //    rfraclow,clow,rfrachigh,chigh,continuous);
+  //rfraclow,clow,rfrachigh,chigh,continuous);
 
   if(ipart==1) {
     discrete = rfraclow*fdaq[nrlow-1] + rfrachigh*fdaq[nrhigh-1];
@@ -402,9 +512,14 @@ Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *di
   fMultSoft=kFALSE;
   
   Char_t fname[1024];
-  sprintf(fname,"%s",gSystem->ExpandPathName(contall));
-  ifstream fincont(fname);
+  snprintf(fname, 1024, "%s",gSystem->ExpandPathName(contall));
+  //PH  ifstream fincont(fname);
+  fstream fincont(fname,ios::in);
+#if defined(__HP_aCC) || defined(__DECCXX)
+  if(!fincont.rdbuf()->is_open()) return -1;
+#else
   if(!fincont.is_open()) return -1;
+#endif
 
   Int_t nn=0; //quarks
   while(fincont>>fxx[nn]>>fcaq[0][nn]>>fcaq[1][nn]>>fcaq[2][nn]>>fcaq[3][nn]>>
@@ -440,9 +555,14 @@ Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *di
   }
   fincont.close();
 
-  sprintf(fname,"%s",gSystem->ExpandPathName(discall));
-  ifstream findisc(fname); 
+  snprintf(fname, 1024, "%s",gSystem->ExpandPathName(discall));
+  //PH  ifstream findisc(fname); 
+  fstream findisc(fname,ios::in); 
+#if defined(__HP_aCC) || defined(__DECCXX)
+  if(!findisc.rdbuf()->is_open()) return -1;
+#else
   if(!findisc.is_open()) return -1;
+#endif
 
   nn=0; //quarks
   while(findisc>>frrr[nn]>>fdaq[nn]) {
@@ -621,11 +741,11 @@ Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xx
 
   // read-in data before first call
   if(!fTablesLoaded){
-    Error("CalcMult","Tables are not loaded.");
+    Error("CalcSingleHard","Tables are not loaded.");
     return -1;
   }
-  if(!fMultSoft){
-    Error("CalcMult","Tables are not loaded for Single Hard Scattering.");
+  if(fMultSoft){
+    Error("CalcSingleHard","Tables are not loaded for Single Hard Scattering.");
     return -1;
   }
 
@@ -688,6 +808,10 @@ Int_t AliQuenchingWeights::CalcMult(Int_t ipart,
                               Double_t w,Double_t qtransp,Double_t length,
                               Double_t &continuous,Double_t &discrete) const
 {
+  // Calculate Multiple Scattering approx. 
+  // weights for given parton type, 
+  // rrrr=0.5*q*L^3 and xxxx=w/wc, wc=0.5*q*L^2
+
   Double_t wc=CalcWC(qtransp,length);
   Double_t rrrr=CalcR(wc,length);
   Double_t xxxx=w/wc;
@@ -698,6 +822,10 @@ Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
                                     Double_t w,Double_t mu,Double_t length,
                                     Double_t &continuous,Double_t &discrete) const
 {
+  // calculate Single Hard approx. 
+  // weights for given parton type, 
+  // rrrr=0.5*mu^2*L^2 and xxxx=w/wc, wc=0.5*mu^2*L
+
   Double_t wcbar=CalcWCbar(mu,length);
   Double_t rrrr=CalcR(wcbar,length);
   Double_t xxxx=w/wcbar;
@@ -706,12 +834,30 @@ Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
 
 Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const 
 { 
-  Double_t R = wc*l*gkConvFmToInvGeV;
-  if(R>gkRMax) {
-    Warning("CalcR","Value of R = %.2f; should be less than %.2f",R,gkRMax);
-    return -R;
+  //calculate r value and 
+  //check if it is less then maximum
+
+  Double_t r = wc*l*fgkConvFmToInvGeV;
+  if(r >= fgkRMax) {
+    Warning("CalcR","Value of r = %.2f; should be less than %.2f", r, fgkRMax);
+    return fgkRMax-1;
+  }  
+  return r;
+}
+
+Double_t AliQuenchingWeights::CalcRk(Double_t k, Double_t I0, Double_t I1) const
+{ 
+  //calculate R value and 
+  //check if it is less then maximum
+
+  Double_t r = fgkRMax-1;
+  if(I0>0)
+    r = 2*I1*I1/I0*k;
+  if(r>=fgkRMax) {
+    Warning("CalcRk","Value of r = %.2f; should be less than %.2f",r,fgkRMax);
+    return fgkRMax-1;
   }  
-  return R;
+  return r;
 }
 
 Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Double_t e) const
@@ -727,34 +873,29 @@ Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Doubl
   }
   if((ipart<1) || (ipart>2)) {
     Fatal("GetELossRandom","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
-    return -1000;
+    return -1000.;
   }
 
-  Int_t l=(Int_t)length;
-  if((length-(Double_t)l)>0.5) l++;
+  Int_t l=GetIndex(length);
   if(l<=0) return 0.;
-  if(l>fLengthMax) l=fLengthMax;
 
   if(fECMethod==kReweight){
-    TH1F *dummy=new TH1F(*fHistos[ipart-1][l-1]);
-    dummy->SetName("dummy");
-    for(Int_t bin=dummy->FindBin(e)+1;bin<=1100;bin++){
-      dummy->SetBinContent(bin,0.);
+    Double_t ret = 2.*e;
+    Int_t ws=0;
+    while(ret>e){
+      ret=fHistos[ipart-1][l-1]->GetRandom(); 
+      if(++ws==1e6){
+       Warning("GetELossRandom",
+                "Stopping reweighting; maximum loss assigned after 1e6 trials.");
+       return e;
+      }
     }
-    dummy->Scale(1./dummy->Integral());
-    Double_t ret=dummy->GetRandom();
-    delete dummy;
-    return ret;
-    //****** !! ALTERNATIVE WAY OF DOING IT !! ***
-    //Double_t ret = 2.*e;
-    //while(ret>e) ret=fHistos[ipart-1][l-1]->GetRandom(); 
-    //return ret;
-    //********************************************
-  } else { //kDefault
-    Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
-    if(ret>e) return e;
     return ret;
   }
+  //kDefault
+  Double_t ret=fHistos[ipart-1][l-1]->GetRandom();
+  if(ret>e) return e;
+  return ret;
 }
 
 Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, Double_t e) const
@@ -766,8 +907,13 @@ Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, Double_t length, D
   return e-loss;
 }
 
-Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, TH1F *hell, Double_t e=1.e6) 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.;
@@ -785,6 +931,265 @@ Double_t AliQuenchingWeights::CalcQuenchedEnergy(Int_t ipart, TH1F *hell, Double
   return e-loss;
 }
 
+Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
+{
+  // return DeltaE for new dynamic version
+  // for given parton type, I0 and I1 value and energy
+  // Dependant on ECM (energy constraint method)
+  // e is used to determine where to set bins to zero.
+
+  // read-in data before first call
+  if(!fTablesLoaded){
+    Fatal("GetELossRandomK","Tables are not loaded.");
+    return -1000.;
+  }
+  if((ipart<1) || (ipart>2)) {
+    Fatal("GetELossRandomK","ipart =%d; but has to be 1 (quark) or 2 (gluon)",ipart);
+    return -1000.;
+  }
+
+  Double_t r=CalcRk(I0,I1);
+  if(r<0.){
+    Fatal("GetELossRandomK","R should not be negative");
+    return -1000.;
+  }
+  Double_t wc=CalcWCk(I1);
+  if(wc<=0.){
+    Fatal("GetELossRandomK","wc should be greater than zero");
+    return -1000.;
+  }
+  if(SampleEnergyLoss(ipart,r)!=0){
+    Fatal("GetELossRandomK","Could not sample energy loss");
+    return -1000.;
+  }
+
+  if(fECMethod==kReweight){
+    Double_t ret = 2.*e;
+    Int_t ws=0;
+    while(ret>e){
+      ret=fHisto->GetRandom(); 
+      if(++ws==1e6){
+       Warning("GetELossRandomK",
+                "Stopping reweighting; maximum loss assigned after 1e6 trials.");
+       return e;
+      }
+    }
+    return ret;
+  }
+
+  //kDefault
+  Double_t ret=fHisto->GetRandom()*wc;
+  if(ret>e) return e;
+  return ret;
+}
+
+Double_t AliQuenchingWeights::CalcQuenchedEnergyK(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
+{
+  //return quenched parton energy
+  //for given parton type, I0 and I1 value and energy
+
+  Double_t loss=GetELossRandomK(ipart,I0,I1,e);
+  return e-loss;
+}
+
+Double_t AliQuenchingWeights::GetELossRandomKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
+{
+  // return DeltaE for new dynamic version
+  // for given parton type, I0 and I1 value and energy
+  // Dependant on ECM (energy constraint method)
+  // e is used to determine where to set bins to zero.
+  // method is optimized and should only be used if 
+  // all parameters are well within the bounds.
+  // read-in data tables before first call 
+
+  Double_t r=CalcRk(I0,I1);
+  if(r<=0.){
+    return 0.;
+  }
+
+  Double_t wc=CalcWCk(I1);
+  if(wc<=0.){
+    return 0.;
+  }
+
+  return GetELossRandomKFastR(ipart,r,wc,e);
+}
+
+Double_t AliQuenchingWeights::GetELossRandomKFastR(Int_t ipart, Double_t r, Double_t wc, Double_t e)
+{
+  // return DeltaE for new dynamic version
+  // for given parton type, R and wc value and energy
+  // Dependant on ECM (energy constraint method)
+  // e is used to determine where to set bins to zero.
+  // method is optimized and should only be used if 
+  // all parameters are well within the bounds.
+  // read-in data tables before first call 
+
+  if(r>=fgkRMax) {
+    r=fgkRMax-1;
+  }  
+  
+  Double_t discrete=0.;
+  Double_t continuous=0.;
+  Int_t bin=1;
+  Double_t xxxx = fHisto->GetBinCenter(bin);
+  if(fMultSoft)
+    CalcMult(ipart,r,xxxx,continuous,discrete);
+  else
+    CalcSingleHard(ipart,r,xxxx,continuous,discrete);
+
+  if(discrete>=1.0) {
+    return 0.; //no energy loss
+  }
+  if (!fHisto) Init();
+  
+  fHisto->SetBinContent(bin,continuous);
+  Int_t kbinmax=fHisto->FindBin(e/wc);
+  if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
+  if(kbinmax==1) return e; //maximum energy loss
+
+  if(fMultSoft) {
+    for(bin=2; bin<=kbinmax; bin++) {
+      xxxx = fHisto->GetBinCenter(bin);
+      CalcMult(ipart,r,xxxx,continuous,discrete);
+      fHisto->SetBinContent(bin,continuous);
+    }
+  } else {
+    for(bin=2; bin<=kbinmax; bin++) {
+      xxxx = fHisto->GetBinCenter(bin);
+      CalcSingleHard(ipart,r,xxxx,continuous,discrete);
+      fHisto->SetBinContent(bin,continuous);
+    }
+  }
+
+  if(fECMethod==kReweight){
+    fHisto->SetBinContent(kbinmax+1,0);
+    fHisto->Fill(0.,discrete*fgkBins/fgkMaxBin);
+  } else if (fECMethod==kReweightCont) {
+    fHisto->SetBinContent(kbinmax+1,0);
+    const Double_t kdelta=fHisto->Integral(1,kbinmax);
+    fHisto->Scale(1./kdelta*(1-discrete));
+    fHisto->Fill(0.,discrete);
+  } else {
+    const Double_t kdelta=fHisto->Integral(1,kbinmax);
+    Double_t val=discrete*fgkBins/fgkMaxBin;
+    fHisto->Fill(0.,val);
+    fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
+  }
+  for(bin=kbinmax+2; bin<=fgkBins; bin++) {
+    fHisto->SetBinContent(bin,0);
+  }
+  //cout << kbinmax << " " << discrete << " " << fHisto->Integral() << endl;
+  Double_t ret=fHisto->GetRandom()*wc;
+  if(ret>e) return e;
+  return ret;
+}
+
+Double_t AliQuenchingWeights::CalcQuenchedEnergyKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e)
+{
+  //return quenched parton energy (fast method)
+  //for given parton type, I0 and I1 value and energy
+
+  Double_t loss=GetELossRandomKFast(ipart,I0,I1,e);
+  return e-loss;
+}
+
+Double_t AliQuenchingWeights::GetDiscreteWeight(Int_t ipart, Double_t I0, Double_t I1)
+{
+  // return discrete weight
+
+  Double_t r=CalcRk(I0,I1);
+  if(r<=0.){
+    return 1.;
+  }
+  return GetDiscreteWeightR(ipart,r);
+}
+
+Double_t AliQuenchingWeights::GetDiscreteWeightR(Int_t ipart, Double_t r)
+{
+  // return discrete weight
+
+  if(r>=fgkRMax) {
+    r=fgkRMax-1;
+  }  
+
+  Double_t discrete=0.;
+  Double_t continuous=0.;
+  Int_t bin=1;
+  Double_t xxxx = fHisto->GetBinCenter(bin);
+  if(fMultSoft)
+    CalcMult(ipart,r,xxxx,continuous,discrete);
+  else
+    CalcSingleHard(ipart,r,xxxx,continuous,discrete);
+  return discrete;
+}
+
+void AliQuenchingWeights::GetZeroLossProb(Double_t &p,Double_t &prw,Double_t &prwcont,
+                                         Int_t ipart,Double_t I0,Double_t I1,Double_t e)
+{
+  //calculate the probabilty that there is no energy
+  //loss for different ways of energy constraint
+  p=1.;prw=1.;prwcont=1.;
+  Double_t r=CalcRk(I0,I1);
+  if(r<=0.){
+    return;
+  }
+  Double_t wc=CalcWCk(I1);
+  if(wc<=0.){
+    return;
+  }
+  GetZeroLossProbR(p,prw,prwcont,ipart,r,wc,e);
+}
+
+void AliQuenchingWeights::GetZeroLossProbR(Double_t &p,Double_t &prw,Double_t &prwcont,
+                                          Int_t ipart, Double_t r,Double_t wc,Double_t e)
+{
+  //calculate the probabilty that there is no energy
+  //loss for different ways of energy constraint
+  if(r>=fgkRMax) {
+    r=fgkRMax-1;
+  }  
+
+  Double_t discrete=0.;
+  Double_t continuous=0.;
+  if (!fHisto) Init();
+  Int_t kbinmax=fHisto->FindBin(e/wc);
+  if(kbinmax>=fgkBins) kbinmax=fgkBins-1;
+  if(fMultSoft) {
+    for(Int_t bin=1; bin<=kbinmax; bin++) {
+      Double_t xxxx = fHisto->GetBinCenter(bin);
+      CalcMult(ipart,r,xxxx,continuous,discrete);
+      fHisto->SetBinContent(bin,continuous);
+    }
+  } else {
+    for(Int_t bin=1; bin<=kbinmax; bin++) {
+      Double_t xxxx = fHisto->GetBinCenter(bin);
+      CalcSingleHard(ipart,r,xxxx,continuous,discrete);
+      fHisto->SetBinContent(bin,continuous);
+    }
+  }
+
+  //non-reweighted P(Delta E = 0)
+  const Double_t kdelta=fHisto->Integral(1,kbinmax);
+  Double_t val=discrete*fgkBins/fgkMaxBin;
+  fHisto->Fill(0.,val);
+  fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta);
+  Double_t hint=fHisto->Integral(1,kbinmax+1);
+  p=fHisto->GetBinContent(1)/hint;
+
+  // reweighted
+  hint=fHisto->Integral(1,kbinmax);
+  prw=fHisto->GetBinContent(1)/hint;
+
+  Double_t xxxx = fHisto->GetBinCenter(1);
+  CalcMult(ipart,r,xxxx,continuous,discrete);
+  fHisto->SetBinContent(1,continuous);
+  hint=fHisto->Integral(1,kbinmax);
+  fHisto->Scale(1./hint*(1-discrete));
+  fHisto->Fill(0.,discrete);
+  prwcont=fHisto->GetBinContent(1);
+}
+
 Int_t AliQuenchingWeights::SampleEnergyLoss() 
 {
   // Has to be called to fill the histograms
@@ -800,14 +1205,14 @@ Int_t AliQuenchingWeights::SampleEnergyLoss()
 
   // read-in data before first call
   if(!fTablesLoaded){
-    Error("CalcMult","Tables are not loaded.");
+    Error("SampleEnergyLoss","Tables are not loaded.");
     return -1;
   }
 
   if(fMultSoft) {
     Int_t lmax=CalcLengthMax(fQTransport);
     if(fLengthMax>lmax){
-      Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,gkRMax);
+      Info("SampleEnergyLoss","Maximum length changed from %d to %d;\nin order to have R < %.f",fLengthMax,lmax,fgkRMax);
       fLengthMax=lmax;
     }
   } else {
@@ -816,8 +1221,8 @@ Int_t AliQuenchingWeights::SampleEnergyLoss()
 
   Reset();
   fHistos=new TH1F**[2];
-  fHistos[0]=new TH1F*[fLengthMax];
-  fHistos[1]=new TH1F*[fLengthMax];
+  fHistos[0]=new TH1F*[4*fLengthMax];
+  fHistos[1]=new TH1F*[4*fLengthMax];
   fLengthMaxOld=fLengthMax; //remember old value in case 
                             //user wants to reset
 
@@ -825,49 +1230,105 @@ Int_t AliQuenchingWeights::SampleEnergyLoss()
   Char_t meddesc[100];
   if(fMultSoft) {
     medvalue=(Int_t)(fQTransport*1000.);
-    sprintf(meddesc,"MS");
+    snprintf(meddesc, 100, "MS");
   } else {
     medvalue=(Int_t)(fMu*1000.);
-    sprintf(meddesc,"SH");
+    snprintf(meddesc, 100, "SH");
   }
 
   for(Int_t ipart=1;ipart<=2;ipart++){
-    for(Int_t l=1;l<=fLengthMax;l++){
-
+    for(Int_t l=1;l<=4*fLengthMax;l++){
       Char_t hname[100];
-      sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
-      Double_t wc = CalcWC(l);
-      fHistos[ipart-1][l-1] = new TH1F(hname,hname,1100,0.,1.1*wc);
+      snprintf(hname, 100, "hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
+      Double_t len=l/4.;
+      Double_t wc = CalcWC(len);
+      fHistos[ipart-1][l-1] = new TH1F(hname,hname,fgkBins,0.,fgkMaxBin*wc);
       fHistos[ipart-1][l-1]->SetXTitle("#Delta E [GeV]");
       fHistos[ipart-1][l-1]->SetYTitle("p(#Delta E)");
       fHistos[ipart-1][l-1]->SetLineColor(4);
 
-      Double_t rrrr = CalcR(wc,l);
+      Double_t rrrr = CalcR(wc,len);
       Double_t discrete=0.;
       // loop on histogram channels
-      for(Int_t bin=1; bin<=1100; bin++) {
+      for(Int_t bin=1; bin<=fgkBins; bin++) {
        Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc;
        Double_t continuous;
-       CalcMult(ipart,rrrr,xxxx,continuous,discrete);
+       if(fMultSoft)
+         CalcMult(ipart,rrrr,xxxx,continuous,discrete);
+       else
+         CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
        fHistos[ipart-1][l-1]->SetBinContent(bin,continuous);
       }
       // add discrete part to distribution
       if(discrete>=1.)
-       for(Int_t bin=2;bin<=1100;bin++) 
+       for(Int_t bin=2;bin<=fgkBins;bin++) 
          fHistos[ipart-1][l-1]->SetBinContent(bin,0.);
       else {
-       Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,1100);
+       Double_t val=discrete/(1.-discrete)*fHistos[ipart-1][l-1]->Integral(1,fgkBins);
        fHistos[ipart-1][l-1]->Fill(0.,val);
       }
-      Double_t hint=fHistos[ipart-1][l-1]->Integral(1,1100);
+      Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins);
       fHistos[ipart-1][l-1]->Scale(1./hint);
     }
   }
   return 0;
 }
 
-const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Int_t l) const
+Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t r)
+{
+  // Sample energy loss directly for one particle type
+  // choses R (safe it and keep it until next call of function)
+
+  // read-in data before first call
+  if(!fTablesLoaded){
+    Error("SampleEnergyLoss","Tables are not loaded.");
+    return -1;
+  }
+
+  Double_t discrete=0.;
+  Double_t continuous=0;;
+  Int_t bin=1;
+  if (!fHisto) Init();
+  Double_t xxxx = fHisto->GetBinCenter(bin);
+  if(fMultSoft)
+    CalcMult(ipart,r,xxxx,continuous,discrete);
+  else
+    CalcSingleHard(ipart,r,xxxx,continuous,discrete);
+
+  if(discrete>=1.) {
+    fHisto->SetBinContent(1,1.);
+    for(bin=2;bin<=fgkBins;bin++) 
+      fHisto->SetBinContent(bin,0.);
+    return 0;
+  }
+
+  fHisto->SetBinContent(bin,continuous);
+  for(bin=2; bin<=fgkBins; bin++) {
+    xxxx = fHisto->GetBinCenter(bin);
+    if(fMultSoft)
+      CalcMult(ipart,r,xxxx,continuous,discrete);
+    else
+      CalcSingleHard(ipart,r,xxxx,continuous,discrete);
+    fHisto->SetBinContent(bin,continuous);
+  }
+
+  Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
+  fHisto->Fill(0.,val);
+  Double_t hint=fHisto->Integral(1,fgkBins);
+  if(hint!=0)
+    fHisto->Scale(1./hint);
+  else {
+    //cout << discrete << " " << hint << " " << continuous << endl;
+    return -1;
+  }
+  return 0;
+}
+
+const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Double_t length) const
 {
+  //return quenching histograms 
+  //for ipart and length
+
   if(!fHistos){
     Fatal("GetELossRandom","Call SampleEnergyLoss method before!");
     return 0;
@@ -877,9 +1338,8 @@ const TH1F* AliQuenchingWeights::GetHisto(Int_t ipart,Int_t l) const
     return 0;
   }
 
+  Int_t l=GetIndex(length);
   if(l<=0) return 0;
-  if(l>fLengthMax) l=fLengthMax;
-
   return fHistos[ipart-1][l-1];
 }
 
@@ -890,34 +1350,30 @@ TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t l
   //        b) mu      = Debye mass (GeV)
   // length = path length in medium (fm)
   // Get from SW tables:
-  // - discrete weight (probability to have NO energy loss)
-  // - continuous weight, as a function of dE
-  // compute up to max dE = 1.1*wc 
-  //   (as an histogram named hContQW_<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;
@@ -939,27 +1395,23 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t
   //        b) mu      = Debye mass (GeV)
   // length = path length in medium (fm)
   // Get from SW tables:
-  // - discrete weight (probability to have NO energy loss)
-  // - continuous weight, as a function of dE
-  // compute up to max dE = 1.1*wc 
-  //   (as an histogram named hContQW_<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})");
@@ -968,8 +1420,8 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t
   histx->SetLineColor(4);
 
   Double_t rrrr = CalcR(wc,length);
-  // loop on histogram channels
-  for(Int_t bin=1; bin<=1100; bin++) {
+  //loop on histogram channels
+  for(Int_t bin=1; bin<=fgkBins; bin++) {
     Double_t xxxx = histx->GetBinCenter(bin);
     Double_t continuous,discrete;
     Int_t ret=0;
@@ -984,8 +1436,63 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t
   return histx;
 }
 
-TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e=1.e6) const
+TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t r) const 
 {
+  // compute P(E) distribution for
+  // given ipart = 1 for quark, 2 for gluon 
+  // and R
+
+  Char_t meddesc[100];
+  if(fMultSoft) {
+    snprintf(meddesc, 100, "MS");
+  } else {
+    snprintf(meddesc, 100, "SH");
+  }
+
+  Char_t hname[100];
+  snprintf(hname, 100, "hQWHistox_%s_%d_%.2f",meddesc,ipart,r);
+  TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
+  histx->SetXTitle("x = #Delta E/#omega_{c}");
+  if(fMultSoft)
+    histx->SetYTitle("p(#Delta E/#omega_{c})");
+  else
+    histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
+  histx->SetLineColor(4);
+
+  Double_t rrrr = r;
+  Double_t continuous=0.,discrete=0.;
+  //loop on histogram channels
+  for(Int_t bin=1; bin<=fgkBins; bin++) {
+    Double_t xxxx = histx->GetBinCenter(bin);
+    Int_t ret=0;
+    if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
+    else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
+    if(ret!=0){
+      delete histx;
+      return 0;
+    };
+    histx->SetBinContent(bin,continuous);
+  }
+
+  //add discrete part to distribution
+  if(discrete>=1.)
+    for(Int_t bin=2;bin<=fgkBins;bin++) 
+      histx->SetBinContent(bin,0.);
+  else {
+    Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins);
+    histx->Fill(0.,val);
+  }
+  Double_t hint=histx->Integral(1,fgkBins);
+  if(hint!=0) histx->Scale(1./hint);
+
+  return histx;
+}
+
+TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
+{
+  // compute energy loss histogram for 
+  // parton type, medium value, length and energy
+
   AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
   if(fMultSoft){
     dummy->SetQTransport(medval);
@@ -999,11 +1506,11 @@ TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_
   Char_t name[100];
   Char_t hname[100];
   if(ipart==1){
-    sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); 
-    sprintf(hname,"hLossQuarks"); 
+    snprintf(name, 100, "Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); 
+    snprintf(hname,100, "hLossQuarks"); 
   } else {
-    sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); 
-    sprintf(hname,"hLossGluons"); 
+    snprintf(name, 100, "Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); 
+    snprintf(hname, 100, "hLossGluons"); 
   }
 
   TH1F *h = new TH1F(hname,name,250,0,250);
@@ -1012,15 +1519,17 @@ TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_
     Double_t loss=dummy->GetELossRandom(ipart,l,e);
     h->Fill(loss);
   }
-
   h->SetStats(kTRUE);
-
   delete dummy;
   return h;
 }
 
-TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e=1.e6) const
+TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
 {
+  // compute energy loss histogram for 
+  // parton type, medium value, 
+  // length distribution and energy
+
   AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
   if(fMultSoft){
     dummy->SetQTransport(medval);
@@ -1034,11 +1543,11 @@ TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *h
   Char_t name[100];
   Char_t hname[100];
   if(ipart==1){
-    sprintf(name,"Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); 
-    sprintf(hname,"hLossQuarks"); 
+    snprintf(name, 100, "Energy Loss Distribution - Quarks;E_{loss} (GeV);#"); 
+    snprintf(hname,100, "hLossQuarks"); 
   } else {
-    sprintf(name,"Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); 
-    sprintf(hname,"hLossGluons"); 
+    snprintf(name,100, "Energy Loss Distribution - Gluons;E_{loss} (GeV);#"); 
+    snprintf(hname, 100, "hLossGluons"); 
   }
 
   TH1F *h = new TH1F(hname,name,250,0,250);
@@ -1047,15 +1556,72 @@ TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *h
     Double_t loss=dummy->GetELossRandom(ipart,hEll,e);
     h->Fill(loss);
   }
-
   h->SetStats(kTRUE);
-
   delete dummy;
   return h;
 }
 
-void AliQuenchingWeights::PlotDiscreteWeights(Int_t len) const
+TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t r) const 
+{
+  // compute energy loss histogram for 
+  // parton type and given R
+
+  TH1F *dummy = ComputeQWHistoX(ipart,r);
+  if(!dummy) return 0;
+
+  Char_t hname[100];
+  snprintf(hname, 100, "hELossHistox_%d_%.2f",ipart,r);
+  TH1F *histx = new TH1F("histxr",hname,fgkBins,0.,fgkMaxBin);
+  for(Int_t i=0;i<100000;i++){
+    //if(i % 1000 == 0) cout << "." << flush;
+    Double_t loss=dummy->GetRandom();
+    histx->Fill(loss);
+  }
+  delete dummy;
+  return histx;
+}
+
+Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,Double_t l) const
+{
+  // compute average energy loss for 
+  // parton type, medium value, length and energy
+
+  TH1F *dummy = ComputeELossHisto(ipart,medval,l);
+  if(!dummy) return 0;
+  Double_t ret=dummy->GetMean();
+  delete dummy;
+  return ret;
+}
+
+Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEll) const
+{
+  // compute average energy loss for 
+  // parton type, medium value, 
+  // length distribution and energy
+
+  TH1F *dummy = ComputeELossHisto(ipart,medval,hEll);
+  if(!dummy) return 0;
+  Double_t ret=dummy->GetMean();
+  delete dummy;
+  return ret;
+}
+
+Double_t  AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t r) const 
 {
+  // compute average energy loss over wc 
+  // for parton type and given R
+
+  TH1F *dummy = ComputeELossHisto(ipart,r);
+  if(!dummy) return 0;
+  Double_t ret=dummy->GetMean();
+  delete dummy;
+  return ret;
+}
+
+void AliQuenchingWeights::PlotDiscreteWeights(Double_t len,Double_t qm) const
+{
+  // plot discrete weights for given length
+
   TCanvas *c; 
   if(fMultSoft) 
     c = new TCanvas("cdiscms","Discrete Weight for Multiple Scattering",0,0,500,400);
@@ -1063,85 +1629,96 @@ void AliQuenchingWeights::PlotDiscreteWeights(Int_t len) const
     c = new TCanvas("cdiscsh","Discrete Weight for Single Hard Scattering",0,0,500,400);
   c->cd();
 
-  TH2F *hframe = new TH2F("hdisc","",2,0,5.1,2,0,1);
+  TH2F *hframe = new TH2F("hdisc","",2,0,qm+.1,2,0,1.25);
   hframe->SetStats(0);
   if(fMultSoft) 
     hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
   else
     hframe->SetXTitle("#mu [GeV]");
-  hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
+  //hframe->SetYTitle("Probability #Delta E = 0 , p_{0}");
+  hframe->SetYTitle("p_{0} (discrete weight)");
   hframe->Draw();
 
-  TGraph *gq=new TGraph(20);
+  Int_t points=(Int_t)qm*4;
+  TGraph *gq=new TGraph(points);
   Int_t i=0;
   if(fMultSoft) {
-    for(Double_t q=0.05;q<=5.05;q+=0.25){
+    for(Double_t q=0.05;q<=qm+.05;q+=0.25){
       Double_t disc,cont;
-      CalcMult(1,1.0, q, len, cont, disc);
+      CalcMult(1,1.0,q,len,cont,disc);
       gq->SetPoint(i,q,disc);i++;
     }
   } else {
-    for(Double_t m=0.05;m<=5.05;m+=0.25){
+    for(Double_t m=0.05;m<=qm+.05;m+=0.25){
       Double_t disc,cont;
-      CalcSingleHard(1,1.0, m, len, cont, disc);
+      CalcSingleHard(1,1.0,m,len,cont, disc);
       gq->SetPoint(i,m,disc);i++;
     }
   }
   gq->SetMarkerStyle(20);
-  gq->Draw("pl");
+  gq->SetMarkerColor(1);
+  gq->SetLineStyle(1);
+  gq->SetLineColor(1);
+  gq->Draw("l");
 
-  TGraph *gg=new TGraph(20);
+  TGraph *gg=new TGraph(points);
   i=0;
   if(fMultSoft){
-    for(Double_t q=0.05;q<=5.05;q+=0.25){
+    for(Double_t q=0.05;q<=qm+.05;q+=0.25){
       Double_t disc,cont;
-      CalcMult(2,1.0, q, 5., cont, disc);
+      CalcMult(2,1.0,q,len,cont,disc);
       gg->SetPoint(i,q,disc);i++;
     }
   } else {
-    for(Double_t m=0.05;m<=5.05;m+=0.25){
+    for(Double_t m=0.05;m<=qm+.05;m+=0.25){
       Double_t disc,cont;
-      CalcSingleHard(2,1.0, m, 5., cont, disc);
+      CalcSingleHard(2,1.0,m,len,cont,disc);
       gg->SetPoint(i,m,disc);i++;
     }
   }
   gg->SetMarkerStyle(24);
-  gg->Draw("pl");
+  gg->SetMarkerColor(2);
+  gg->SetLineStyle(2);
+  gg->SetLineColor(2);
+  gg->Draw("l");
 
   TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
   l1a->SetFillStyle(0);
   l1a->SetBorderSize(0);
   Char_t label[100];
-  sprintf(label,"L = %d fm",len);
+  snprintf(label, 100, "L = %.1f fm",len);
   l1a->AddEntry(gq,label,"");
-  l1a->AddEntry(gq,"quark","pl");
-  l1a->AddEntry(gg,"gluon","pl");
+  l1a->AddEntry(gq,"quark projectile","l");
+  l1a->AddEntry(gg,"gluon projectile","l");
   l1a->Draw();
 
   c->Update();
 }
 
-void AliQuenchingWeights::PlotContWeights(Int_t itype,Int_t ell) const
+void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
 {
+  // plot continous weights for 
+  // given parton type and length
+
   Float_t medvals[3];
   Char_t title[1024];
   Char_t name[1024];
   if(fMultSoft) {
     if(itype==1)
-      sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
+      snprintf(title, 1024, "Cont. Weight for Multiple Scattering - Quarks");
     else if(itype==2)
-      sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
+      snprintf(title, 1024, "Cont. Weight for Multiple Scattering - Gluons");
     else return;
     medvals[0]=4;medvals[1]=1;medvals[2]=0.5;
-    sprintf(name,"ccont-ms-%d",itype);
+    snprintf(name, 1024, "ccont-ms-%d",itype);
   } else {
     if(itype==1)
-      sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
+      snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Quarks");
     else if(itype==2)
-      sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
+      snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Gluons");
     else return;
     medvals[0]=2;medvals[1]=1;medvals[2]=0.5;
-    sprintf(name,"ccont-ms-%d",itype);
+    snprintf(name, 1024, "ccont-ms-%d",itype);
   }
 
   TCanvas *c = new TCanvas(name,title,0,0,500,400);
@@ -1165,21 +1742,21 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Int_t ell) const
   l1a->SetFillStyle(0);
   l1a->SetBorderSize(0);
   Char_t label[100];
-  sprintf(label,"L = %d fm",ell);
+  snprintf(label, 100, "L = %.1f fm",ell);
   l1a->AddEntry(h1,label,"");
   if(fMultSoft) {
-    sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
+    snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[0]);
     l1a->AddEntry(h1,label,"pl");
-    sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
+    snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[1]);
     l1a->AddEntry(h2,label,"pl");
-    sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
-    l1a->AddEntry(h3,label,"pl");
+    snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medvals[2]);
+    l1a->AddEntry(h3, label,"pl");
   } else {
-    sprintf(label,"#mu = %.1f GeV",medvals[0]);
+    snprintf(label, 100, "#mu = %.1f GeV",medvals[0]);
     l1a->AddEntry(h1,label,"pl");
-    sprintf(label,"#mu = %.1f GeV",medvals[1]);
+    snprintf(label, 100, "#mu = %.1f GeV",medvals[1]);
     l1a->AddEntry(h2,label,"pl");
-    sprintf(label,"#mu = %.1f GeV",medvals[2]);
+    snprintf(label, 100, "#mu = %.1f GeV",medvals[2]);
     l1a->AddEntry(h3,label,"pl");
   }
   l1a->Draw();
@@ -1187,24 +1764,27 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Int_t ell) const
   c->Update();
 }
 
-void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const
+void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
 {
+  // plot continous weights for 
+  // given parton type and medium value
+
   Char_t title[1024];
   Char_t name[1024];
   if(fMultSoft) {
     if(itype==1)
-      sprintf(title,"Cont. Weight for Multiple Scattering - Quarks");
+      snprintf(title,1024, "Cont. Weight for Multiple Scattering - Quarks");
     else if(itype==2)
-      sprintf(title,"Cont. Weight for Multiple Scattering - Gluons");
+      snprintf(title,1024, "Cont. Weight for Multiple Scattering - Gluons");
     else return;
-    sprintf(name,"ccont2-ms-%d",itype);
+    snprintf(name,1024, "ccont2-ms-%d",itype);
   } else {
     if(itype==1)
-      sprintf(title,"Cont. Weight for Single Hard Scattering - Quarks");
+      snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Quarks");
     else if(itype==2)
-      sprintf(title,"Cont. Weight for Single Hard Scattering - Gluons");
+      snprintf(title, 1024, "Cont. Weight for Single Hard Scattering - Gluons");
     else return;
-    sprintf(name,"ccont2-sh-%d",itype);
+    snprintf(name, 1024, "ccont2-sh-%d",itype);
   }
   TCanvas *c = new TCanvas(name,title,0,0,500,400);
   c->cd();
@@ -1228,9 +1808,9 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const
   l1a->SetBorderSize(0);
   Char_t label[100];
   if(fMultSoft)
-    sprintf(label,"#hat{q} = %.1f GeV^{2}/fm",medval);
+    snprintf(label, 100, "#hat{q} = %.1f GeV^{2}/fm",medval);
   else
-    sprintf(label,"#mu = %.1f GeV",medval);
+    snprintf(label, 100, "#mu = %.1f GeV",medval);
 
   l1a->AddEntry(h1,label,"");
   l1a->AddEntry(h1,"L = 8 fm","pl");
@@ -1241,26 +1821,29 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const
   c->Update();
 }
 
-void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e=1.e6)  const
+void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t qm,Double_t e) const
 {
+  // plot average energy loss for given length
+  // and parton energy 
+
   if(!fTablesLoaded){
-    Error("CalcMult","Tables are not loaded.");
+    Error("PlotAvgELoss","Tables are not loaded.");
     return;
   }
 
   Char_t title[1024];
   Char_t name[1024];
   if(fMultSoft){ 
-    sprintf(title,"Average Loss for Multiple Scattering");
-    sprintf(name,"cavgelossms");
+    snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
+    snprintf(name, 1024, "cavgelossms");
   } else {
-    sprintf(title,"Average Loss for Single Hard Scattering");
-    sprintf(name,"cavgelosssh");
+    snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering");
+    snprintf(name, 1024, "cavgelosssh");
   }
 
   TCanvas *c = new TCanvas(name,title,0,0,500,400);
   c->cd();
-  TH2F *hframe = new TH2F("avgloss",title,2,0,5.1,2,0,100);
+  TH2F *hframe = new TH2F("avgloss","",2,0,qm+.1,2,0,100);
   hframe->SetStats(0);
   if(fMultSoft) 
     hframe->SetXTitle("#hat{q} [GeV^{2}/fm]");
@@ -1271,54 +1854,71 @@ void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e=1.e6)  const
 
   TGraph *gq=new TGraph(20);
   Int_t i=0;
-  for(Double_t v=0.05;v<=5.05;v+=0.25){
+  for(Double_t v=0.05;v<=qm+.05;v+=0.25){
     TH1F *dummy=ComputeELossHisto(1,v,len,e);
     Double_t avgloss=dummy->GetMean();
     gq->SetPoint(i,v,avgloss);i++;
     delete dummy;
   }
-  gq->SetMarkerStyle(20);
+  gq->SetMarkerStyle(21);
   gq->Draw("pl");
 
-  TGraph *gg=new TGraph(20);
+  Int_t points=(Int_t)qm*4;
+  TGraph *gg=new TGraph(points);
   i=0;
-  for(Double_t v=0.05;v<=5.05;v+=0.25){
+  for(Double_t v=0.05;v<=qm+.05;v+=0.25){
     TH1F *dummy=ComputeELossHisto(2,v,len,e);
     Double_t avgloss=dummy->GetMean();
     gg->SetPoint(i,v,avgloss);i++;
     delete dummy;
   }
-  gg->SetMarkerStyle(24);
+  gg->SetMarkerStyle(20);
+  gg->SetMarkerColor(2);
   gg->Draw("pl");
 
-  TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
+  TGraph *gratio=new TGraph(points);
+  for(i=0;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=1.e6) const
+void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
 {
+  // plot average energy loss for given
+  // length distribution and parton energy
+
   if(!fTablesLoaded){
-    Error("CalcMult","Tables are not loaded.");
+    Error("PlotAvgELossVs","Tables are not loaded.");
     return;
   }
 
   Char_t title[1024];
   Char_t name[1024];
   if(fMultSoft){ 
-    sprintf(title,"Average Loss for Multiple Scattering");
-    sprintf(name,"cavgelossms2");
+    snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
+    snprintf(name, 1024, "cavgelossms2");
   } else {
-    sprintf(title,"Average Loss for Single Hard Scattering");
-    sprintf(name,"cavgelosssh2");
+    snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering");
+    snprintf(name, 1024, "cavgelosssh2");
   }
 
   TCanvas *c = new TCanvas(name,title,0,0,500,400);
@@ -1354,34 +1954,135 @@ void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e=1.e6) const
   gg->SetMarkerStyle(24);
   gg->Draw("pl");
 
+  TGraph *gratio=new TGraph(20);
+  for(i=0;i<20;i++){
+    Double_t x,y,x2,y2;
+    gg->GetPoint(i,x,y);
+    gq->GetPoint(i,x2,y2);
+    if(y2>0)
+      gratio->SetPoint(i,x,y/y2*10/2.25);
+    else gratio->SetPoint(i,x,0);
+  }
+  gratio->SetLineStyle(4);
+  //gratio->Draw();
+
   TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
   l1a->SetFillStyle(0);
   l1a->SetBorderSize(0);
   Char_t label[100];
-  sprintf(label,"<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);
@@ -1418,7 +2119,7 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Int_t len) const
   l1a->SetFillStyle(0);
   l1a->SetBorderSize(0);
   Char_t label[100];
-  sprintf(label,"L = %d fm",len);
+  snprintf(label, 100, "L = %.1f fm",len);
   l1a->AddEntry(gq,label,"");
   l1a->AddEntry(gq,"quark","pl");
   l1a->AddEntry(gg,"gluon","pl");
@@ -1429,19 +2130,22 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Int_t len) const
 
 void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
 {
+  // plot relative energy loss for given
+  // length distribution and parton energy versus pt
+
   if(!fTablesLoaded){
-    Error("CalcMult","Tables are not loaded.");
+    Error("PlotAvgELossVsPt","Tables are not loaded.");
     return;
   }
 
   Char_t title[1024];
   Char_t name[1024];
   if(fMultSoft){
-    sprintf(title,"Relative Loss for Multiple Scattering");
-    sprintf(name,"cavgelossvsptms2");
+    snprintf(title, 1024, "Relative Energy Loss for Multiple Scattering");
+    snprintf(name,  1024, "cavgelossvsptms2");
   } else {
-    sprintf(title,"Relative Loss for Single Hard Scattering");
-    sprintf(name,"cavgelossvsptsh2");
+    snprintf(title, 1024, "Relative Energy Loss for Single Hard Scattering");
+    snprintf(name,  1024, "cavgelossvsptsh2");
   }
   TCanvas *c = new TCanvas(name,title,0,0,500,400);
   c->cd();
@@ -1477,7 +2181,7 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
   l1a->SetFillStyle(0);
   l1a->SetBorderSize(0);
   Char_t label[100];
-  sprintf(label,"<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");
@@ -1486,3 +2190,13 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
   c->Update();
 }
 
+Int_t AliQuenchingWeights::GetIndex(Double_t len) const
+{
+  //get the index according to length
+  if(len>fLengthMax) len=fLengthMax;
+
+  Int_t l=Int_t(len/0.25);
+  if((len-l*0.25)>0.125) l++;
+  return l;
+}
+