Fix for #84631: Lambda_c->(Lambda(1520)->p,K),pi branching ratio
[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
c0d429ef 174 Double_t polar[3] = {0,0,0};
175 Double_t origin[3] = {0,0,0};
176 Double_t p[4] = {0};
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();
c0d429ef 333 p[3] = iparticle->Energy();
e67e0fb7 334 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
335 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
336 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
337
21391258 338 tof = fTime + kconv*iparticle->T();
1a52e0ed 339
340 imo = -1;
341 TParticle* mother = 0;
342 if (hasMother) {
343 imo = iparticle->GetFirstMother();
8507138f 344 mother = (TParticle *) fParticles.At(imo);
cd19c617 345 imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
1a52e0ed 346 } // if has mother
347
9c97cedd 348
349
1a52e0ed 350 Bool_t tFlag = (fTrackIt && (ks == 1));
c0d429ef 351 PushTrack(tFlag, imo, kf,
352 p[0], p[1], p[2], p[3],
353 origin[0], origin[1], origin[2], tof,
354 polar[0], polar[1], polar[2],
355 kPNoProcess, nt, 1., ks);
1a52e0ed 356 KeepTrack(nt);
357 newPos[i] = nt;
358 } // if selected
359 } // particle loop
360 delete[] newPos;
361 delete[] pSelected;
1a52e0ed 362 if (nc>0) {
363 jev += nc;
364 if (jev >= fNpart || fNpart == -1) {
1a52e0ed 365 break;
366 }
1a52e0ed 367 }
368 } // event loop
369 MakeHeader();
370 SetHighWaterMark(nt);
371}
372
1a52e0ed 373//______________________________________________________________________________
374Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
375{
376//
377// Looks recursively if one of the daughters has been selected
378//
379// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
380 Int_t imin = -1;
381 Int_t imax = -1;
382 Int_t i;
383 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
384 Bool_t selected = kFALSE;
385 if (hasDaughters) {
386 imin = iparticle->GetFirstDaughter();
387 imax = iparticle->GetLastDaughter();
388 for (i = imin; i <= imax; i++){
8507138f 389 TParticle * jparticle = (TParticle *) fParticles.At(i);
1a52e0ed 390 Int_t ip = jparticle->GetPdgCode();
391 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
392 selected=kTRUE; break;
393 }
394 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
395 }
396 } else {
397 return kFALSE;
398 }
399 return selected;
400}
401
402
403
404//______________________________________________________________________________
405Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
406{
407// Select flavor of particle
408// 0: all
409// 4: charm and beauty
410// 5: beauty
411 Bool_t res = 0;
412
413 if (fFlavor == 0) {
414 res = kTRUE;
415 } else {
416 Int_t ifl = TMath::Abs(pid/100);
417 if (ifl > 10) ifl/=10;
418 res = (fFlavor == ifl);
419 }
420//
421// This part if gamma writing is inhibited
422 if (fNoGammas)
423 res = res && (pid != kGamma && pid != kPi0);
424//
425 return res;
426}
427
428//______________________________________________________________________________
429Bool_t AliGenDPMjet::Stable(TParticle* particle)
430{
431// Return true for a stable particle
432//
433
434// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
435 if (particle->GetStatusCode() == 1) return kTRUE;
436 else return kFALSE;
437
438}
439
1a52e0ed 440//______________________________________________________________________________
441void AliGenDPMjet::MakeHeader()
442{
4eb21a7b 443 printf("MakeHeader %13.3f \n", fDPMjet->GetBImpac());
1a52e0ed 444// Builds the event header, to be called after each event
445 AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
446 ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
447 ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
448 ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
3c349853 449 ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetProjParticipants(),
450 fDPMjet->GetTargParticipants());
451 ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
452 // Bookkeeping for kinematic bias
1a52e0ed 453 ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
3c349853 454 // Event Vertex
f2ba02a4 455 header->SetPrimaryVertex(fVertex);
21391258 456 header->SetInteractionTime(fTime);
1a52e0ed 457 gAlice->SetGenEventHeader(header);
d2593e7a 458 AddHeader(header);
1a52e0ed 459}
460
9dcc0beb 461void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
462{
463 // Add header to container or runloader
464 if (fContainer) {
465 fContainer->AddHeader(header);
466 } else {
33c3c91a 467 AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
9dcc0beb 468 }
469}
1a52e0ed 470
471
472//______________________________________________________________________________
ba758f5a 473AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /*rhs*/)
1a52e0ed 474{
475// Assignment operator
476 return *this;
477}
478
479
d2593e7a 480void AliGenDPMjet::FinishRun()
481{
482 // Print run statistics
483 fDPMjet->Dt_Dtuout();
484}
485
486
1a52e0ed 487
488//______________________________________________________________________________