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