dad713bc320af47597952f981ac5bc99ebb0ace6
[u/mrichter/AliRoot.git] / MFT / AliGenParamPionsKaons.cxx
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
35 ClassImp(AliGenParamPionsKaons)
36
37 //====================================================================================================================================================
38
39 AliGenParamPionsKaons::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
55 AliGenParamPionsKaons::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
76 void 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
162 void 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
180 void 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 //====================================================================================================================================================