Added function to get P(Delta E=0)

diff --git a/FASTSIM/AliQuenchingWeights.cxx b/FASTSIM/AliQuenchingWeights.cxx
index 297fec1a41061c9eb9658bab3905f091a2615c16..e747ec7e9e8de0384185d943a186f1994f3c7a7d 100644 (file)
--- a/FASTSIM/AliQuenchingWeights.cxx
+++ b/FASTSIM/AliQuenchingWeights.cxx
@@ -72,7 +72,6 @@ AliQuenchingWeights::AliQuenchingWeights()
   SetQTransport();
   SetK();
   fECMethod=kDefault; 
-  //SetECMethod();
   SetLengthMax();
   fLengthMaxOld=0;
   fInstanceNumber=fgCounter++;
@@ -135,8 +134,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) 
@@ -1010,7 +1010,7 @@ Double_t AliQuenchingWeights::GetELossRandomKFastR(Int_t ipart, Double_t R, Doub
   if(R>=fgkRMax) {
     R=fgkRMax-1;
   }  
-
+  
   Double_t discrete=0.;
   Double_t continuous=0.;
   Int_t bin=1;
@@ -1020,12 +1020,13 @@ Double_t AliQuenchingWeights::GetELossRandomKFastR(Int_t ipart, Double_t R, Doub
   else
     CalcSingleHard(ipart,R,xxxx,continuous,discrete);
 
-  if(discrete>0.999) {
+  if(discrete>=1.0) {
     return 0.; //no energy loss
   }
 
   fHisto->SetBinContent(bin,continuous);
-  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) {
@@ -1042,19 +1043,24 @@ Double_t AliQuenchingWeights::GetELossRandomKFastR(Int_t ipart, Double_t R, Doub
     }
   }
 
-  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->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(Int_t 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;
@@ -1069,6 +1075,98 @@ 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 &prw_cont,
+                                         Int_t ipart,Double_t I0,Double_t I1,Double_t e)
+{
+  p=1.;prw=1.;prw_cont=1.;
+  Double_t R=CalcRk(I0,I1);
+  if(R<=0.){
+    return;
+  }
+  Double_t wc=CalcWCk(I1);
+  if(wc<=0.){
+    return;
+  }
+  GetZeroLossProbR(p,prw,prw_cont,ipart,R,wc,e);
+}
+
+void AliQuenchingWeights::GetZeroLossProbR(Double_t &p,Double_t &prw,Double_t &prw_cont,
+                                          Int_t ipart, Double_t R,Double_t wc,Double_t e)
+{
+  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);
+  prw_cont=fHisto->GetBinContent(1);
+}
+
 Int_t AliQuenchingWeights::SampleEnergyLoss() 
 {
   // Has to be called to fill the histograms

diff --git a/FASTSIM/AliQuenchingWeights.h b/FASTSIM/AliQuenchingWeights.h
index ca8af435104dde64f9c728a2fbca38a8f2bd2f6c..9906205f9da94120929449235f1afbc48890a44e 100644 (file)
--- a/FASTSIM/AliQuenchingWeights.h
+++ b/FASTSIM/AliQuenchingWeights.h
@@ -22,7 +22,7 @@ class TH1F;
 
 class AliQuenchingWeights : public TObject {
  public:
-  enum kECMethod {kDefault=0,kReweight=1};
+  enum kECMethod {kDefault=0,kReweight=1,kReweightCont=2};
 
   AliQuenchingWeights();
   AliQuenchingWeights(const AliQuenchingWeights& a);
@@ -42,6 +42,13 @@ class AliQuenchingWeights : public TObject {
   Double_t GetELossRandomKFastR(Int_t ipart, Double_t R, Double_t wc, Double_t e=1.e10);
   Double_t CalcQuenchedEnergyKFast(Int_t ipart, Double_t I0, Double_t I1, Double_t e);
 
+  Double_t GetDiscreteWeight(Int_t ipart, Double_t I0, Double_t I1);
+  Double_t GetDiscreteWeightR(Int_t ipart, Double_t R);
+  void GetZeroLossProb(Double_t &p,Double_t &prw,Double_t &prw_cont,
+                      Int_t ipart,Double_t I0,Double_t I1,Double_t e=1.e10);
+  void GetZeroLossProbR(Double_t &p,Double_t &prw, Double_t &prw_cont,
+                       Int_t ipart,Double_t R,Double_t wc,Double_t e=1.e10);
+
   //multiple soft scattering approximation
   Int_t InitMult(const Char_t *contall="$(ALICE_ROOT)/FASTSIM/data/cont_mult.all",
                  const Char_t *discall="$(ALICE_ROOT)/FASTSIM/data/disc_mult.all"); 
@@ -147,7 +154,7 @@ class AliQuenchingWeights : public TObject {
   Int_t fInstanceNumber; //instance number of class
 
   Bool_t fMultSoft;     //approximation type
-  Bool_t fECMethod;     //energy constraint method
+  kECMethod fECMethod;     //energy constraint method
   Double_t fQTransport; //transport coefficient [GeV^2/fm]]
   Double_t fMu;         //Debye screening mass
   Double_t fK;          //proportional constant [fm]