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