X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=FASTSIM%2FAliQuenchingWeights.cxx;h=9392a996552b3308683e8a7c0cc0fe85f5b64c79;hb=be33cb74f70a0817a64b2490d36619946e1db698;hp=05c006eef41a0c5dcadf9df3f23adc353c827d2a;hpb=facee35a354abc9e40902cc2a8ff9185a93b6164;p=u%2Fmrichter%2FAliRoot.git diff --git a/FASTSIM/AliQuenchingWeights.cxx b/FASTSIM/AliQuenchingWeights.cxx index 05c006eef41..9392a996552 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" @@ -53,38 +53,46 @@ const Double_t AliQuenchingWeights::fgkConvFmToInvGeV = 1./0.197; const Double_t AliQuenchingWeights::fgkRMax = 1.e6; // hist binning -const Int_t AliQuenchingWeights::fgBins = 1300; -const Double_t AliQuenchingWeights::fgMaxBin = 1.3; +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(); - SetK(); - fECMethod=kReweight; //this is to force printout - SetECMethod(); - SetLengthMax(); - fLengthMaxOld=0; - fInstanceNumber=fgCounter++; - Char_t name[100]; - sprintf(name,"hhistoqw_%d",fInstanceNumber); - fHisto = new TH1F(name,"",fgBins,0.,fgMaxBin); - for(Int_t bin=1;bin<=fgBins;bin++) + fHisto = new TH1F(Form("hhistoqw_%d",fInstanceNumber),"",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 @@ -100,8 +108,8 @@ AliQuenchingWeights::AliQuenchingWeights(const AliQuenchingWeights& a) fInstanceNumber=fgCounter++; Char_t name[100]; sprintf(name,"hhistoqw_%d",fInstanceNumber); - fHisto = new TH1F(name,"",fgBins,0.,fgMaxBin); - for(Int_t bin=1;bin<=fgBins;bin++) + 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 @@ -132,12 +140,12 @@ void AliQuenchingWeights::SetECMethod(kECMethod type) { //set energy constraint method - if(fECMethod==type) return; 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) @@ -727,7 +735,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; } @@ -817,15 +825,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 @@ -833,14 +841,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 @@ -867,9 +875,9 @@ Double_t AliQuenchingWeights::GetELossRandom(Int_t ipart, Double_t length, Doubl Int_t ws=0; while(ret>e){ ret=fHistos[ipart-1][l-1]->GetRandom(); - if(++ws==1e5){ - Warning("GetELossRandomK", - "Aborted reweighting; maximum loss assigned after 1e5 trials."); + if(++ws==1e6){ + Warning("GetELossRandom", + "Aborted reweighting; maximum loss assigned after 1e6 trials."); return e; } } @@ -931,18 +939,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){ - cout << ipart << " " << R << endl; + if(SampleEnergyLoss(ipart,r)!=0){ Fatal("GetELossRandomK","Could not sample energy loss"); return -1000.; } @@ -952,9 +959,9 @@ Double_t AliQuenchingWeights::GetELossRandomK(Int_t ipart, Double_t I0, Double_t Int_t ws=0; while(ret>e){ ret=fHisto->GetRandom(); - if(++ws==1e5){ + if(++ws==1e6){ Warning("GetELossRandomK", - "Aborted reweighting; maximum loss assigned after 1e5 trials."); + "Aborted reweighting; maximum loss assigned after 1e6 trials."); return e; } } @@ -976,6 +983,203 @@ Double_t AliQuenchingWeights::CalcQuenchedEnergyK(Int_t ipart, Double_t I0, Doub 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 + } + + 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.; + + 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 @@ -1028,7 +1232,7 @@ Int_t AliQuenchingWeights::SampleEnergyLoss() sprintf(hname,"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,fgBins,0.,fgMaxBin*wc); + 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); @@ -1036,7 +1240,7 @@ Int_t AliQuenchingWeights::SampleEnergyLoss() Double_t rrrr = CalcR(wc,len); Double_t discrete=0.; // loop on histogram channels - for(Int_t bin=1; bin<=fgBins; bin++) { + for(Int_t bin=1; bin<=fgkBins; bin++) { Double_t xxxx = fHistos[ipart-1][l-1]->GetBinCenter(bin)/wc; Double_t continuous; if(fMultSoft) @@ -1047,20 +1251,20 @@ Int_t AliQuenchingWeights::SampleEnergyLoss() } // add discrete part to distribution if(discrete>=1.) - for(Int_t bin=2;bin<=fgBins;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,fgBins); + 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,fgBins); + Double_t hint=fHistos[ipart-1][l-1]->Integral(1,fgkBins); fHistos[ipart-1][l-1]->Scale(1./hint); } } 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) @@ -1076,35 +1280,34 @@ Int_t AliQuenchingWeights::SampleEnergyLoss(Int_t ipart, Double_t R) 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>=1.) { fHisto->SetBinContent(1,1.); - for(Int_t bin=2;bin<=fgBins;bin++) + for(bin=2;bin<=fgkBins;bin++) fHisto->SetBinContent(bin,0.); return 0; } fHisto->SetBinContent(bin,continuous); - for(Int_t bin=2; bin<=fgBins; 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); } - // add discrete part to distribution - Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgBins); + + Double_t val=discrete/(1.-discrete)*fHisto->Integral(1,fgkBins); fHisto->Fill(0.,val); - - Double_t hint=fHisto->Integral(1,fgBins); + Double_t hint=fHisto->Integral(1,fgkBins); if(hint!=0) fHisto->Scale(1./hint); else { - cout << discrete << " " << hint << " " << continuous << endl; + //cout << discrete << " " << hint << " " << continuous << endl; return -1; } return 0; @@ -1152,14 +1355,14 @@ TH1F* AliQuenchingWeights::ComputeQWHisto(Int_t ipart,Double_t medval,Double_t l sprintf(hname,"hContQWHisto_%s_%d_%d_%d",meddesc,ipart, (Int_t)(medval*1000.),(Int_t)length); - TH1F *hist = new TH1F("hist",hname,fgBins,0.,fgMaxBin*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<=fgBins; bin++) { + for(Int_t bin=1; bin<=fgkBins; bin++) { Double_t xxxx = hist->GetBinCenter(bin)/wc; Double_t continuous,discrete; Int_t ret=0; @@ -1197,7 +1400,7 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t sprintf(hname,"hContQWHistox_%s_%d_%d_%d",meddesc,ipart, (Int_t)(medval*1000.),(Int_t)length); - TH1F *histx = new TH1F("histx",hname,fgBins,0.,fgMaxBin); + 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})"); @@ -1207,7 +1410,7 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t medval,Double_t Double_t rrrr = CalcR(wc,length); //loop on histogram channels - for(Int_t bin=1; bin<=fgBins; bin++) { + for(Int_t bin=1; bin<=fgkBins; bin++) { Double_t xxxx = histx->GetBinCenter(bin); Double_t continuous,discrete; Int_t ret=0; @@ -1222,7 +1425,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 @@ -1236,8 +1439,8 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t R) const } Char_t hname[100]; - sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,R); - TH1F *histx = new TH1F("histx",hname,fgBins,0.,fgMaxBin); + sprintf(hname,"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})"); @@ -1245,10 +1448,10 @@ 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<=fgBins; bin++) { + 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); @@ -1262,13 +1465,13 @@ TH1F* AliQuenchingWeights::ComputeQWHistoX(Int_t ipart,Double_t R) const //add discrete part to distribution if(discrete>=1.) - for(Int_t bin=2;bin<=fgBins;bin++) + for(Int_t bin=2;bin<=fgkBins;bin++) histx->SetBinContent(bin,0.); else { - Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgBins); + Double_t val=discrete/(1.-discrete)*histx->Integral(1,fgkBins); histx->Fill(0.,val); } - Double_t hint=histx->Integral(1,fgBins); + Double_t hint=histx->Integral(1,fgkBins); if(hint!=0) histx->Scale(1./hint); return histx; @@ -1347,17 +1550,17 @@ 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); - TH1F *histx = new TH1F("histxr",hname,fgBins,0.,fgMaxBin); + sprintf(hname,"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(); @@ -1392,19 +1595,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 @@ -1415,50 +1618,58 @@ 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); @@ -1466,8 +1677,8 @@ void AliQuenchingWeights::PlotDiscreteWeights(Double_t len) const Char_t label[100]; sprintf(label,"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(); @@ -1599,7 +1810,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 @@ -1621,7 +1832,7 @@ void AliQuenchingWeights::PlotAvgELoss(Double_t len,Double_t e) const 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]"); @@ -1632,28 +1843,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); @@ -1663,14 +1876,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); 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(); @@ -1731,7 +1944,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); @@ -1740,7 +1953,7 @@ 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); @@ -1750,7 +1963,7 @@ void AliQuenchingWeights::PlotAvgELoss(TH1F *hEll,Double_t e) const 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(); @@ -1809,7 +2022,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); @@ -1968,10 +2181,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; }