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