X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliGenHIJINGpara.cxx;h=2a462bdf88c1c0da60e572e24d4e3194bc8fe39f;hb=2ad38d566a92536912bbf772ba5d7732cb2e95a2;hp=17fae20c87a23625da0f2e929e0f270a4420ee1b;hpb=edbba5db3169dd442b975bb9824947fc573ff7c3;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenHIJINGpara.cxx b/EVGEN/AliGenHIJINGpara.cxx index 17fae20c87a..2a462bdf88c 100644 --- a/EVGEN/AliGenHIJINGpara.cxx +++ b/EVGEN/AliGenHIJINGpara.cxx @@ -13,75 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.19 2003/01/14 10:50:18 alibrary -Cleanup of STEER coding conventions - -Revision 1.18 2002/12/11 11:58:11 morsch -Bug in formula for pi0 energy for decay corrected. - -Revision 1.17 2002/12/10 17:44:57 morsch -Correct mother child relation for pi0. - -Revision 1.16 2002/11/28 11:46:15 morsch -Don't track pi0 if already decayed. - -Revision 1.15 2002/11/28 11:38:53 morsch -Typo corrected. - -Revision 1.14 2002/11/26 17:12:36 morsch -Decay pi0 if requested. - -Revision 1.13 2002/10/14 14:55:35 hristov -Merging the VirtualMC branch to the main development branch (HEAD) - -Revision 1.11.4.1 2002/07/24 08:56:28 alibrary -Updating EVGEN on TVirtulaMC - -Revision 1.12 2002/06/19 06:56:34 hristov -Memory leak corrected - -Revision 1.11 2002/03/20 10:21:13 hristov -Set fPtMax to 15 GeV in order to avoid some numerical problems - -Revision 1.10 2001/10/15 16:44:46 morsch -- Possibility for vertex distribution truncation. -- Write mc header with vertex position. - -Revision 1.9 2001/07/27 17:09:36 morsch -Use local SetTrack, KeepTrack and SetHighWaterMark methods -to delegate either to local stack or to stack owned by AliRun. -(Piotr Skowronski, A.M.) - -Revision 1.8 2001/07/20 11:03:58 morsch -Issue warning message if used outside allowed eta range (-8 to 8). - -Revision 1.7 2001/07/17 12:41:01 morsch -- Calculation of fraction of event corresponding to selected pt-range corrected -(R. Turrisi) -- Parent weight corrected. - -Revision 1.6 2001/05/16 14:57:10 alibrary -New files for folders and Stack - -Revision 1.5 2000/12/21 16:24:06 morsch -Coding convention clean-up - -Revision 1.4 2000/11/30 07:12:50 alibrary -Introducing new Rndm and QA classes - -Revision 1.3 2000/10/02 21:28:06 fca -Removal of useless dependecies via forward declarations - -Revision 1.2 2000/07/11 18:24:55 fca -Coding convention corrections + few minor bug fixes - -Revision 1.1 2000/06/09 20:20:30 morsch -Same class as previously in AliSimpleGen.cxx -All coding rule violations except RS3 corrected (AM) - -*/ +/* $Id$ */ // Parameterisation of pi and K, eta and pt distributions // used for the ALICE TDRs. @@ -109,21 +41,26 @@ All coding rule violations except RS3 corrected (AM) #include #include #include +#include #include #include +#include +#include #include "AliConst.h" #include "AliDecayer.h" -#include "AliDecayerPythia.h" #include "AliGenEventHeader.h" #include "AliGenHIJINGpara.h" #include "AliRun.h" ClassImp(AliGenHIJINGpara) -AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para) + +AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para): + AliGenerator(para) { -// copy constructor +// Copy constructor + para.Copy(*this); } //_____________________________________________________________________________ @@ -134,23 +71,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; } @@ -284,14 +222,16 @@ void AliGenHIJINGpara::Init() TMath::Max((Double_t)fThetaMin/2,1.e-10))); fPtpi = new TF1("ptpi",&ptpi,0,20,0); fPtka = new TF1("ptka",&ptka,0,20,0); + fPtpi->SetNpx(1000); + fPtka->SetNpx(1000); fETApic = new TF1("etapic",&etapic,etaMin,etaMax,0); fETAkac = new TF1("etakac",&etakac,etaMin,etaMax,0); - 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); @@ -311,8 +251,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); @@ -326,7 +274,8 @@ void AliGenHIJINGpara::Init() } // // - if (fPi0Decays) fDecayer = new AliDecayerPythia(); + if (fPi0Decays && gMC) + fDecayer = gMC->GetDecayer(); } @@ -346,7 +295,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; @@ -359,39 +308,40 @@ void AliGenHIJINGpara::Generate() for (j=0;j<3;j++) origin[j]=fOrigin[j]; if(fVertexSmear == kPerEvent) { - Float_t dv[3]; - dv[2] = 1.e10; - while(TMath::Abs(dv[2]) > fCutVertexZ*fOsigma[2]) { - Rndm(random,6); - for (j=0; j < 3; j++) { - dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())* - TMath::Sqrt(-2*TMath::Log(random[2*j+1])); - } - } - for (j=0; j < 3; j++) origin[j] += dv[j]; + Vertex(); + for (j=0; j < 3; j++) origin[j] = fVertex[j]; } // if kPerEvent TArrayF eventVertex; eventVertex.Set(3); eventVertex[0] = origin[0]; eventVertex[1] = origin[1]; eventVertex[2] = origin[2]; - + for(i=0;iGetRandom())); if(thetafThetaMax) 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(ptotfPMax) continue; @@ -405,14 +355,22 @@ 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 - SetTrack(0,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight); + PushTrack(0,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight); KeepTrack(fNt); DecayPi0(origin, p); } else { - SetTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight); + PushTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,fNt,fParentWeight); KeepTrack(fNt); } @@ -429,12 +387,6 @@ void AliGenHIJINGpara::Generate() gAlice->SetGenEventHeader(header); } -AliGenHIJINGpara& AliGenHIJINGpara::operator=(const AliGenHIJINGpara& rhs) -{ -// Assignment operator - return *this; -} - void AliGenHIJINGpara::SetPtRange(Float_t ptmin, Float_t ptmax) { AliGenerator::SetPtRange(ptmin, ptmax); } @@ -469,8 +421,30 @@ void AliGenHIJINGpara::DecayPi0(Float_t* orig, Float_t * p) p[2] = iParticle->Pz(); Int_t part = iParticle->GetPdgCode(); - SetTrack(fTrackIt, fNt, part, p, orig, polar, 0, kPDecay, nt, fParentWeight); + PushTrack(fTrackIt, fNt, part, p, orig, polar, 0, kPDecay, nt, fParentWeight); KeepTrack(nt); } fNt = nt; } + +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)"); + +}