Coding rule violations.
[u/mrichter/AliRoot.git] / FASTSIM / AliQuenchingWeights.cxx
index 5fad38c..b8e0c2a 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
-//
-//  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]             
+
+// 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            
 //
-//----------------------------------------------------------------------------
+//=================== 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"
 
@@ -47,39 +52,64 @@ const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197;
 // maximum value of R
 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::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)
 {
   //default constructor 
 
-  fTablesLoaded=kFALSE;
-  fMultSoft=kTRUE; 
-  fHistos=0;
-  SetMu();
-  SetQTransport();
-  SetECMethod();
-  SetLengthMax();
-  fLengthMaxOld=0;
-  fInstanceNumber=fgCounter++;
 }
 
 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 
+  // 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=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
 }
@@ -87,6 +117,16 @@ 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()
@@ -94,7 +134,7 @@ 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];
   }
@@ -110,8 +150,9 @@ 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
+  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) 
@@ -123,7 +164,7 @@ Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
   fMultSoft=kTRUE;
   
   Char_t fname[1024];
-  sprintf(fname,"%s",gSystem->ExpandPathName(contall));
+  snprintf(fname,1024, "%s",gSystem->ExpandPathName(contall));
   //PH  ifstream fincont(fname);
   fstream fincont(fname,ios::in);
 #if defined(__HP_aCC) || defined(__DECCXX)
@@ -159,7 +200,7 @@ Int_t AliQuenchingWeights::InitMult(const Char_t *contall,const Char_t *discall)
   }
   fincont.close();
 
-  sprintf(fname,"%s",gSystem->ExpandPathName(discall));
+  snprintf(fname,1024, "%s",gSystem->ExpandPathName(discall));
   //PH  ifstream findisc(fname); 
   fstream findisc(fname,ios::in); 
 #if defined(__HP_aCC) || defined(__DECCXX)
@@ -375,7 +416,7 @@ Int_t AliQuenchingWeights::CalcMult(Int_t ipart, Double_t rrrr,Double_t xxxx,
   continuous=0.;
   discrete=0.;
 
-  // read-in data before first call
+  //read-in data before first call
   if(!fTablesLoaded){
     Error("CalcMult","Tables are not loaded.");
     return -1;
@@ -449,7 +490,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];
@@ -469,7 +510,7 @@ Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *di
   fMultSoft=kFALSE;
   
   Char_t fname[1024];
-  sprintf(fname,"%s",gSystem->ExpandPathName(contall));
+  snprintf(fname, 1024, "%s",gSystem->ExpandPathName(contall));
   //PH  ifstream fincont(fname);
   fstream fincont(fname,ios::in);
 #if defined(__HP_aCC) || defined(__DECCXX)
@@ -512,7 +553,7 @@ Int_t AliQuenchingWeights::InitSingleHard(const Char_t *contall,const Char_t *di
   }
   fincont.close();
 
-  sprintf(fname,"%s",gSystem->ExpandPathName(discall));
+  snprintf(fname, 1024, "%s",gSystem->ExpandPathName(discall));
   //PH  ifstream findisc(fname); 
   fstream findisc(fname,ios::in); 
 #if defined(__HP_aCC) || defined(__DECCXX)
@@ -698,11 +739,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;
   }
 
@@ -791,15 +832,30 @@ Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
 
 Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const 
 { 
+  //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 = wc*l*fgkConvFmToInvGeV;
-  if(R>fgkRMax) {
-    Warning("CalcR","Value of R = %.2f; should be less than %.2f",R,fgkRMax);
-    return -R;
+  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
@@ -815,34 +871,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
@@ -878,6 +929,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
@@ -893,7 +1203,7 @@ 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;
   }
 
@@ -909,8 +1219,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
 
@@ -918,48 +1228,101 @@ 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
@@ -973,9 +1336,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];
 }
 
@@ -986,34 +1348,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;
@@ -1035,27 +1393,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})");
@@ -1064,8 +1418,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;
@@ -1080,10 +1434,62 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t
   return histx;
 }
 
+TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t r) const 
+{
+  // compute P(E) distribution for
+  // given ipart = 1 for quark, 2 for gluon 
+  // and R
+
+  Char_t meddesc[100];
+  if(fMultSoft) {
+    snprintf(meddesc, 100, "MS");
+  } else {
+    snprintf(meddesc, 100, "SH");
+  }
+
+  Char_t hname[100];
+  snprintf(hname, 100, "hQWHistox_%s_%d_%.2f",meddesc,ipart,r);
+  TH1F *histx = new TH1F("histx",hname,fgkBins,0.,fgkMaxBin);
+  histx->SetXTitle("x = #Delta E/#omega_{c}");
+  if(fMultSoft)
+    histx->SetYTitle("p(#Delta E/#omega_{c})");
+  else
+    histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
+  histx->SetLineColor(4);
+
+  Double_t rrrr = r;
+  Double_t continuous=0.,discrete=0.;
+  //loop on histogram channels
+  for(Int_t bin=1; bin<=fgkBins; bin++) {
+    Double_t xxxx = histx->GetBinCenter(bin);
+    Int_t ret=0;
+    if(fMultSoft) ret=CalcMult(ipart,rrrr,xxxx,continuous,discrete);
+    else ret=CalcSingleHard(ipart,rrrr,xxxx,continuous,discrete);
+    if(ret!=0){
+      delete histx;
+      return 0;
+    };
+    histx->SetBinContent(bin,continuous);
+  }
+
+  //add discrete part to distribution
+  if(discrete>=1.)
+    for(Int_t bin=2;bin<=fgkBins;bin++) 
+      histx->SetBinContent(bin,0.);
+  else {
+    Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins);
+    histx->Fill(0.,val);
+  }
+  Double_t hint=histx->Integral(1,fgkBins);
+  if(hint!=0) histx->Scale(1./hint);
+
+  return histx;
+}
+
 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_t l,Double_t e) const
 {
-  //compute energy loss histogram for 
-  //parton type, medium value, length and energy
+  // compute energy loss histogram for 
+  // parton type, medium value, length and energy
 
   AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
   if(fMultSoft){
@@ -1098,11 +1504,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);
@@ -1111,18 +1517,16 @@ TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,Double_
     Double_t loss=dummy->GetELossRandom(ipart,l,e);
     h->Fill(loss);
   }
-
   h->SetStats(kTRUE);
-
   delete dummy;
   return h;
 }
 
 TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *hEll,Double_t e) const
 {
-  //compute energy loss histogram for 
-  //parton type, medium value, 
-  //length distribution and energy
+  // compute energy loss histogram for 
+  // parton type, medium value, 
+  // length distribution and energy
 
   AliQuenchingWeights *dummy=new AliQuenchingWeights(*this);
   if(fMultSoft){
@@ -1137,11 +1541,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);
@@ -1150,16 +1554,71 @@ 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
+  // plot discrete weights for given length
 
   TCanvas *c; 
   if(fMultSoft) 
@@ -1168,88 +1627,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
+  // 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);
@@ -1273,21 +1740,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();
@@ -1295,27 +1762,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
+  // 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();
@@ -1339,9 +1806,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");
@@ -1352,29 +1819,29 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t medval) const
   c->Update();
 }
 
-void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e)  const
+void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t qm,Double_t e) const
 {
-  //plot average energy loss for given length
-  //and parton energy 
+  // 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]");
@@ -1385,34 +1852,48 @@ void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e)  const
 
   TGraph *gq=new TGraph(20);
   Int_t i=0;
-  for(Double_t v=0.05;v<=5.05;v+=0.25){
+  for(Double_t v=0.05;v<=qm+.05;v+=0.25){
     TH1F *dummy=ComputeELossHisto(1,v,len,e);
     Double_t avgloss=dummy->GetMean();
     gq->SetPoint(i,v,avgloss);i++;
     delete dummy;
   }
-  gq->SetMarkerStyle(20);
+  gq->SetMarkerStyle(21);
   gq->Draw("pl");
 
-  TGraph *gg=new TGraph(20);
+  Int_t points=(Int_t)qm*4;
+  TGraph *gg=new TGraph(points);
   i=0;
-  for(Double_t v=0.05;v<=5.05;v+=0.25){
+  for(Double_t v=0.05;v<=qm+.05;v+=0.25){
     TH1F *dummy=ComputeELossHisto(2,v,len,e);
     Double_t avgloss=dummy->GetMean();
     gg->SetPoint(i,v,avgloss);i++;
     delete dummy;
   }
-  gg->SetMarkerStyle(24);
+  gg->SetMarkerStyle(20);
+  gg->SetMarkerColor(2);
   gg->Draw("pl");
 
-  TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
+  TGraph *gratio=new TGraph(points);
+  for(i=0;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();
@@ -1420,22 +1901,22 @@ void AliQuenchingWeights::PlotAvgELoss(Int_t len,Double_t e)  const
 
 void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
 {
-  //plot average energy loss for given
-  //length distribution and parton energy
+  // 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);
@@ -1471,37 +1952,135 @@ void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
   gg->SetMarkerStyle(24);
   gg->Draw("pl");
 
+  TGraph *gratio=new TGraph(20);
+  for(i=0;i<20;i++){
+    Double_t x,y,x2,y2;
+    gg->GetPoint(i,x,y);
+    gq->GetPoint(i,x2,y2);
+    if(y2>0)
+      gratio->SetPoint(i,x,y/y2*10/2.25);
+    else gratio->SetPoint(i,x,0);
+  }
+  gratio->SetLineStyle(4);
+  //gratio->Draw();
+
   TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
   l1a->SetFillStyle(0);
   l1a->SetBorderSize(0);
   Char_t label[100];
-  sprintf(label,"<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 relative energy loss for given
-  //length and parton energy versus pt
+  // 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);
@@ -1538,7 +2117,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");
@@ -1549,22 +2128,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
+  // 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();
@@ -1600,7 +2179,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");
@@ -1609,3 +2188,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;
+}
+