Bug in storage manager making possible to delete permanent event fixed
[u/mrichter/AliRoot.git] / MFT / AliGenParamPionsKaons.cxx
CommitLineData
8e94b347 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
35ClassImp(AliGenParamPionsKaons)
36
37//====================================================================================================================================================
38
39AliGenParamPionsKaons::AliGenParamPionsKaons():
40 AliGenerator(),
41 fGeneratePion(kTRUE),
42 fGenerateKaon(kTRUE),
43 fPtVsRapidityPrimaryPosPions(0x0),
44 fPtVsRapidityPrimaryNegPions(0x0),
45 fPtVsRapidityPrimaryPosKaons(0x0),
46 fPtVsRapidityPrimaryNegKaons(0x0),
47 fHistPdgCode(0x0) {
48
49 // Default constructor
50
51}
52
53//====================================================================================================================================================
54
55AliGenParamPionsKaons::AliGenParamPionsKaons(Int_t nPart, Char_t *inputFile):
56 AliGenerator(nPart),
57 fGeneratePion(kTRUE),
58 fGenerateKaon(kTRUE),
59 fPtVsRapidityPrimaryPosPions(0x0),
60 fPtVsRapidityPrimaryNegPions(0x0),
61 fPtVsRapidityPrimaryPosKaons(0x0),
62 fPtVsRapidityPrimaryNegKaons(0x0),
63 fHistPdgCode(0x0) {
64
65 // Standard constructor
66
67 fName = "ParamPionsKaons";
68 fTitle = "Parametric pions and kaons generator";
69
70 LoadInputHistos(inputFile);
71
72}
73
74//====================================================================================================================================================
75
76void AliGenParamPionsKaons::Generate() {
77
78 // Generate one trigger
79
80 Double_t polar[3]= {0,0,0};
81
82 Double_t origin[3];
83 Double_t p[3];
84
85 Double_t mass=0., pt=0., rap=0., mom=0., energy=0, mt=0., phi=0., time=0.;
86 Int_t nt;
87 Int_t pdgCode;
88
89 for (Int_t j=0; j<3; j++) origin[j] = fOrigin[j];
90 time = fTimeOrigin;
91 if (fVertexSmear==kPerEvent) {
92 Vertex();
93 for (Int_t j=0; j<3; j++) origin[j] = fVertex[j];
94 time = fTime;
95 }
96
97 Int_t nPartGenerated = 0;
98
99 while (nPartGenerated<fNpart) {
100
101 pdgCode = TMath::Nint(fHistPdgCode->GetRandom());
102 if (TMath::Abs(pdgCode)==321 && !fGenerateKaon) continue;
103 if (TMath::Abs(pdgCode)==211 && !fGeneratePion) continue;
104
105 mass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
106
107 switch (pdgCode) {
108 case 211:
109 fPtVsRapidityPrimaryPosPions->GetRandom2(pt, rap);
110 break;
111 case -211:
112 fPtVsRapidityPrimaryNegPions->GetRandom2(pt, rap);
113 break;
114 case 321:
115 fPtVsRapidityPrimaryPosKaons->GetRandom2(pt, rap);
116 break;
117 case -321:
118 fPtVsRapidityPrimaryNegKaons->GetRandom2(pt, rap);
119 break;
120 }
121
122 mt = TMath::Sqrt(mass*mass + pt*pt);
123 energy = mt * TMath::CosH(rap);
124 mom = TMath::Sqrt(energy*energy - mass*mass);
125
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;
129
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);
134
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);
140
141 nPartGenerated++;
142
143 }
144
145 AliGenEventHeader* header = new AliGenEventHeader("ParamPionsKaons");
146 header->SetPrimaryVertex(fVertex);
147 header->SetNProduced(fNpart);
148 header->SetInteractionTime(fTime);
149
150 // Passes header either to the container or to gAlice
151 if (fContainer) {
152 fContainer->AddHeader(header);
153 }
154 else {
155 gAlice->SetGenEventHeader(header);
156 }
157
158}
159
160//====================================================================================================================================================
161
162void AliGenParamPionsKaons::Init() {
163
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");
173
174 AliPDG::AddParticlesToPdgDataBase();
175
176}
177
178//====================================================================================================================================================
179
180void AliGenParamPionsKaons::LoadInputHistos(Char_t *inputFile) {
181
182 TFile *fileIn = new TFile(inputFile);
183
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");
189
190 myPtVsRapidityPrimaryPosPions -> SetName("myPtVsRapidityPrimaryPosPions");
191 myPtVsRapidityPrimaryNegPions -> SetName("myPtVsRapidityPrimaryNegPions");
192 myPtVsRapidityPrimaryPosKaons -> SetName("myPtVsRapidityPrimaryPosKaons");
193 myPtVsRapidityPrimaryNegKaons -> SetName("myPtVsRapidityPrimaryNegKaons");
194 myHistPdgCode -> SetName("myHistPdgCode");
195
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));
206 }
207 }
208
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));
219 }
220 }
221
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));
232 }
233 }
234
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));
245 }
246 }
247
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));
254 }
255
256// fHistPdgCode = new TH1D("fHistPdgCode", "fHistPdgCode", 10, 0., 10);
257// fHistPdgCode -> Fill(3.);
258
259}
260
261//====================================================================================================================================================