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