Adding variable comments
[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.),
f97d4a8e 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),
f97d4a8e 57 fNoGammas(0),
58 fLHC(0),
9c97cedd 59 fPi0Decay(1),
717765ce 60 fDecayAll(0),
f97d4a8e 61 fGenImpPar(0.),
62 fProcess(kDpmMb)
1a52e0ed 63{
64// Constructor
e7c989e4 65 fEnergyCMS = 5500.;
f2ba02a4 66 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed 67}
68
69
70//______________________________________________________________________________
71AliGenDPMjet::AliGenDPMjet(Int_t npart)
f97d4a8e 72 :AliGenMC(npart),
73 fBeamEn(2750.),
f97d4a8e 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),
f97d4a8e 84 fNoGammas(0),
85 fLHC(0),
9c97cedd 86 fPi0Decay(1),
717765ce 87 fDecayAll(0),
f97d4a8e 88 fGenImpPar(0.),
89 fProcess(kDpmMb)
1a52e0ed 90{
91// Default PbPb collisions at 5. 5 TeV
92//
e7c989e4 93 fEnergyCMS = 5500.;
1a52e0ed 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.),
f97d4a8e 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),
f97d4a8e 115 fNoGammas(0),
116 fLHC(0),
9c97cedd 117 fPi0Decay(1),
717765ce 118 fDecayAll(0),
f97d4a8e 119 fGenImpPar(0.),
120 fProcess(kDpmMb)
ba758f5a 121{
f97d4a8e 122 // Dummy copy constructor
e7c989e4 123 fEnergyCMS = 5500.;
ba758f5a 124}
1a52e0ed 125
126//______________________________________________________________________________
127AliGenDPMjet::~AliGenDPMjet()
128{
129// Destructor
1a52e0ed 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
cd19c617 147 if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
ba758f5a 148
1a52e0ed 149 fDPMjet->SetfFCentr(fICentr);
e67e0fb7 150 fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
f97d4a8e 151 fDPMjet->SetPi0Decay(fPi0Decay);
717765ce 152 fDPMjet->SetDecayAll(fDecayAll);
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};
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();
9c97cedd 188
1a52e0ed 189 fTrials++;
190
8507138f 191 fDPMjet->ImportParticles(&fParticles,"All");
1a52e0ed 192 if (fLHC) Boost();
193
194 // Temporaneo
195 fGenImpPar = fDPMjet->GetBImpac();
196
8507138f 197 Int_t np = fParticles.GetEntriesFast();
1a52e0ed 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
1a52e0ed 210// First select parent particles
211
212 for (i = 0; i<np; i++) {
8507138f 213 TParticle *iparticle = (TParticle *) fParticles.At(i);
1a52e0ed 214
215// Is this a parent particle ?
216
217 if (Stable(iparticle)) continue;
218
219 Bool_t selected = kTRUE;
220 Bool_t hasSelectedDaughters = kFALSE;
221
222 kf = iparticle->GetPdgCode();
cd19c617 223 if (kf == 92 || kf == 99999) continue;
1a52e0ed 224 ks = iparticle->GetStatusCode();
225// No initial state partons
226 if (ks==21) continue;
227
228 if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
229 SelectFlavor(kf);
9c97cedd 230
231
1a52e0ed 232 hasSelectedDaughters = DaughtersSelection(iparticle);
233
9c97cedd 234
1a52e0ed 235// Put particle on the stack if it is either selected or
236// it is the mother of at least one seleted particle
237
238 if (selected || hasSelectedDaughters) {
239 nc++;
240 pSelected[i] = 1;
241 } // selected
242 } // particle loop parents
243
244// Now select the final state particles
245
246
247 for (i=0; i<np; i++) {
8507138f 248 TParticle *iparticle = (TParticle *) fParticles.At(i);
1a52e0ed 249
250// Is this a final state particle ?
251
252 if (!Stable(iparticle)) continue;
253
254 Bool_t selected = kTRUE;
255 kf = iparticle->GetPdgCode();
256 ks = iparticle->GetStatusCode();
257
258// --------------------------------------------------------------------------
259// Count spectator neutrons and protons (ks == 13, 14)
260 if(ks == 13 || ks == 14){
261 if(kf == kNeutron) fSpecn += 1;
262 if(kf == kProton) fSpecp += 1;
263 }
264// --------------------------------------------------------------------------
265
266 if (!fSelectAll) {
267 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
268 if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
269 }
270
271// Put particle on the stack if selected
272
273 if (selected) {
274 nc++;
275 pSelected[i] = 1;
276 } // selected
277 } // particle loop final state
278
279// Write particles to stack
280
281 for (i = 0; i<np; i++) {
8507138f 282 TParticle * iparticle = (TParticle *) fParticles.At(i);
1a52e0ed 283 Bool_t hasMother = (iparticle->GetFirstMother()>=0);
1a52e0ed 284 if (pSelected[i]) {
285
286 kf = iparticle->GetPdgCode();
1a52e0ed 287 ks = iparticle->GetStatusCode();
288
289 p[0] = iparticle->Px();
290 p[1] = iparticle->Py();
291 p[2] = iparticle->Pz();
e67e0fb7 292 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
293 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
294 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
295
1a52e0ed 296 tof = kconv*iparticle->T();
297
298 imo = -1;
299 TParticle* mother = 0;
300 if (hasMother) {
301 imo = iparticle->GetFirstMother();
8507138f 302 mother = (TParticle *) fParticles.At(imo);
cd19c617 303 imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
1a52e0ed 304 } // if has mother
305
9c97cedd 306
307
1a52e0ed 308 Bool_t tFlag = (fTrackIt && (ks == 1));
642f15cf 309 PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
1a52e0ed 310 KeepTrack(nt);
311 newPos[i] = nt;
312 } // if selected
313 } // particle loop
314 delete[] newPos;
315 delete[] pSelected;
1a52e0ed 316 if (nc>0) {
317 jev += nc;
318 if (jev >= fNpart || fNpart == -1) {
1a52e0ed 319 break;
320 }
1a52e0ed 321 }
322 } // event loop
323 MakeHeader();
324 SetHighWaterMark(nt);
325}
326
1a52e0ed 327//______________________________________________________________________________
328Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
329{
330//
331// Looks recursively if one of the daughters has been selected
332//
333// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
334 Int_t imin = -1;
335 Int_t imax = -1;
336 Int_t i;
337 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
338 Bool_t selected = kFALSE;
339 if (hasDaughters) {
340 imin = iparticle->GetFirstDaughter();
341 imax = iparticle->GetLastDaughter();
342 for (i = imin; i <= imax; i++){
8507138f 343 TParticle * jparticle = (TParticle *) fParticles.At(i);
1a52e0ed 344 Int_t ip = jparticle->GetPdgCode();
345 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
346 selected=kTRUE; break;
347 }
348 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
349 }
350 } else {
351 return kFALSE;
352 }
353 return selected;
354}
355
356
357
358//______________________________________________________________________________
359Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
360{
361// Select flavor of particle
362// 0: all
363// 4: charm and beauty
364// 5: beauty
365 Bool_t res = 0;
366
367 if (fFlavor == 0) {
368 res = kTRUE;
369 } else {
370 Int_t ifl = TMath::Abs(pid/100);
371 if (ifl > 10) ifl/=10;
372 res = (fFlavor == ifl);
373 }
374//
375// This part if gamma writing is inhibited
376 if (fNoGammas)
377 res = res && (pid != kGamma && pid != kPi0);
378//
379 return res;
380}
381
382//______________________________________________________________________________
383Bool_t AliGenDPMjet::Stable(TParticle* particle)
384{
385// Return true for a stable particle
386//
387
388// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
389 if (particle->GetStatusCode() == 1) return kTRUE;
390 else return kFALSE;
391
392}
393
1a52e0ed 394//______________________________________________________________________________
395void AliGenDPMjet::MakeHeader()
396{
397// Builds the event header, to be called after each event
398 AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
399 ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
400 ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
401 ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
1a52e0ed 402 ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
403 fDPMjet->GetfIt());
9dcc0beb 404 ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
1a52e0ed 405// Bookkeeping for kinematic bias
406 ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
407// Event Vertex
f2ba02a4 408 header->SetPrimaryVertex(fVertex);
1a52e0ed 409 gAlice->SetGenEventHeader(header);
d2593e7a 410 AddHeader(header);
1a52e0ed 411}
412
9dcc0beb 413void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
414{
415 // Add header to container or runloader
416 if (fContainer) {
417 fContainer->AddHeader(header);
418 } else {
33c3c91a 419 AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
9dcc0beb 420 }
421}
1a52e0ed 422
423
424//______________________________________________________________________________
ba758f5a 425AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /*rhs*/)
1a52e0ed 426{
427// Assignment operator
428 return *this;
429}
430
431
d2593e7a 432void AliGenDPMjet::FinishRun()
433{
434 // Print run statistics
435 fDPMjet->Dt_Dtuout();
436}
437
438
1a52e0ed 439
440//______________________________________________________________________________