X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=FASTSIM%2FAliQuenchingWeights.cxx;h=b8e0c2a5dfce1d60a739a372b9bdab16823612fb;hb=15a2657d65de145ed24fe7662470198336f1d1d7;hp=4df0597dcdb2f58ef6eb33521bec371c0e7a68c4;hpb=8ab8044efb1a89fc60c431cb3a86e45304ecbc4f;p=u%2Fmrichter%2FAliRoot.git diff --git a/FASTSIM/AliQuenchingWeights.cxx b/FASTSIM/AliQuenchingWeights.cxx index 4df0597dcdb..b8e0c2a5dfc 100644 --- a/FASTSIM/AliQuenchingWeights.cxx +++ b/FASTSIM/AliQuenchingWeights.cxx @@ -15,13 +15,12 @@ /* $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 @@ -41,6 +40,7 @@ #include #include #include +#include #include #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,86 +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 } - + if (!fHisto) Init(); + fHisto->SetBinContent(bin,continuous); -#if 1 /* fast new version*/ - const Int_t kbinmax=fHisto->FindBin(e/wc); + 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<=kbinmax; 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<=kbinmax; 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); } } - const Double_t kdelta=fHisto->Integral(1,kbinmax); - - Double_t val=discrete*fgkBins/fgkMaxBin; - fHisto->Fill(0.,val); - - if(fECMethod==kReweight) + if(fECMethod==kReweight){ fHisto->SetBinContent(kbinmax+1,0); - else - fHisto->SetBinContent(kbinmax+1,(1-discrete)*fgkBins/fgkMaxBin-kdelta); - - for(Int_t bin=kbinmax+2; bin<=fgkBins; bin++) { - fHisto->SetBinContent(bin,0); - } - -#else - if(fMultSoft) { - for(Int_t bin=2; bin<=fgkBins; bin++) { - xxxx = fHisto->GetBinCenter(bin); - CalcMult(ipart,R,xxxx,continuous,discrete); - fHisto->SetBinContent(bin,continuous); - } + 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 { - for(Int_t bin=2; bin<=fgkBins; bin++) { - xxxx = fHisto->GetBinCenter(bin); - CalcSingleHard(ipart,R,xxxx,continuous,discrete); - fHisto->SetBinContent(bin,continuous); - } + 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); } - - // 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); - } + for(bin=kbinmax+2; bin<=fgkBins; bin++) { + fHisto->SetBinContent(bin,0); } -#endif - + //cout << kbinmax << " " << discrete << " " << fHisto->Integral() << endl; Double_t ret=fHisto->GetRandom()*wc; if(ret>e) return e; return ret; @@ -1078,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 @@ -1118,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); @@ -1162,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) @@ -1176,26 +1286,27 @@ 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); } @@ -1243,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); @@ -1288,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); @@ -1323,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 @@ -1331,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) @@ -1346,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++) { @@ -1393,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); @@ -1430,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); @@ -1448,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; @@ -1493,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 @@ -1516,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(); @@ -1584,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); @@ -1621,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(); @@ -1652,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(); @@ -1687,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"); @@ -1700,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 @@ -1713,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]"); @@ -1733,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;iGetPoint(i,x,y); gq->GetPoint(i,x2,y2); @@ -1764,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(); @@ -1791,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); @@ -1832,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); @@ -1841,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," = %.2f fm",hEll->GetMean()); + snprintf(label, 100, " = %.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(); @@ -1870,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; } @@ -1910,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); @@ -1926,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"); @@ -1955,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); @@ -1996,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"); @@ -2018,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(); @@ -2058,7 +2179,7 @@ void AliQuenchingWeights::PlotAvgELossVsPt(Double_t medval,TH1F *hEll) const l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; - sprintf(label," = %.2f fm",hEll->GetMean()); + snprintf(label, 100, " = %.2f fm",hEll->GetMean()); l1a->AddEntry(gq,label,""); l1a->AddEntry(gq,"quark","pl"); l1a->AddEntry(gg,"gluon","pl"); @@ -2069,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; }