Magnetic field units were changed to Tesla when Reve was moved into ROOT.
[u/mrichter/AliRoot.git] / TDPMjet / AliGenDPMjet.cxx
CommitLineData
1a52e0ed 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// Generator using DPMJET as an external generator
18// The main DPMJET options are accessable for the user through this interface.
19// Uses the TDPMjet implementation of TGenerator.
20
21#include <TDPMjet.h>
22#include <TRandom.h>
23#include <TArrayI.h>
24#include <TParticle.h>
25#include <TGraph.h>
26#include <TDatabasePDG.h>
27#include <TParticlePDG.h>
28#include <TParticleClassPDG.h>
29#include <TPDGCode.h>
30#include <TLorentzVector.h>
9dcc0beb 31#include "AliRunLoader.h"
1a52e0ed 32#include "AliGenDPMjet.h"
33#include "AliGenDPMjetEventHeader.h"
1a52e0ed 34#include "AliRun.h"
f2ba02a4 35#include "AliDpmJetRndm.h"
9dcc0beb 36#include "AliHeader.h"
37#include "AliStack.h"
38#include "AliMC.h"
1a52e0ed 39
f97d4a8e 40ClassImp(AliGenDPMjet)
1a52e0ed 41
42//______________________________________________________________________________
43AliGenDPMjet::AliGenDPMjet()
f97d4a8e 44 :AliGenMC(),
45 fBeamEn(2750.),
46 fEnergyCMS(5500.),
47 fMinImpactParam(0.),
48 fMaxImpactParam(5.),
49 fICentr(0),
50 fSelectAll(0),
51 fFlavor(0),
52 fTrials(0),
53 fSpectators(1),
54 fSpecn(0),
55 fSpecp(0),
56 fDPMjet(0),
57 fParticles(0),
58 fNoGammas(0),
59 fLHC(0),
60 fPi0Decay(0),
61 fGenImpPar(0.),
62 fProcess(kDpmMb)
1a52e0ed 63{
64// Constructor
f2ba02a4 65 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed 66}
67
68
69//______________________________________________________________________________
70AliGenDPMjet::AliGenDPMjet(Int_t npart)
f97d4a8e 71 :AliGenMC(npart),
72 fBeamEn(2750.),
73 fEnergyCMS(5500.),
74 fMinImpactParam(0.),
75 fMaxImpactParam(5.),
76 fICentr(0),
77 fSelectAll(0),
78 fFlavor(0),
79 fTrials(0),
80 fSpectators(1),
81 fSpecn(0),
82 fSpecp(0),
83 fDPMjet(0),
84 fParticles(new TClonesArray("TParticle",10000)),
85 fNoGammas(0),
86 fLHC(0),
87 fPi0Decay(0),
88 fGenImpPar(0.),
89 fProcess(kDpmMb)
1a52e0ed 90{
91// Default PbPb collisions at 5. 5 TeV
92//
93 fName = "DPMJET";
94 fTitle= "Particle Generator using DPMJET";
1a52e0ed 95 SetTarget();
96 SetProjectile();
f2ba02a4 97 fVertex.Set(3);
f2ba02a4 98 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed 99}
100
ba758f5a 101AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
f97d4a8e 102 :AliGenMC(),
103 fBeamEn(2750.),
104 fEnergyCMS(5500.),
105 fMinImpactParam(0.),
106 fMaxImpactParam(5.),
107 fICentr(0),
108 fSelectAll(0),
109 fFlavor(0),
110 fTrials(0),
111 fSpectators(1),
112 fSpecn(0),
113 fSpecp(0),
114 fDPMjet(0),
115 fParticles(0),
116 fNoGammas(0),
117 fLHC(0),
118 fPi0Decay(0),
119 fGenImpPar(0.),
120 fProcess(kDpmMb)
ba758f5a 121{
f97d4a8e 122 // Dummy copy constructor
ba758f5a 123}
1a52e0ed 124
125//______________________________________________________________________________
126AliGenDPMjet::~AliGenDPMjet()
127{
128// Destructor
129 delete fParticles;
130}
1a52e0ed 131//______________________________________________________________________________
132void AliGenDPMjet::Init()
133{
134// Initialization
135
e67e0fb7 136 SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
1a52e0ed 137 fBeamEn,fEnergyCMS));
138
35c87559 139 fDPMjet=(TDPMjet*) fMCEvGen;
1a52e0ed 140 //
141 // **** Flag to force central production
142 // fICentr=1. central production forced
143 // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam
144 // fICentr<-99 -> fraction of x-sec. = XSFRAC
145 // fICentr=-1. -> evaporation/fzc suppressed
146 // fICentr<-1. -> evaporation/fzc suppressed
ba758f5a 147 if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1);
148
1a52e0ed 149 fDPMjet->SetfFCentr(fICentr);
e67e0fb7 150 fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
f97d4a8e 151 fDPMjet->SetPi0Decay(fPi0Decay);
1a52e0ed 152//
153// Initialize DPMjet
154//
155 fDPMjet->Initialize();
1a52e0ed 156}
157
158
159//______________________________________________________________________________
160void AliGenDPMjet::Generate()
161{
162// Generate one event
163
164 Float_t polar[3] = {0,0,0};
165 Float_t origin[3] = {0,0,0};
166 Float_t origin0[3] = {0,0,0};
e67e0fb7 167 Float_t p[3];
1a52e0ed 168 Float_t tof;
169
170// converts from mm/c to s
171 const Float_t kconv = 0.001/2.999792458e8;
172 Int_t nt = 0;
173 Int_t jev = 0;
e67e0fb7 174 Int_t kf, ks, imo;
1a52e0ed 175 kf = 0;
1a52e0ed 176 fTrials = 0;
e67e0fb7 177 // Set collision vertex position
178 if (fVertexSmear == kPerEvent) Vertex();
1a52e0ed 179
180 while(1)
181 {
182// Generate one event
183// --------------------------------------------------------------------------
184 fSpecn = 0;
185 fSpecp = 0;
186// --------------------------------------------------------------------------
187 fDPMjet->GenerateEvent();
188 fTrials++;
189
190 fDPMjet->ImportParticles(fParticles,"All");
191 if (fLHC) Boost();
192
193 // Temporaneo
194 fGenImpPar = fDPMjet->GetBImpac();
195
196 Int_t np = fParticles->GetEntriesFast();
197 printf("\n **************************************************%d\n",np);
198 Int_t nc=0;
199 if (np==0) continue;
200 Int_t i;
201 Int_t* newPos = new Int_t[np];
202 Int_t* pSelected = new Int_t[np];
203
204 for (i = 0; i<np; i++) {
205 newPos[i] = i;
206 pSelected[i] = 0;
207 }
208
209// Get event vertex
210
f2ba02a4 211 fVertex[0] = origin0[0];
212 fVertex[1] = origin0[1];
213 fVertex[2] = origin0[2];
1a52e0ed 214
215// First select parent particles
216
217 for (i = 0; i<np; i++) {
218 TParticle *iparticle = (TParticle *) fParticles->At(i);
1a52e0ed 219
220// Is this a parent particle ?
221
222 if (Stable(iparticle)) continue;
223
224 Bool_t selected = kTRUE;
225 Bool_t hasSelectedDaughters = kFALSE;
226
227 kf = iparticle->GetPdgCode();
228 if (kf == 92) continue;
229 ks = iparticle->GetStatusCode();
230// No initial state partons
231 if (ks==21) continue;
232
233 if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
234 SelectFlavor(kf);
235 hasSelectedDaughters = DaughtersSelection(iparticle);
236
237// Put particle on the stack if it is either selected or
238// it is the mother of at least one seleted particle
239
240 if (selected || hasSelectedDaughters) {
241 nc++;
242 pSelected[i] = 1;
243 } // selected
244 } // particle loop parents
245
246// Now select the final state particles
247
248
249 for (i=0; i<np; i++) {
250 TParticle *iparticle = (TParticle *) fParticles->At(i);
1a52e0ed 251
252// Is this a final state particle ?
253
254 if (!Stable(iparticle)) continue;
255
256 Bool_t selected = kTRUE;
257 kf = iparticle->GetPdgCode();
258 ks = iparticle->GetStatusCode();
259
260// --------------------------------------------------------------------------
261// Count spectator neutrons and protons (ks == 13, 14)
262 if(ks == 13 || ks == 14){
263 if(kf == kNeutron) fSpecn += 1;
264 if(kf == kProton) fSpecp += 1;
265 }
266// --------------------------------------------------------------------------
267
268 if (!fSelectAll) {
269 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
270 if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
271 }
272
273// Put particle on the stack if selected
274
275 if (selected) {
276 nc++;
277 pSelected[i] = 1;
278 } // selected
279 } // particle loop final state
280
281// Write particles to stack
282
283 for (i = 0; i<np; i++) {
284 TParticle * iparticle = (TParticle *) fParticles->At(i);
285 Bool_t hasMother = (iparticle->GetFirstMother()>=0);
1a52e0ed 286 if (pSelected[i]) {
287
288 kf = iparticle->GetPdgCode();
1a52e0ed 289 ks = iparticle->GetStatusCode();
290
291 p[0] = iparticle->Px();
292 p[1] = iparticle->Py();
293 p[2] = iparticle->Pz();
e67e0fb7 294 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
295 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
296 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
297
1a52e0ed 298 tof = kconv*iparticle->T();
299
300 imo = -1;
301 TParticle* mother = 0;
302 if (hasMother) {
303 imo = iparticle->GetFirstMother();
304 mother = (TParticle *) fParticles->At(imo);
305 imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
306 } // if has mother
307
1a52e0ed 308 Bool_t tFlag = (fTrackIt && (ks == 1));
309
642f15cf 310 PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
1a52e0ed 311 KeepTrack(nt);
312 newPos[i] = nt;
313 } // if selected
314 } // particle loop
315 delete[] newPos;
316 delete[] pSelected;
1a52e0ed 317 if (nc>0) {
318 jev += nc;
319 if (jev >= fNpart || fNpart == -1) {
1a52e0ed 320 break;
321 }
1a52e0ed 322 }
323 } // event loop
324 MakeHeader();
325 SetHighWaterMark(nt);
326}
327
1a52e0ed 328//______________________________________________________________________________
329Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
330{
331//
332// Looks recursively if one of the daughters has been selected
333//
334// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
335 Int_t imin = -1;
336 Int_t imax = -1;
337 Int_t i;
338 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
339 Bool_t selected = kFALSE;
340 if (hasDaughters) {
341 imin = iparticle->GetFirstDaughter();
342 imax = iparticle->GetLastDaughter();
343 for (i = imin; i <= imax; i++){
344 TParticle * jparticle = (TParticle *) fParticles->At(i);
345 Int_t ip = jparticle->GetPdgCode();
346 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
347 selected=kTRUE; break;
348 }
349 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
350 }
351 } else {
352 return kFALSE;
353 }
354 return selected;
355}
356
357
358
359//______________________________________________________________________________
360Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
361{
362// Select flavor of particle
363// 0: all
364// 4: charm and beauty
365// 5: beauty
366 Bool_t res = 0;
367
368 if (fFlavor == 0) {
369 res = kTRUE;
370 } else {
371 Int_t ifl = TMath::Abs(pid/100);
372 if (ifl > 10) ifl/=10;
373 res = (fFlavor == ifl);
374 }
375//
376// This part if gamma writing is inhibited
377 if (fNoGammas)
378 res = res && (pid != kGamma && pid != kPi0);
379//
380 return res;
381}
382
383//______________________________________________________________________________
384Bool_t AliGenDPMjet::Stable(TParticle* particle)
385{
386// Return true for a stable particle
387//
388
389// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
390 if (particle->GetStatusCode() == 1) return kTRUE;
391 else return kFALSE;
392
393}
394
1a52e0ed 395//______________________________________________________________________________
396void AliGenDPMjet::MakeHeader()
397{
398// Builds the event header, to be called after each event
399 AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
400 ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
401 ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
402 ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
1a52e0ed 403 ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
404 fDPMjet->GetfIt());
9dcc0beb 405 ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
1a52e0ed 406// Bookkeeping for kinematic bias
407 ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
408// Event Vertex
f2ba02a4 409 header->SetPrimaryVertex(fVertex);
1a52e0ed 410 gAlice->SetGenEventHeader(header);
9dcc0beb 411 AddHeader(header);
1a52e0ed 412}
413
9dcc0beb 414void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
415{
416 // Add header to container or runloader
417 if (fContainer) {
418 fContainer->AddHeader(header);
419 } else {
420 AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
421 }
422}
1a52e0ed 423
424
425//______________________________________________________________________________
ba758f5a 426AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /*rhs*/)
1a52e0ed 427{
428// Assignment operator
429 return *this;
430}
431
432
433
434//______________________________________________________________________________