Coding rule violations.
[u/mrichter/AliRoot.git] / FASTSIM / AliQuenchingWeights.cxx
index ce1b27f..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
@@ -29,8 +28,8 @@
 //
 //=================== 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)
+// 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>
@@ -41,6 +40,7 @@
 #include <TGraph.h>
 #include <TROOT.h>
 #include <TSystem.h>
+#include <TString.h>
 #include <TLegend.h>
 #include "AliQuenchingWeights.h"
 
@@ -61,30 +61,36 @@ 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();
-  SetK();
-  fECMethod=kDefault; 
-  //SetECMethod();
-  SetLengthMax();
-  fLengthMaxOld=0;
-  fInstanceNumber=fgCounter++;
-  Char_t name[100];
-  sprintf(name,"hhistoqw_%d",fInstanceNumber);
-  fHisto = new TH1F(name,"",fgkBins,0.,fgkMaxBin);
-  for(Int_t bin=1;bin<=fgkBins;bin++) 
-    fHisto->SetBinContent(bin,0.);
 }
 
 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 
 
@@ -99,7 +105,7 @@ AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a)
   fLengthMax=a.GetLengthMax();
   fInstanceNumber=fgCounter++;
   Char_t name[100];
-  sprintf(name,"hhistoqw_%d",fInstanceNumber);
+  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.);
@@ -114,6 +120,15 @@ AliQuenchingWeights::~AliQuenchingWeights()
   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
@@ -135,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) 
@@ -148,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)
@@ -184,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)
@@ -494,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)
@@ -537,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)
@@ -726,7 +742,7 @@ Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart, Double_t rrrr,Double_t xx
     Error("CalcSingleHard","Tables are not loaded.");
     return -1;
   }
-  if(!fMultSoft){
+  if(fMultSoft){
     Error("CalcSingleHard","Tables are not loaded for Single Hard Scattering.");
     return -1;
   }
@@ -816,15 +832,15 @@ Int_t AliQuenchingWeights::CalcSingleHard(Int_t ipart,
 
 Double_t AliQuenchingWeights::CalcR(Double_t wc, Double_t l) const 
 { 
-  //calculate R value and 
+  //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 = wc*l*fgkConvFmToInvGeV;
+  if(r >= fgkRMax) {
+    Warning("CalcR","Value of r = %.2f; should be less than %.2f", r, fgkRMax);
+    return fgkRMax-1;
   }  
-  return R;
+  return r;
 }
 
 Double_t AliQuenchingWeights::CalcRk(Double_t k, Double_t I0, Double_t I1) const
@@ -832,14 +848,14 @@ 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;
+  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 -R;
+    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
@@ -868,7 +884,7 @@ Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Doubl
       ret=fHistos[ipart-1][l-1]->GetRandom(); 
       if(++ws==1e6){
        Warning("GetELossRandom",
-                "Aborted reweighting; maximum loss assigned after 1e6 trials.");
+                "Stopping reweighting; maximum loss assigned after 1e6 trials.");
        return e;
       }
     }
@@ -930,17 +946,17 @@ Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t
     return -1000.;
   }
 
-  Double_t R=CalcRk(I0,I1);
-  if(R<0){
+  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){
+  if(wc<=0.){
     Fatal("GetELossRandomK","wc should be greater than zero");
     return -1000.;
   }
-  if(SampleEnergyLoss(ipart,R)!=0){
+  if(SampleEnergyLoss(ipart,r)!=0){
     Fatal("GetELossRandomK","Could not sample energy loss");
     return -1000.;
   }
@@ -952,7 +968,7 @@ Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t
       ret=fHisto->GetRandom(); 
       if(++ws==1e6){
        Warning("GetELossRandomK",
-                "Aborted reweighting; maximum loss assigned after 1e6 trials.");
+                "Stopping reweighting; maximum loss assigned after 1e6 trials.");
        return e;
       }
     }
@@ -984,56 +1000,84 @@ Double_t AliQuenchingWeights::GetELossRandomKFast(Int_t ipart, Double_t I0, Doub
   // all parameters are well within the bounds.
   // read-in data tables before first call 
 
-  Double_t R=CalcRk(I0,I1);
-  if(R<=0){
+  Double_t r=CalcRk(I0,I1);
+  if(r<=0.){
     return 0.;
   }
 
   Double_t wc=CalcWCk(I1);
-  if(wc<=0){
+  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;
+  Double_t continuous=0.;
   Int_t bin=1;
   Double_t xxxx = fHisto->GetBinCenter(bin);
   if(fMultSoft)
-    CalcMult(ipart,R,xxxx,continuous,discrete);
+    CalcMult(ipart,r,xxxx,continuous,discrete);
   else
-    CalcSingleHard(ipart,R,xxxx,continuous,discrete);
+    CalcSingleHard(ipart,r,xxxx,continuous,discrete);
 
-  if(discrete>0.999) {
-    return 0.0;
+  if(discrete>=1.0) {
+    return 0.; //no energy loss
   }
-
-  //set first bin
+  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(Int_t bin=2; bin<=fgkBins; bin++) {
+    for(bin=2; bin<=kbinmax; bin++) {
       xxxx = fHisto->GetBinCenter(bin);
-      CalcMult(ipart,R,xxxx,continuous,discrete);
+      CalcMult(ipart,r,xxxx,continuous,discrete);
       fHisto->SetBinContent(bin,continuous);
     }
   } else {
-    for(Int_t bin=2; bin<=fgkBins; bin++) {
+    for(bin=2; bin<=kbinmax; bin++) {
       xxxx = fHisto->GetBinCenter(bin);
-      CalcSingleHard(ipart,R,xxxx,continuous,discrete);
+      CalcSingleHard(ipart,r,xxxx,continuous,discrete);
       fHisto->SetBinContent(bin,continuous);
     }
   }
 
-  // add discrete part to distribution
-  Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins);
-  fHisto->Fill(0.,val);
-
-  if(fECMethod==kReweight) {
-    for(Int_t bin=fHisto->FindBin(e/wc)+1; bin<=fgkBins; bin++) {
-      fHisto->SetBinContent(bin,0);
-    }
+  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;
@@ -1048,6 +1092,102 @@ Double_t AliQuenchingWeights::CalcQuenchedEnergyKFast(Int_t ipart, Double_t I0,
   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
@@ -1088,16 +1228,16 @@ 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<=4*fLengthMax;l++){
       Char_t hname[100];
-      sprintf(hname,"hDisc-ContQW_%s_%d_%d_%d_%d",meddesc,fInstanceNumber,ipart,medvalue,l);
+      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);
@@ -1132,7 +1272,7 @@ Int_t AliQuenchingWeights::SampleEnergyLoss()
   return 0;
 }
 
-Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t R)
+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)
@@ -1146,28 +1286,30 @@ Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t R)
   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);
+    CalcMult(ipart,r,xxxx,continuous,discrete);
   else
-    CalcSingleHard(ipart,R,xxxx,continuous,discrete);
+    CalcSingleHard(ipart,r,xxxx,continuous,discrete);
 
   if(discrete>=1.) {
     fHisto->SetBinContent(1,1.);
-    for(Int_t bin=2;bin<=fgkBins;bin++) 
+    for(bin=2;bin<=fgkBins;bin++) 
       fHisto->SetBinContent(bin,0.);
     return 0;
   }
 
   fHisto->SetBinContent(bin,continuous);
-  for(Int_t bin=2; bin<=fgkBins; bin++) {
+  for(bin=2; bin<=fgkBins; bin++) {
     xxxx = fHisto->GetBinCenter(bin);
     if(fMultSoft)
-      CalcMult(ipart,R,xxxx,continuous,discrete);
+      CalcMult(ipart,r,xxxx,continuous,discrete);
     else
-      CalcSingleHard(ipart,R,xxxx,continuous,discrete);
+      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);
@@ -1212,14 +1354,14 @@ TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t l
   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,fgkBins,0.,fgkMaxBin*wc);
@@ -1257,14 +1399,14 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t
   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,fgkBins,0.,fgkMaxBin);
@@ -1292,7 +1434,7 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t
   return histx;
 }
 
-TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t R) const 
+TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t r) const 
 {
   // compute P(E) distribution for
   // given ipart = 1 for quark, 2 for gluon 
@@ -1300,13 +1442,13 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t R) const
 
   Char_t meddesc[100];
   if(fMultSoft) {
-    sprintf(meddesc,"MS");
+    snprintf(meddesc, 100, "MS");
   } else {
-    sprintf(meddesc,"SH");
+    snprintf(meddesc, 100, "SH");
   }
 
   Char_t hname[100];
-  sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,R);
+  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)
@@ -1315,7 +1457,7 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t R) const
     histx->SetYTitle("p(#Delta E/#bar#omega_{c})");
   histx->SetLineColor(4);
 
-  Double_t rrrr = R;
+  Double_t rrrr = r;
   Double_t continuous=0.,discrete=0.;
   //loop on histogram channels
   for(Int_t bin=1; bin<=fgkBins; bin++) {
@@ -1362,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);
@@ -1399,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);
@@ -1417,16 +1559,16 @@ TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t medval,TH1F *h
   return h;
 }
 
-TH1F* AliQuenchingWeights::ComputeELossHisto(Int_t ipart,Double_t R) 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);
+  TH1F *dummy = ComputeQWHistoX(ipart,r);
   if(!dummy) return 0;
 
   Char_t hname[100];
-  sprintf(hname,"hELossHistox_%d_%.2f",ipart,R);
+  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;
@@ -1462,19 +1604,19 @@ Double_t AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t medval,TH1F *hEl
   return ret;
 }
 
-Double_t  AliQuenchingWeights::GetMeanELoss(Int_t ipart,Double_t R) const 
+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);
+  TH1F *dummy = ComputeELossHisto(ipart,r);
   if(!dummy) return 0;
   Double_t ret=dummy->GetMean();
   delete dummy;
   return ret;
 }
 
-void AliQuenchingWeights::PlotDiscreteWeights(Double_t len) const
+void AliQuenchingWeights::PlotDiscreteWeights(Double_t len,Double_t qm) const
 {
   // plot discrete weights for given length
 
@@ -1485,59 +1627,67 @@ void AliQuenchingWeights::PlotDiscreteWeights(Double_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);
       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);
       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,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,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 = %.1f 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();
@@ -1553,20 +1703,20 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
   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);
@@ -1590,21 +1740,21 @@ void AliQuenchingWeights::PlotContWeights(Int_t itype,Double_t ell) const
   l1a->SetFillStyle(0);
   l1a->SetBorderSize(0);
   Char_t label[100];
-  sprintf(label,"L = %.1f 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();
@@ -1621,18 +1771,18 @@ void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
   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();
@@ -1656,9 +1806,9 @@ void AliQuenchingWeights::PlotContWeightsVsL(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");
@@ -1669,7 +1819,7 @@ void AliQuenchingWeights::PlotContWeightsVsL(Int_t itype,Double_t medval) const
   c->Update();
 }
 
-void AliQuenchingWeights::PlotAvgELoss(Double_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 
@@ -1682,16 +1832,16 @@ void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t e)  const
   Char_t title[1024];
   Char_t name[1024];
   if(fMultSoft){ 
-    sprintf(title,"Average Energy Loss for Multiple Scattering");
-    sprintf(name,"cavgelossms");
+    snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
+    snprintf(name, 1024, "cavgelossms");
   } else {
-    sprintf(title,"Average Energy 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]");
@@ -1702,28 +1852,30 @@ void AliQuenchingWeights::PlotAvgELoss(Double_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");
 
-  TGraph *gratio=new TGraph(20);
-  for(Int_t i=0;i<20;i++){
+  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);
@@ -1733,14 +1885,14 @@ void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t e)  const
   }
   gratio->SetLineStyle(4);
   gratio->Draw();
-  TLegend *l1a = new TLegend(0.5,0.6,.95,0.8);
+  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 = %.1f 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();
 
@@ -1760,11 +1912,11 @@ void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
   Char_t title[1024];
   Char_t name[1024];
   if(fMultSoft){ 
-    sprintf(title,"Average Energy Loss for Multiple Scattering");
-    sprintf(name,"cavgelossms2");
+    snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
+    snprintf(name, 1024, "cavgelossms2");
   } else {
-    sprintf(title,"Average Energy 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);
@@ -1801,7 +1953,7 @@ void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
   gg->Draw("pl");
 
   TGraph *gratio=new TGraph(20);
-  for(Int_t i=0;i<20;i++){
+  for(i=0;i<20;i++){
     Double_t x,y,x2,y2;
     gg->GetPoint(i,x,y);
     gq->GetPoint(i,x2,y2);
@@ -1810,17 +1962,17 @@ void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const
     else gratio->SetPoint(i,x,0);
   }
   gratio->SetLineStyle(4);
-  gratio->Draw();
+  //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->AddEntry(gratio,"gluon/quark/2.25*10","pl");
   l1a->Draw();
 
   c->Update();
@@ -1839,12 +1991,12 @@ void AliQuenchingWeights::PlotAvgELossVsL(Double_t e)  const
   Char_t name[1024];
   Float_t medval;
   if(fMultSoft){ 
-    sprintf(title,"Average Energy Loss for Multiple Scattering");
-    sprintf(name,"cavgelosslms");
+    snprintf(title, 1024, "Average Energy Loss for Multiple Scattering");
+    snprintf(name, 1024,  "cavgelosslms");
     medval=fQTransport;
   } else {
-    sprintf(title,"Average Energy Loss for Single Hard Scattering");
-    sprintf(name,"cavgelosslsh");
+    snprintf(title, 1024, "Average Energy Loss for Single Hard Scattering");
+    snprintf(name, 1024,  "cavgelosslsh");
     medval=fMu;
   }
 
@@ -1879,7 +2031,7 @@ void AliQuenchingWeights::PlotAvgELossVsL(Double_t e)  const
   gg->Draw("pl");
 
   TGraph *gratio=new TGraph((Int_t)fLengthMax*4);
-  for(Int_t i=0;i<=(Int_t)fLengthMax*4;i++){
+  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);
@@ -1895,9 +2047,9 @@ void AliQuenchingWeights::PlotAvgELossVsL(Double_t e)  const
   l1a->SetBorderSize(0);
   Char_t label[100];
   if(fMultSoft) 
-    sprintf(label,"#hat{q} = %.2f GeV^{2}/fm",medval);
+    snprintf(label, 100, "#hat{q} = %.2f GeV^{2}/fm",medval);
   else
-    sprintf(label,"#mu = %.2f GeV",medval);
+    snprintf(label, 100, "#mu = %.2f GeV",medval);
   l1a->AddEntry(gq,label,"");
   l1a->AddEntry(gq,"quark","pl");
   l1a->AddEntry(gg,"gluon","pl");
@@ -1924,11 +2076,11 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
   Char_t title[1024];
   Char_t name[1024];
   if(fMultSoft){
-    sprintf(title,"Relative Energy Loss for Multiple Scattering");
-    sprintf(name,"cavgelossvsptms");
+    snprintf(title, 1024, "Relative Energy Loss for Multiple Scattering");
+    snprintf(name, 1024, "cavgelossvsptms");
   } else {
-    sprintf(title,"Relative Energy 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);
@@ -1965,7 +2117,7 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,Double_t len) const
   l1a->SetFillStyle(0);
   l1a->SetBorderSize(0);
   Char_t label[100];
-  sprintf(label,"L = %.1f fm",len);
+  snprintf(label, 100, "L = %.1f fm",len);
   l1a->AddEntry(gq,label,"");
   l1a->AddEntry(gq,"quark","pl");
   l1a->AddEntry(gg,"gluon","pl");
@@ -1987,11 +2139,11 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
   Char_t title[1024];
   Char_t name[1024];
   if(fMultSoft){
-    sprintf(title,"Relative Energy Loss for Multiple Scattering");
-    sprintf(name,"cavgelossvsptms2");
+    snprintf(title, 1024, "Relative Energy Loss for Multiple Scattering");
+    snprintf(name,  1024, "cavgelossvsptms2");
   } else {
-    sprintf(title,"Relative Energy 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();
@@ -2027,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");
@@ -2038,10 +2190,11 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const
 
 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.5)>0.125) l++;
+  if((len-l*0.25)>0.125) l++;
   return l;
 }