Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
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 : **************************************************************************/
15 :
16 : //====================================================================================================================================================
17 : //
18 : // Parametric generator of primary pions and kaons
19 : //
20 : // Contact author: antonio.uras@cern.ch
21 : //
22 : //====================================================================================================================================================
23 :
24 : #include "TPDGCode.h"
25 :
26 : #include "AliConst.h"
27 : #include "AliRun.h"
28 : #include "AliGenEventHeader.h"
29 : #include "TDatabasePDG.h"
30 : #include "AliPDG.h"
31 : #include "TFile.h"
32 : #include "TROOT.h"
33 : #include "AliGenParamPionsKaons.h"
34 : #include "TVector3.h"
35 :
36 12 : ClassImp(AliGenParamPionsKaons)
37 :
38 : //====================================================================================================================================================
39 :
40 : AliGenParamPionsKaons::AliGenParamPionsKaons():
41 0 : AliGenerator(),
42 0 : fGeneratePion(kTRUE),
43 0 : fGenerateKaon(kTRUE),
44 0 : fPtVsRapidityPrimaryPosPions(0x0),
45 0 : fPtVsRapidityPrimaryNegPions(0x0),
46 0 : fPtVsRapidityPrimaryPosKaons(0x0),
47 0 : fPtVsRapidityPrimaryNegKaons(0x0),
48 0 : fHistPdgCode(0x0) {
49 :
50 : // Default constructor
51 :
52 0 : }
53 :
54 : //====================================================================================================================================================
55 :
56 : AliGenParamPionsKaons::AliGenParamPionsKaons(Int_t nPart, Char_t *inputFile):
57 0 : AliGenerator(nPart),
58 0 : fGeneratePion(kTRUE),
59 0 : fGenerateKaon(kTRUE),
60 0 : fPtVsRapidityPrimaryPosPions(0x0),
61 0 : fPtVsRapidityPrimaryNegPions(0x0),
62 0 : fPtVsRapidityPrimaryPosKaons(0x0),
63 0 : fPtVsRapidityPrimaryNegKaons(0x0),
64 0 : fHistPdgCode(0x0) {
65 :
66 : // Standard constructor
67 :
68 0 : fName = "ParamPionsKaons";
69 0 : fTitle = "Parametric pions and kaons generator";
70 :
71 0 : LoadInputHistos(inputFile);
72 :
73 0 : }
74 :
75 : //====================================================================================================================================================
76 :
77 : void AliGenParamPionsKaons::Generate() {
78 :
79 : // Generate one trigger
80 :
81 : Double_t polar[3]= {0,0,0};
82 :
83 0 : Double_t origin[3];
84 0 : Double_t p[3];
85 :
86 0 : Double_t mass=0., pt=0., rap=0., mom=0., energy=0, mt=0., phi=0., time=0.;
87 0 : Int_t nt;
88 : Int_t pdgCode;
89 : Double_t theta = 0.;
90 :
91 0 : for (Int_t j=0; j<3; j++) origin[j] = fOrigin[j];
92 0 : time = fTimeOrigin;
93 0 : if (fVertexSmear==kPerEvent) {
94 0 : Vertex();
95 0 : for (Int_t j=0; j<3; j++) origin[j] = fVertex[j];
96 0 : time = fTime;
97 0 : }
98 :
99 : Int_t nPartGenerated = 0;
100 :
101 0 : while (nPartGenerated<fNpart) {
102 :
103 0 : pdgCode = TMath::Nint(fHistPdgCode->GetRandom());
104 0 : if (TMath::Abs(pdgCode)==321 && !fGenerateKaon) continue;
105 0 : if (TMath::Abs(pdgCode)==211 && !fGeneratePion) continue;
106 :
107 0 : switch (pdgCode) {
108 : case 211:
109 0 : fPtVsRapidityPrimaryPosPions->GetRandom2(pt, rap);
110 0 : break;
111 : case -211:
112 0 : fPtVsRapidityPrimaryNegPions->GetRandom2(pt, rap);
113 0 : break;
114 : case 321:
115 0 : fPtVsRapidityPrimaryPosKaons->GetRandom2(pt, rap);
116 0 : break;
117 : case -321:
118 0 : fPtVsRapidityPrimaryNegKaons->GetRandom2(pt, rap);
119 0 : break;
120 : }
121 :
122 0 : mass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
123 :
124 0 : mt = TMath::Sqrt(mass*mass + pt*pt);
125 0 : energy = mt * TMath::CosH(rap);
126 0 : mom = TMath::Sqrt(energy*energy - mass*mass);
127 :
128 0 : if (TestBit(kYRange)) if (rap<fYMin || rap>fYMax) continue;
129 0 : if (TestBit(kMomentumRange)) if (mom<fPMin || rap>fPMax) continue;
130 0 : if (TestBit(kPtRange)) if (pt<fPtMin || pt>fPtMax) continue;
131 :
132 0 : phi = fPhiMin + gRandom->Rndm()*(fPhiMax-fPhiMin);
133 0 : p[0] = pt*TMath::Cos(phi);
134 0 : p[1] = pt*TMath::Sin(phi);
135 0 : p[2] = mt*TMath::SinH(rap);
136 :
137 0 : TVector3 pv = TVector3(p);
138 0 : theta = pv.Theta();
139 :
140 0 : if (TestBit(kThetaRange)) if (theta<fThetaMin || theta>fThetaMax) continue;
141 :
142 0 : PushTrack(1, -1, Int_t(pdgCode),
143 0 : p[0],p[1],p[2],energy,
144 0 : origin[0],origin[1],origin[2],Double_t(time),
145 : polar[0],polar[1],polar[2],
146 : kPPrimary, nt, 1., 1);
147 :
148 0 : nPartGenerated++;
149 :
150 0 : }
151 :
152 0 : AliGenEventHeader* header = new AliGenEventHeader("ParamPionsKaons");
153 0 : header->SetPrimaryVertex(fVertex);
154 0 : header->SetNProduced(fNpart);
155 0 : header->SetInteractionTime(fTime);
156 :
157 : // Passes header either to the container or to gAlice
158 0 : if (fContainer) {
159 0 : fContainer->AddHeader(header);
160 0 : }
161 : else {
162 0 : gAlice->SetGenEventHeader(header);
163 : }
164 :
165 0 : }
166 :
167 : //====================================================================================================================================================
168 :
169 : void AliGenParamPionsKaons::Init() {
170 :
171 : // Initialisation, check consistency of selected ranges
172 0 : if (TestBit(kPtRange) && TestBit(kMomentumRange))
173 0 : Fatal("Init","You should not set the momentum range and the pt range at the same time!\n");
174 0 : if ((!TestBit(kPtRange)) && (!TestBit(kMomentumRange)))
175 0 : Fatal("Init","You should set either the momentum or the pt range!\n");
176 0 : if ((TestBit(kYRange) && TestBit(kThetaRange)) || (TestBit(kYRange) && TestBit(kEtaRange)) || (TestBit(kEtaRange) && TestBit(kThetaRange)))
177 0 : Fatal("Init","You should only set the range of one of these variables: y, eta or theta\n");
178 0 : if ((!TestBit(kYRange)) && (!TestBit(kEtaRange)) && (!TestBit(kThetaRange)))
179 0 : Fatal("Init","You should set the range of one of these variables: y, eta or theta\n");
180 :
181 0 : AliPDG::AddParticlesToPdgDataBase();
182 :
183 0 : }
184 :
185 : //====================================================================================================================================================
186 :
187 : void AliGenParamPionsKaons::LoadInputHistos(Char_t *inputFile) {
188 :
189 0 : TFile *fileIn = new TFile(inputFile);
190 :
191 0 : TH2D *myPtVsRapidityPrimaryPosPions = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryPosPions");
192 0 : TH2D *myPtVsRapidityPrimaryNegPions = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryNegPions");
193 0 : TH2D *myPtVsRapidityPrimaryPosKaons = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryPosKaons");
194 0 : TH2D *myPtVsRapidityPrimaryNegKaons = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryNegKaons");
195 0 : TH1D *myHistPdgCode = (TH1D*) fileIn->Get("fHistPdgCode");
196 :
197 0 : myPtVsRapidityPrimaryPosPions -> SetName("myPtVsRapidityPrimaryPosPions");
198 0 : myPtVsRapidityPrimaryNegPions -> SetName("myPtVsRapidityPrimaryNegPions");
199 0 : myPtVsRapidityPrimaryPosKaons -> SetName("myPtVsRapidityPrimaryPosKaons");
200 0 : myPtVsRapidityPrimaryNegKaons -> SetName("myPtVsRapidityPrimaryNegKaons");
201 0 : myHistPdgCode -> SetName("myHistPdgCode");
202 :
203 0 : fPtVsRapidityPrimaryPosPions = new TH2D("fPtVsRapidityPrimaryPosPions","",
204 0 : myPtVsRapidityPrimaryPosPions->GetXaxis()->GetNbins(),
205 0 : myPtVsRapidityPrimaryPosPions->GetXaxis()->GetXmin(),
206 0 : myPtVsRapidityPrimaryPosPions->GetXaxis()->GetXmax(),
207 0 : myPtVsRapidityPrimaryPosPions->GetYaxis()->GetNbins(),
208 0 : myPtVsRapidityPrimaryPosPions->GetYaxis()->GetXmin(),
209 0 : myPtVsRapidityPrimaryPosPions->GetYaxis()->GetXmax());
210 0 : for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryPosPions->GetXaxis()->GetNbins(); iBinX++) {
211 0 : for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryPosPions->GetYaxis()->GetNbins(); iBinY++) {
212 0 : fPtVsRapidityPrimaryPosPions->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryPosPions->GetBinContent(iBinX+1,iBinY+1));
213 : }
214 : }
215 :
216 0 : fPtVsRapidityPrimaryNegPions = new TH2D("fPtVsRapidityPrimaryNegPions","",
217 0 : myPtVsRapidityPrimaryNegPions->GetXaxis()->GetNbins(),
218 0 : myPtVsRapidityPrimaryNegPions->GetXaxis()->GetXmin(),
219 0 : myPtVsRapidityPrimaryNegPions->GetXaxis()->GetXmax(),
220 0 : myPtVsRapidityPrimaryNegPions->GetYaxis()->GetNbins(),
221 0 : myPtVsRapidityPrimaryNegPions->GetYaxis()->GetXmin(),
222 0 : myPtVsRapidityPrimaryNegPions->GetYaxis()->GetXmax());
223 0 : for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryNegPions->GetXaxis()->GetNbins(); iBinX++) {
224 0 : for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryNegPions->GetYaxis()->GetNbins(); iBinY++) {
225 0 : fPtVsRapidityPrimaryNegPions->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryNegPions->GetBinContent(iBinX+1,iBinY+1));
226 : }
227 : }
228 :
229 0 : fPtVsRapidityPrimaryPosKaons = new TH2D("fPtVsRapidityPrimaryPosKaons","",
230 0 : myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetNbins(),
231 0 : myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetXmin(),
232 0 : myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetXmax(),
233 0 : myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetNbins(),
234 0 : myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetXmin(),
235 0 : myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetXmax());
236 0 : for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetNbins(); iBinX++) {
237 0 : for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetNbins(); iBinY++) {
238 0 : fPtVsRapidityPrimaryPosKaons->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryPosKaons->GetBinContent(iBinX+1,iBinY+1));
239 : }
240 : }
241 :
242 0 : fPtVsRapidityPrimaryNegKaons = new TH2D("fPtVsRapidityPrimaryNegKaons","",
243 0 : myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetNbins(),
244 0 : myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetXmin(),
245 0 : myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetXmax(),
246 0 : myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetNbins(),
247 0 : myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetXmin(),
248 0 : myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetXmax());
249 0 : for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetNbins(); iBinX++) {
250 0 : for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetNbins(); iBinY++) {
251 0 : fPtVsRapidityPrimaryNegKaons->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryNegKaons->GetBinContent(iBinX+1,iBinY+1));
252 : }
253 : }
254 :
255 0 : fHistPdgCode = new TH1D("fHistPdgCode","",
256 0 : myHistPdgCode->GetXaxis()->GetNbins(),
257 0 : myHistPdgCode->GetXaxis()->GetXmin(),
258 0 : myHistPdgCode->GetXaxis()->GetXmax());
259 0 : for (Int_t iBinX=0; iBinX<myHistPdgCode->GetXaxis()->GetNbins(); iBinX++) {
260 0 : fHistPdgCode->SetBinContent(iBinX+1, myHistPdgCode->GetBinContent(iBinX+1));
261 : }
262 :
263 : // fHistPdgCode = new TH1D("fHistPdgCode", "fHistPdgCode", 10, 0., 10);
264 : // fHistPdgCode -> Fill(3.);
265 :
266 0 : }
267 :
268 : //====================================================================================================================================================
|
