1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //====================================================================================================================================================
18 // Parametric generator of primary pions and kaons
20 // Contact author: antonio.uras@cern.ch
22 //====================================================================================================================================================
28 #include "AliGenEventHeader.h"
29 #include "TDatabasePDG.h"
33 #include "AliGenParamPionsKaons.h"
35 ClassImp(AliGenParamPionsKaons)
37 //====================================================================================================================================================
39 AliGenParamPionsKaons::AliGenParamPionsKaons():
43 fPtVsRapidityPrimaryPosPions(0x0),
44 fPtVsRapidityPrimaryNegPions(0x0),
45 fPtVsRapidityPrimaryPosKaons(0x0),
46 fPtVsRapidityPrimaryNegKaons(0x0),
49 // Default constructor
53 //====================================================================================================================================================
55 AliGenParamPionsKaons::AliGenParamPionsKaons(Int_t nPart, Char_t *inputFile):
59 fPtVsRapidityPrimaryPosPions(0x0),
60 fPtVsRapidityPrimaryNegPions(0x0),
61 fPtVsRapidityPrimaryPosKaons(0x0),
62 fPtVsRapidityPrimaryNegKaons(0x0),
65 // Standard constructor
67 fName = "ParamPionsKaons";
68 fTitle = "Parametric pions and kaons generator";
70 LoadInputHistos(inputFile);
74 //====================================================================================================================================================
76 void AliGenParamPionsKaons::Generate() {
78 // Generate one trigger
80 Double_t polar[3]= {0,0,0};
85 Double_t mass=0., pt=0., rap=0., mom=0., energy=0, mt=0., phi=0., time=0.;
89 for (Int_t j=0; j<3; j++) origin[j] = fOrigin[j];
91 if (fVertexSmear==kPerEvent) {
93 for (Int_t j=0; j<3; j++) origin[j] = fVertex[j];
97 Int_t nPartGenerated = 0;
99 while (nPartGenerated<fNpart) {
101 pdgCode = TMath::Nint(fHistPdgCode->GetRandom());
102 if (TMath::Abs(pdgCode)==321 && !fGenerateKaon) continue;
103 if (TMath::Abs(pdgCode)==211 && !fGeneratePion) continue;
105 mass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
109 fPtVsRapidityPrimaryPosPions->GetRandom2(pt, rap);
112 fPtVsRapidityPrimaryNegPions->GetRandom2(pt, rap);
115 fPtVsRapidityPrimaryPosKaons->GetRandom2(pt, rap);
118 fPtVsRapidityPrimaryNegKaons->GetRandom2(pt, rap);
122 mt = TMath::Sqrt(mass*mass + pt*pt);
123 energy = mt * TMath::CosH(rap);
124 mom = TMath::Sqrt(energy*energy - mass*mass);
126 if (TestBit(kYRange)) if (rap<fYMin || rap>fYMax) continue;
127 if (TestBit(kMomentumRange)) if (mom<fPMin || rap>fPMax) continue;
128 if (TestBit(kPtRange)) if (pt<fPtMin || pt>fPtMax) continue;
130 phi = fPhiMin + gRandom->Rndm()*(fPhiMax-fPhiMin);
131 p[0] = pt*TMath::Cos(phi);
132 p[1] = pt*TMath::Sin(phi);
133 p[2] = mt*TMath::SinH(rap);
135 PushTrack(1, -1, Int_t(pdgCode),
136 p[0],p[1],p[2],energy,
137 origin[0],origin[1],origin[2],Double_t(time),
138 polar[0],polar[1],polar[2],
139 kPPrimary, nt, 1., 1);
145 AliGenEventHeader* header = new AliGenEventHeader("ParamPionsKaons");
146 header->SetPrimaryVertex(fVertex);
147 header->SetNProduced(fNpart);
148 header->SetInteractionTime(fTime);
150 // Passes header either to the container or to gAlice
152 fContainer->AddHeader(header);
155 gAlice->SetGenEventHeader(header);
160 //====================================================================================================================================================
162 void AliGenParamPionsKaons::Init() {
164 // Initialisation, check consistency of selected ranges
165 if (TestBit(kPtRange) && TestBit(kMomentumRange))
166 Fatal("Init","You should not set the momentum range and the pt range at the same time!\n");
167 if ((!TestBit(kPtRange)) && (!TestBit(kMomentumRange)))
168 Fatal("Init","You should set either the momentum or the pt range!\n");
169 if ((TestBit(kYRange) && TestBit(kThetaRange)) || (TestBit(kYRange) && TestBit(kEtaRange)) || (TestBit(kEtaRange) && TestBit(kThetaRange)))
170 Fatal("Init","You should only set the range of one of these variables: y, eta or theta\n");
171 if ((!TestBit(kYRange)) && (!TestBit(kEtaRange)) && (!TestBit(kThetaRange)))
172 Fatal("Init","You should set the range of one of these variables: y, eta or theta\n");
174 AliPDG::AddParticlesToPdgDataBase();
178 //====================================================================================================================================================
180 void AliGenParamPionsKaons::LoadInputHistos(Char_t *inputFile) {
182 TFile *fileIn = new TFile(inputFile);
184 TH2D *myPtVsRapidityPrimaryPosPions = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryPosPions");
185 TH2D *myPtVsRapidityPrimaryNegPions = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryNegPions");
186 TH2D *myPtVsRapidityPrimaryPosKaons = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryPosKaons");
187 TH2D *myPtVsRapidityPrimaryNegKaons = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryNegKaons");
188 TH1D *myHistPdgCode = (TH1D*) fileIn->Get("fHistPdgCode");
190 myPtVsRapidityPrimaryPosPions -> SetName("myPtVsRapidityPrimaryPosPions");
191 myPtVsRapidityPrimaryNegPions -> SetName("myPtVsRapidityPrimaryNegPions");
192 myPtVsRapidityPrimaryPosKaons -> SetName("myPtVsRapidityPrimaryPosKaons");
193 myPtVsRapidityPrimaryNegKaons -> SetName("myPtVsRapidityPrimaryNegKaons");
194 myHistPdgCode -> SetName("myHistPdgCode");
196 fPtVsRapidityPrimaryPosPions = new TH2D("fPtVsRapidityPrimaryPosPions","",
197 myPtVsRapidityPrimaryPosPions->GetXaxis()->GetNbins(),
198 myPtVsRapidityPrimaryPosPions->GetXaxis()->GetXmin(),
199 myPtVsRapidityPrimaryPosPions->GetXaxis()->GetXmax(),
200 myPtVsRapidityPrimaryPosPions->GetYaxis()->GetNbins(),
201 myPtVsRapidityPrimaryPosPions->GetYaxis()->GetXmin(),
202 myPtVsRapidityPrimaryPosPions->GetYaxis()->GetXmax());
203 for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryPosPions->GetXaxis()->GetNbins(); iBinX++) {
204 for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryPosPions->GetYaxis()->GetNbins(); iBinY++) {
205 fPtVsRapidityPrimaryPosPions->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryPosPions->GetBinContent(iBinX+1,iBinY+1));
209 fPtVsRapidityPrimaryNegPions = new TH2D("fPtVsRapidityPrimaryNegPions","",
210 myPtVsRapidityPrimaryNegPions->GetXaxis()->GetNbins(),
211 myPtVsRapidityPrimaryNegPions->GetXaxis()->GetXmin(),
212 myPtVsRapidityPrimaryNegPions->GetXaxis()->GetXmax(),
213 myPtVsRapidityPrimaryNegPions->GetYaxis()->GetNbins(),
214 myPtVsRapidityPrimaryNegPions->GetYaxis()->GetXmin(),
215 myPtVsRapidityPrimaryNegPions->GetYaxis()->GetXmax());
216 for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryNegPions->GetXaxis()->GetNbins(); iBinX++) {
217 for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryNegPions->GetYaxis()->GetNbins(); iBinY++) {
218 fPtVsRapidityPrimaryNegPions->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryNegPions->GetBinContent(iBinX+1,iBinY+1));
222 fPtVsRapidityPrimaryPosKaons = new TH2D("fPtVsRapidityPrimaryPosKaons","",
223 myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetNbins(),
224 myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetXmin(),
225 myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetXmax(),
226 myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetNbins(),
227 myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetXmin(),
228 myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetXmax());
229 for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetNbins(); iBinX++) {
230 for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetNbins(); iBinY++) {
231 fPtVsRapidityPrimaryPosKaons->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryPosKaons->GetBinContent(iBinX+1,iBinY+1));
235 fPtVsRapidityPrimaryNegKaons = new TH2D("fPtVsRapidityPrimaryNegKaons","",
236 myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetNbins(),
237 myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetXmin(),
238 myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetXmax(),
239 myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetNbins(),
240 myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetXmin(),
241 myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetXmax());
242 for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetNbins(); iBinX++) {
243 for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetNbins(); iBinY++) {
244 fPtVsRapidityPrimaryNegKaons->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryNegKaons->GetBinContent(iBinX+1,iBinY+1));
248 fHistPdgCode = new TH1D("fHistPdgCode","",
249 myHistPdgCode->GetXaxis()->GetNbins(),
250 myHistPdgCode->GetXaxis()->GetXmin(),
251 myHistPdgCode->GetXaxis()->GetXmax());
252 for (Int_t iBinX=0; iBinX<myHistPdgCode->GetXaxis()->GetNbins(); iBinX++) {
253 fHistPdgCode->SetBinContent(iBinX+1, myHistPdgCode->GetBinContent(iBinX+1));
256 // fHistPdgCode = new TH1D("fHistPdgCode", "fHistPdgCode", 10, 0., 10);
257 // fHistPdgCode -> Fill(3.);
261 //====================================================================================================================================================