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;
SetMu();
SetQTransport();
SetK();
- fECMethod=kReweight; //this is to force printout
- SetECMethod();
+ fECMethod=kDefault;
+ //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(name,"",fgkBins,0.,fgkMaxBin);
+ for(Int_t bin=1;bin<=fgkBins;bin++)
fHisto->SetBinContent(bin,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(name,"",fgkBins,0.,fgkMaxBin);
+ for(Int_t bin=1;bin<=fgkBins;bin++)
fHisto->SetBinContent(bin,0.);
//Missing in the class is the pathname
{
//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.");
return -1000.;
}
if(SampleEnergyLoss(ipart,R)!=0){
- cout << ipart << " " << R << endl;
Fatal("GetELossRandomK","Could not sample energy loss");
return -1000.;
}
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.;
+ }
+
+ 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>0.999) {
+ return 0.0;
+ }
+
+ //set first bin
+ fHisto->SetBinContent(bin,continuous);
+
+ if(fMultSoft) {
+ for(Int_t bin=2; bin<=fgkBins; bin++) {
+ xxxx = fHisto->GetBinCenter(bin);
+ CalcMult(ipart,R,xxxx,continuous,discrete);
+ fHisto->SetBinContent(bin,continuous);
+ }
+ } else {
+ for(Int_t bin=2; bin<=fgkBins; bin++) {
+ xxxx = fHisto->GetBinCenter(bin);
+ 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);
+ }
+ }
+
+ 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;
+}
+
Int_t AliQuenchingWeights::SampleEnergyLoss()
{
// Has to be called to fill the histograms
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);
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)
}
// 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);
}
}
if(discrete>=1.) {
fHisto->SetBinContent(1,1.);
- for(Int_t bin=2;bin<=fgBins;bin++)
+ for(Int_t bin=2;bin<=fgkBins;bin++)
fHisto->SetBinContent(bin,0.);
return 0;
}
fHisto->SetBinContent(bin,continuous);
- for(Int_t bin=2; bin<=fgBins; bin++) {
+ for(Int_t bin=2; bin<=fgkBins; bin++) {
xxxx = fHisto->GetBinCenter(bin);
if(fMultSoft)
CalcMult(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;
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;
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})");
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;
Char_t hname[100];
sprintf(hname,"hQWHistox_%s_%d_%.2f",meddesc,ipart,R);
- 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})");
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);
//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;
Char_t hname[100];
sprintf(hname,"hELossHistox_%d_%.2f",ipart,R);
- TH1F *histx = new TH1F("histxr",hname,fgBins,0.,fgMaxBin);
+ 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();