AliRoot Core  a565103 (a565103)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliGenParamPionsKaons.cxx
Go to the documentation of this file.
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 
37 
38 //====================================================================================================================================================
39 
41  AliGenerator(),
42  fGeneratePion(kTRUE),
43  fGenerateKaon(kTRUE),
44  fPtVsRapidityPrimaryPosPions(0x0),
45  fPtVsRapidityPrimaryNegPions(0x0),
46  fPtVsRapidityPrimaryPosKaons(0x0),
47  fPtVsRapidityPrimaryNegKaons(0x0),
48  fHistPdgCode(0x0) {
49 
50  // Default constructor
51 
52 }
53 
54 //====================================================================================================================================================
55 
56 AliGenParamPionsKaons::AliGenParamPionsKaons(Int_t nPart, Char_t *inputFile):
57  AliGenerator(nPart),
58  fGeneratePion(kTRUE),
59  fGenerateKaon(kTRUE),
60  fPtVsRapidityPrimaryPosPions(0x0),
61  fPtVsRapidityPrimaryNegPions(0x0),
62  fPtVsRapidityPrimaryPosKaons(0x0),
63  fPtVsRapidityPrimaryNegKaons(0x0),
64  fHistPdgCode(0x0) {
65 
66  // Standard constructor
67 
68  fName = "ParamPionsKaons";
69  fTitle = "Parametric pions and kaons generator";
70 
71  LoadInputHistos(inputFile);
72 
73 }
74 
75 //====================================================================================================================================================
76 
78 
79  // Generate one trigger
80 
81  Double_t polar[3]= {0,0,0};
82 
83  Double_t origin[3];
84  Double_t p[3];
85 
86  Double_t mass=0., pt=0., rap=0., mom=0., energy=0, mt=0., phi=0., time=0.;
87  Int_t nt;
88  Int_t pdgCode;
89  Double_t theta = 0.;
90 
91  for (Int_t j=0; j<3; j++) origin[j] = fOrigin[j];
92  time = fTimeOrigin;
93  if (fVertexSmear==kPerEvent) {
94  Vertex();
95  for (Int_t j=0; j<3; j++) origin[j] = fVertex[j];
96  time = fTime;
97  }
98 
99  Int_t nPartGenerated = 0;
100 
101  while (nPartGenerated<fNpart) {
102 
103  pdgCode = TMath::Nint(fHistPdgCode->GetRandom());
104  if (TMath::Abs(pdgCode)==321 && !fGenerateKaon) continue;
105  if (TMath::Abs(pdgCode)==211 && !fGeneratePion) continue;
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  mass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
123 
124  mt = TMath::Sqrt(mass*mass + pt*pt);
125  energy = mt * TMath::CosH(rap);
126  mom = TMath::Sqrt(energy*energy - mass*mass);
127 
128  if (TestBit(kYRange)) if (rap<fYMin || rap>fYMax) continue;
129  if (TestBit(kMomentumRange)) if (mom<fPMin || rap>fPMax) continue;
130  if (TestBit(kPtRange)) if (pt<fPtMin || pt>fPtMax) continue;
131 
132  phi = fPhiMin + gRandom->Rndm()*(fPhiMax-fPhiMin);
133  p[0] = pt*TMath::Cos(phi);
134  p[1] = pt*TMath::Sin(phi);
135  p[2] = mt*TMath::SinH(rap);
136 
137  TVector3 pv = TVector3(p);
138  theta = pv.Theta();
139 
140  if (TestBit(kThetaRange)) if (theta<fThetaMin || theta>fThetaMax) continue;
141 
142  PushTrack(1, -1, Int_t(pdgCode),
143  p[0],p[1],p[2],energy,
144  origin[0],origin[1],origin[2],Double_t(time),
145  polar[0],polar[1],polar[2],
146  kPPrimary, nt, 1., 1);
147 
148  nPartGenerated++;
149 
150  }
151 
152  AliGenEventHeader* header = new AliGenEventHeader("ParamPionsKaons");
153  header->SetPrimaryVertex(fVertex);
154  header->SetNProduced(fNpart);
155  header->SetInteractionTime(fTime);
156 
157  // Passes header either to the container or to gAlice
158  if (fContainer) {
159  fContainer->AddHeader(header);
160  }
161  else {
162  gAlice->SetGenEventHeader(header);
163  }
164 
165 }
166 
167 //====================================================================================================================================================
168 
170 
171  // Initialisation, check consistency of selected ranges
172  if (TestBit(kPtRange) && TestBit(kMomentumRange))
173  Fatal("Init","You should not set the momentum range and the pt range at the same time!\n");
174  if ((!TestBit(kPtRange)) && (!TestBit(kMomentumRange)))
175  Fatal("Init","You should set either the momentum or the pt range!\n");
176  if ((TestBit(kYRange) && TestBit(kThetaRange)) || (TestBit(kYRange) && TestBit(kEtaRange)) || (TestBit(kEtaRange) && TestBit(kThetaRange)))
177  Fatal("Init","You should only set the range of one of these variables: y, eta or theta\n");
178  if ((!TestBit(kYRange)) && (!TestBit(kEtaRange)) && (!TestBit(kThetaRange)))
179  Fatal("Init","You should set the range of one of these variables: y, eta or theta\n");
180 
181  AliPDG::AddParticlesToPdgDataBase();
182 
183 }
184 
185 //====================================================================================================================================================
186 
187 void AliGenParamPionsKaons::LoadInputHistos(Char_t *inputFile) {
188 
189  TFile *fileIn = new TFile(inputFile);
190 
191  TH2D *myPtVsRapidityPrimaryPosPions = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryPosPions");
192  TH2D *myPtVsRapidityPrimaryNegPions = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryNegPions");
193  TH2D *myPtVsRapidityPrimaryPosKaons = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryPosKaons");
194  TH2D *myPtVsRapidityPrimaryNegKaons = (TH2D*) fileIn->Get("fPtVsRapidityPrimaryNegKaons");
195  TH1D *myHistPdgCode = (TH1D*) fileIn->Get("fHistPdgCode");
196 
197  myPtVsRapidityPrimaryPosPions -> SetName("myPtVsRapidityPrimaryPosPions");
198  myPtVsRapidityPrimaryNegPions -> SetName("myPtVsRapidityPrimaryNegPions");
199  myPtVsRapidityPrimaryPosKaons -> SetName("myPtVsRapidityPrimaryPosKaons");
200  myPtVsRapidityPrimaryNegKaons -> SetName("myPtVsRapidityPrimaryNegKaons");
201  myHistPdgCode -> SetName("myHistPdgCode");
202 
203  fPtVsRapidityPrimaryPosPions = new TH2D("fPtVsRapidityPrimaryPosPions","",
204  myPtVsRapidityPrimaryPosPions->GetXaxis()->GetNbins(),
205  myPtVsRapidityPrimaryPosPions->GetXaxis()->GetXmin(),
206  myPtVsRapidityPrimaryPosPions->GetXaxis()->GetXmax(),
207  myPtVsRapidityPrimaryPosPions->GetYaxis()->GetNbins(),
208  myPtVsRapidityPrimaryPosPions->GetYaxis()->GetXmin(),
209  myPtVsRapidityPrimaryPosPions->GetYaxis()->GetXmax());
210  for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryPosPions->GetXaxis()->GetNbins(); iBinX++) {
211  for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryPosPions->GetYaxis()->GetNbins(); iBinY++) {
212  fPtVsRapidityPrimaryPosPions->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryPosPions->GetBinContent(iBinX+1,iBinY+1));
213  }
214  }
215 
216  fPtVsRapidityPrimaryNegPions = new TH2D("fPtVsRapidityPrimaryNegPions","",
217  myPtVsRapidityPrimaryNegPions->GetXaxis()->GetNbins(),
218  myPtVsRapidityPrimaryNegPions->GetXaxis()->GetXmin(),
219  myPtVsRapidityPrimaryNegPions->GetXaxis()->GetXmax(),
220  myPtVsRapidityPrimaryNegPions->GetYaxis()->GetNbins(),
221  myPtVsRapidityPrimaryNegPions->GetYaxis()->GetXmin(),
222  myPtVsRapidityPrimaryNegPions->GetYaxis()->GetXmax());
223  for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryNegPions->GetXaxis()->GetNbins(); iBinX++) {
224  for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryNegPions->GetYaxis()->GetNbins(); iBinY++) {
225  fPtVsRapidityPrimaryNegPions->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryNegPions->GetBinContent(iBinX+1,iBinY+1));
226  }
227  }
228 
229  fPtVsRapidityPrimaryPosKaons = new TH2D("fPtVsRapidityPrimaryPosKaons","",
230  myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetNbins(),
231  myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetXmin(),
232  myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetXmax(),
233  myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetNbins(),
234  myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetXmin(),
235  myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetXmax());
236  for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryPosKaons->GetXaxis()->GetNbins(); iBinX++) {
237  for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryPosKaons->GetYaxis()->GetNbins(); iBinY++) {
238  fPtVsRapidityPrimaryPosKaons->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryPosKaons->GetBinContent(iBinX+1,iBinY+1));
239  }
240  }
241 
242  fPtVsRapidityPrimaryNegKaons = new TH2D("fPtVsRapidityPrimaryNegKaons","",
243  myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetNbins(),
244  myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetXmin(),
245  myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetXmax(),
246  myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetNbins(),
247  myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetXmin(),
248  myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetXmax());
249  for (Int_t iBinX=0; iBinX<myPtVsRapidityPrimaryNegKaons->GetXaxis()->GetNbins(); iBinX++) {
250  for (Int_t iBinY=0; iBinY<myPtVsRapidityPrimaryNegKaons->GetYaxis()->GetNbins(); iBinY++) {
251  fPtVsRapidityPrimaryNegKaons->SetBinContent(iBinX+1, iBinY+1, myPtVsRapidityPrimaryNegKaons->GetBinContent(iBinX+1,iBinY+1));
252  }
253  }
254 
255  fHistPdgCode = new TH1D("fHistPdgCode","",
256  myHistPdgCode->GetXaxis()->GetNbins(),
257  myHistPdgCode->GetXaxis()->GetXmin(),
258  myHistPdgCode->GetXaxis()->GetXmax());
259  for (Int_t iBinX=0; iBinX<myHistPdgCode->GetXaxis()->GetNbins(); iBinX++) {
260  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 }
267 
268 //====================================================================================================================================================
Float_t fThetaMax
Definition: fastMUONSim.C:57
ClassImp(AliGenParamPionsKaons) AliGenParamPionsKaons
AliRun * gAlice
virtual void LoadInputHistos(Char_t *inputFile)