X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliGenHIJINGpara.cxx;h=3596b4deb8b27d9ca3a4b1b52b6965ea4df7dfa5;hb=418e091f31ef997ccabfd6dfc4448dfc6e495c11;hp=dc4b01a96756262839c6616c9d1a86cfac48ad1d;hpb=0af12c00a39ec44e138214fa2c156f236a68f52a;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenHIJINGpara.cxx b/EVGEN/AliGenHIJINGpara.cxx index dc4b01a9675..3596b4deb8b 100644 --- a/EVGEN/AliGenHIJINGpara.cxx +++ b/EVGEN/AliGenHIJINGpara.cxx @@ -38,28 +38,26 @@ /////////////////////////////////////////////////////////////////// #include +#include #include #include #include -#include +#include #include +#include +#include #include #include "AliConst.h" #include "AliDecayer.h" #include "AliGenEventHeader.h" #include "AliGenHIJINGpara.h" +#include "AliLog.h" #include "AliRun.h" ClassImp(AliGenHIJINGpara) -AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para): - AliGenerator(para) -{ -// Copy constructor - para.Copy(*this); -} //_____________________________________________________________________________ static Double_t ptpi(Double_t *px, Double_t *) @@ -69,23 +67,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; } @@ -159,40 +158,46 @@ static Double_t etakac( Double_t *py, Double_t *) //_____________________________________________________________________________ AliGenHIJINGpara::AliGenHIJINGpara() - :AliGenerator() + :AliGenerator(), + fNt(-1), + fNpartProd(0), + fPi0Decays(kFALSE), + fPtWgtPi(0.), + fPtWgtKa(0.), + fPtpi(0), + fPtka(0), + fETApic(0), + fETAkac(0), + fDecayer(0) { // // Default constructor // - fPtpi = 0; - fPtka = 0; - fETApic = 0; - fETAkac = 0; - fDecayer = 0; - fNt = -1; SetCutVertexZ(); SetPtRange(); - SetPi0Decays(); } //_____________________________________________________________________________ AliGenHIJINGpara::AliGenHIJINGpara(Int_t npart) - :AliGenerator(npart) + :AliGenerator(npart), + fNt(-1), + fNpartProd(npart), + fPi0Decays(kFALSE), + fPtWgtPi(0.), + fPtWgtKa(0.), + fPtpi(0), + fPtka(0), + fETApic(0), + fETAkac(0), + fDecayer(0) { // // Standard constructor // fName="HIJINGpara"; fTitle="HIJING Parametrisation Particle Generator"; - fPtpi = 0; - fPtka = 0; - fETApic = 0; - fETAkac = 0; - fDecayer = 0; - fNt = -1; SetCutVertexZ(); SetPtRange(); - SetPi0Decays(); } //_____________________________________________________________________________ @@ -218,15 +223,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); @@ -247,30 +258,36 @@ void AliGenHIJINGpara::Init() Float_t phiFrac = (fPhiMax-fPhiMin)/2/TMath::Pi(); - fParentWeight = Float_t(fNpart)/(intETASel*ptFrac*phiFrac); + fParentWeight = (intETASel*ptFrac*phiFrac) / Float_t(fNpart); if (fAnalog != 0) { fPtWgtPi = (fPtMax - fPtMin) / fPtpi->Integral(0., 20.); fPtWgtKa = (fPtMax - fPtMin) / fPtka->Integral(0., 20.); - fParentWeight = Float_t(fNpart)/(intETASel*phiFrac); + fParentWeight = (intETASel*phiFrac) / Float_t(fNpart); } - printf("%s: The number of particles in the selected kinematic region corresponds to %f percent of a full event\n ", - ClassName(),100.*fParentWeight); + AliInfo(Form("The number of particles in the selected kinematic region corresponds to %f percent of a full event", + 100./ fParentWeight)); // Issue warning message if etaMin or etaMax are outside the alowed range // of the parametrization if (etaMin < -8.001 || etaMax > 8.001) { - printf("\n \n WARNING FROM AliGenHIJINGPara !"); - printf("\n YOU ARE USING THE PARAMETERISATION OUTSIDE "); - printf("\n THE ALLOWED PSEUDORAPIDITY RANGE (-8. - 8.)"); - printf("\n YOUR LIMITS: %f %f \n \n ", etaMin, etaMax); + AliWarning("\nYOU ARE USING THE PARAMETERISATION OUTSIDE "); + AliWarning("THE ALLOWED PSEUDORAPIDITY RANGE (-8. - 8.)"); + AliWarning(Form("YOUR LIMITS: %f %f \n ", etaMin, etaMax)); } // // if (fPi0Decays && gMC) fDecayer = gMC->GetDecayer(); + + if (fPi0Decays) + { + fDecayer->SetForceDecay(kNeutralPion); + fDecayer->Init(); + } + } @@ -350,7 +367,7 @@ void AliGenHIJINGpara::Generate() TMath::Sqrt(-2*TMath::Log(random[2*j+1])); } } - + if (fAnalog == 0) { wgt = fParentWeight; } else { @@ -361,11 +378,13 @@ void AliGenHIJINGpara::Generate() if (part == kPi0 && fPi0Decays){ // // Decay pi0 if requested - PushTrack(0,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight); + PushTrack(0,-1,part,p,origin,polar,0,kPPrimary,fNt,wgt); KeepTrack(fNt); DecayPi0(origin, p); } else { - PushTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight); + // printf("fNt %d", fNt); + PushTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,fNt,wgt); + KeepTrack(fNt); } @@ -379,6 +398,7 @@ void AliGenHIJINGpara::Generate() AliGenEventHeader* header = new AliGenEventHeader("HIJINGparam"); // Event Vertex header->SetPrimaryVertex(eventVertex); + header->SetNProduced(fNpartProd); gAlice->SetGenEventHeader(header); } @@ -407,7 +427,8 @@ void AliGenHIJINGpara::DecayPi0(Float_t* orig, Float_t * p) // Float_t polar[3] = {0., 0., 0.}; Int_t np = fDecayer->ImportParticles(particles); - Int_t nt; + fNpartProd += (np-1); + Int_t nt = 0; for (Int_t i = 1; i < np; i++) { TParticle* iParticle = (TParticle *) particles->At(i); @@ -422,7 +443,18 @@ void AliGenHIJINGpara::DecayPi0(Float_t* orig, Float_t * p) fNt = nt; } -void AliGenHIJINGpara::Copy(AliGenHIJINGpara &) const +void AliGenHIJINGpara::Draw( const char * /*opt*/) { - Fatal("Copy","Not implemented!\n"); + // + // 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)"); + }