adaptation to AliConfig::GetDeafultFolderName

[u/mrichter/AliRoot.git] / EVGEN / AliGenHIJINGpara.cxx

diff --git a/EVGEN/AliGenHIJINGpara.cxx b/EVGEN/AliGenHIJINGpara.cxx
index 76174048c90fd04927dd8e40d36aef37136669f3..1d12581281aa101eac78ce3b3a630a3d2ebfabcc 100644 (file)
--- a/EVGEN/AliGenHIJINGpara.cxx
+++ b/EVGEN/AliGenHIJINGpara.cxx
@@ -41,8 +41,12 @@
 #include <TClonesArray.h>
 #include <TDatabasePDG.h>
 #include <TF1.h>
+#include <TH1.h>
 #include <TParticle.h>
 #include <TPDGCode.h>
+#include <TCanvas.h>
+#include <TVirtualMC.h>
+#include <TROOT.h>
 
 #include "AliConst.h"
 #include "AliDecayer.h"
@@ -68,23 +72,24 @@ static Double_t ptpi(Double_t *px, Double_t *)
   //     POWER LAW FOR PT > 500 MEV
   //     MT SCALING BELOW (T=160 MEV)
   //
-  const Double_t kp0 = 1.3;
-  const Double_t kxn = 8.28;
-  const Double_t kxlim=0.5;
-  const Double_t kt=0.160;
-  const Double_t kxmpi=0.139;
-  const Double_t kb=1.;
+  const Double_t kp0    = 1.3;
+  const Double_t kxn    = 8.28;
+  const Double_t kxlim  = 0.5;
+  const Double_t kt     = 0.160;
+  const Double_t kxmpi  = 0.139;
+  const Double_t kb     = 1.;
   Double_t y, y1, xmpi2, ynorm, a;
-  Double_t x=*px;
+  Double_t x = *px;
   //
-  y1=TMath::Power(kp0/(kp0+kxlim),kxn);
-  xmpi2=kxmpi*kxmpi;
-  ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
-  a=ynorm/y1;
+  y1 = TMath::Power(kp0 / (kp0 + kxlim), kxn);
+  xmpi2 = kxmpi * kxmpi;
+  ynorm = kb * (TMath::Exp(-sqrt(kxlim * kxlim + xmpi2) / kt ));
+  a = ynorm / y1;
   if (x > kxlim)
-    y=a*TMath::Power(kp0/(kp0+x),kxn);
+    y = a * TMath::Power(kp0 / (kp0 + x), kxn);
   else
-    y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
+    y = kb* TMath::Exp(-sqrt(x * x + xmpi2) / kt);
+  
   return y*x;
 }
 
@@ -217,15 +222,21 @@ void AliGenHIJINGpara::Init()
     Float_t etaMax = -TMath::Log(TMath::Tan(
        TMath::Max((Double_t)fThetaMin/2,1.e-10)));
     fPtpi   = new TF1("ptpi",&ptpi,0,20,0);
+    gROOT->GetListOfFunctions()->Remove(fPtpi);
     fPtka   = new TF1("ptka",&ptka,0,20,0);
+    gROOT->GetListOfFunctions()->Remove(fPtka);
+    fPtpi->SetNpx(1000);
+    fPtka->SetNpx(1000);
     fETApic = new TF1("etapic",&etapic,etaMin,etaMax,0);
+    gROOT->GetListOfFunctions()->Remove(fETApic);
     fETAkac = new TF1("etakac",&etakac,etaMin,etaMax,0);
+    gROOT->GetListOfFunctions()->Remove(fETAkac);
 
-    TF1 etaPic0("etapic",&etapic,-7,7,0);
-    TF1 etaKac0("etakac",&etakac,-7,7,0);
+    TF1 etaPic0("etaPic0",&etapic,-7,7,0);
+    TF1 etaKac0("etaKac0",&etakac,-7,7,0);
 
-    TF1 ptPic0("ptpi",&ptpi,0.,15.,0);
-    TF1 ptKac0("ptka",&ptka,0.,15.,0);
+    TF1 ptPic0("ptPic0",&ptpi,0.,15.,0);
+    TF1 ptKac0("ptKac0",&ptka,0.,15.,0);
 
     Float_t intETApi  = etaPic0.Integral(-0.5, 0.5);
     Float_t intETAka  = etaKac0.Integral(-0.5, 0.5);
@@ -245,8 +256,16 @@ void AliGenHIJINGpara::Init()
 //  Fraction of events corresponding to the selected phi-range    
     Float_t phiFrac    = (fPhiMax-fPhiMin)/2/TMath::Pi();
 
+    
     fParentWeight = Float_t(fNpart)/(intETASel*ptFrac*phiFrac);
     
+    if (fAnalog != 0) {
+       fPtWgtPi = (fPtMax - fPtMin) / fPtpi->Integral(0., 20.);
+       fPtWgtKa = (fPtMax - fPtMin) / fPtka->Integral(0., 20.);
+       fParentWeight = Float_t(fNpart)/(intETASel*phiFrac);
+    }
+    
+    
     printf("%s: The number of particles in the selected kinematic region corresponds to %f percent of a full event\n ", 
           ClassName(),100.*fParentWeight);
 
@@ -281,7 +300,7 @@ void AliGenHIJINGpara::Generate()
     const Int_t kKaons[4] = {kK0Long, kK0Short, kKPlus, kKMinus};
     //
     Float_t origin[3];
-    Float_t pt, pl, ptot;
+    Float_t pt, pl, ptot, wgt;
     Float_t phi, theta;
     Float_t p[3];
     Int_t i, part, j;
@@ -302,23 +321,32 @@ void AliGenHIJINGpara::Generate()
     eventVertex[0] = origin[0];
     eventVertex[1] = origin[1];
     eventVertex[2] = origin[2];
-
+     
     for(i=0;i<fNpart;i++) {
        while(1) {
-           Rndm(random,3);
+           Rndm(random,4);
            if(random[0]<kBorne) {
                part=kPions[Int_t (random[1]*3)];
                ptf=fPtpi;
                etaf=fETApic;
+               wgt = fPtWgtPi;
            } else {
                part=kKaons[Int_t (random[1]*4)];
                ptf=fPtka;
                etaf=fETAkac;
+               wgt = fPtWgtKa;
            }
            phi=fPhiMin+random[2]*(fPhiMax-fPhiMin);
            theta=2*TMath::ATan(TMath::Exp(-etaf->GetRandom()));
            if(theta<fThetaMin || theta>fThetaMax) continue;
-           pt=ptf->GetRandom();
+        
+           if (fAnalog == 0) { 
+               pt = ptf->GetRandom();
+           } else {
+               pt = fPtMin + random[3] * (fPtMax - fPtMin);
+           }
+           
+           
            pl=pt/TMath::Tan(theta);
            ptot=TMath::Sqrt(pt*pt+pl*pl);
            if(ptot<fPMin || ptot>fPMax) continue;
@@ -332,6 +360,14 @@ void AliGenHIJINGpara::Generate()
                        TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
                }
            }
+           
+           if (fAnalog == 0) { 
+               wgt = fParentWeight;
+           } else {
+               wgt *= (fParentWeight * ptf->Eval(pt));
+           }
+           
+           
            if (part == kPi0 && fPi0Decays){
 //
 //          Decay pi0 if requested
@@ -396,7 +432,24 @@ void AliGenHIJINGpara::DecayPi0(Float_t* orig, Float_t * p)
     fNt = nt;
 }
 
-void AliGenHIJINGpara::Copy(AliGenHIJINGpara &) const
+void AliGenHIJINGpara::Copy(TObject &) const
 {
   Fatal("Copy","Not implemented!\n");
 }
+
+
+void AliGenHIJINGpara::Draw( const char * /*opt*/)
+{
+    //
+    // Draw the pT and y Distributions
+    //
+     TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
+     c0->Divide(2,1);
+     c0->cd(1);
+     fPtpi->Draw();
+     fPtpi->GetHistogram()->SetXTitle("p_{T} (GeV)");     
+     c0->cd(2);
+     fPtka->Draw();
+     fPtka->GetHistogram()->SetXTitle("p_{T} (GeV)");     
+
+}