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