added a couple of functions
[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>
31
32#include "AliGenDPMjet.h"
33#include "AliGenDPMjetEventHeader.h"
34#include "AliPythia.h"
35#include "AliRun.h"
f2ba02a4 36#include "AliDpmJetRndm.h"
1a52e0ed 37
38 ClassImp(AliGenDPMjet)
39
40
41//______________________________________________________________________________
42AliGenDPMjet::AliGenDPMjet()
43 :AliGenMC()
44{
45// Constructor
46 fParticles = 0;
47 fDPMjet = 0;
f2ba02a4 48 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed 49}
50
51
52//______________________________________________________________________________
53AliGenDPMjet::AliGenDPMjet(Int_t npart)
54 :AliGenMC(npart)
55{
56// Default PbPb collisions at 5. 5 TeV
57//
58 fName = "DPMJET";
59 fTitle= "Particle Generator using DPMJET";
60
1a52e0ed 61 SetEnergyCMS();
62 SetTarget();
63 SetProjectile();
64 SetCentral();
65 SetImpactParameterRange();
66 SetBoostLHC();
67
68 fKeep = 0;
69 fDecaysOff = 1;
70 fEvaluate = 0;
71 fSelectAll = 0;
72 fFlavor = 0;
73 fSpectators = 1;
f2ba02a4 74 fVertex.Set(3);
1a52e0ed 75
76 fParticles = new TClonesArray("TParticle",10000);
77
78 // Set random number generator
1a52e0ed 79 fDPMjet = 0;
80 // Instance AliPythia
81 AliPythia::Instance();
f2ba02a4 82 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed 83}
84
ba758f5a 85AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
86:AliGenMC()
87{
88}
1a52e0ed 89
90//______________________________________________________________________________
91AliGenDPMjet::~AliGenDPMjet()
92{
93// Destructor
94 delete fParticles;
95}
96
97
98//______________________________________________________________________________
99void AliGenDPMjet::Init()
100{
101// Initialization
102
e67e0fb7 103 SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
1a52e0ed 104 fBeamEn,fEnergyCMS));
105
35c87559 106 fDPMjet=(TDPMjet*) fMCEvGen;
1a52e0ed 107 //
108 // **** Flag to force central production
109 // fICentr=1. central production forced
110 // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam
111 // fICentr<-99 -> fraction of x-sec. = XSFRAC
112 // fICentr=-1. -> evaporation/fzc suppressed
113 // fICentr<-1. -> evaporation/fzc suppressed
ba758f5a 114 if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1);
115
1a52e0ed 116 fDPMjet->SetfFCentr(fICentr);
e67e0fb7 117 fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
1a52e0ed 118
119//
120// Initialize DPMjet
121//
122 fDPMjet->Initialize();
1a52e0ed 123}
124
125
126//______________________________________________________________________________
127void AliGenDPMjet::Generate()
128{
129// Generate one event
130
131 Float_t polar[3] = {0,0,0};
132 Float_t origin[3] = {0,0,0};
133 Float_t origin0[3] = {0,0,0};
e67e0fb7 134 Float_t p[3];
1a52e0ed 135 Float_t tof;
136
137// converts from mm/c to s
138 const Float_t kconv = 0.001/2.999792458e8;
139 Int_t nt = 0;
140 Int_t jev = 0;
e67e0fb7 141 Int_t kf, ks, imo;
1a52e0ed 142 kf = 0;
1a52e0ed 143 fTrials = 0;
e67e0fb7 144 // Set collision vertex position
145 if (fVertexSmear == kPerEvent) Vertex();
1a52e0ed 146
147 while(1)
148 {
149// Generate one event
150// --------------------------------------------------------------------------
151 fSpecn = 0;
152 fSpecp = 0;
153// --------------------------------------------------------------------------
154 fDPMjet->GenerateEvent();
155 fTrials++;
156
157 fDPMjet->ImportParticles(fParticles,"All");
158 if (fLHC) Boost();
159
160 // Temporaneo
161 fGenImpPar = fDPMjet->GetBImpac();
162
163 Int_t np = fParticles->GetEntriesFast();
164 printf("\n **************************************************%d\n",np);
165 Int_t nc=0;
166 if (np==0) continue;
167 Int_t i;
168 Int_t* newPos = new Int_t[np];
169 Int_t* pSelected = new Int_t[np];
170
171 for (i = 0; i<np; i++) {
172 newPos[i] = i;
173 pSelected[i] = 0;
174 }
175
176// Get event vertex
177
f2ba02a4 178 fVertex[0] = origin0[0];
179 fVertex[1] = origin0[1];
180 fVertex[2] = origin0[2];
1a52e0ed 181
182// First select parent particles
183
184 for (i = 0; i<np; i++) {
185 TParticle *iparticle = (TParticle *) fParticles->At(i);
1a52e0ed 186
187// Is this a parent particle ?
188
189 if (Stable(iparticle)) continue;
190
191 Bool_t selected = kTRUE;
192 Bool_t hasSelectedDaughters = kFALSE;
193
194 kf = iparticle->GetPdgCode();
195 if (kf == 92) continue;
196 ks = iparticle->GetStatusCode();
197// No initial state partons
198 if (ks==21) continue;
199
200 if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
201 SelectFlavor(kf);
202 hasSelectedDaughters = DaughtersSelection(iparticle);
203
204// Put particle on the stack if it is either selected or
205// it is the mother of at least one seleted particle
206
207 if (selected || hasSelectedDaughters) {
208 nc++;
209 pSelected[i] = 1;
210 } // selected
211 } // particle loop parents
212
213// Now select the final state particles
214
215
216 for (i=0; i<np; i++) {
217 TParticle *iparticle = (TParticle *) fParticles->At(i);
1a52e0ed 218
219// Is this a final state particle ?
220
221 if (!Stable(iparticle)) continue;
222
223 Bool_t selected = kTRUE;
224 kf = iparticle->GetPdgCode();
225 ks = iparticle->GetStatusCode();
226
227// --------------------------------------------------------------------------
228// Count spectator neutrons and protons (ks == 13, 14)
229 if(ks == 13 || ks == 14){
230 if(kf == kNeutron) fSpecn += 1;
231 if(kf == kProton) fSpecp += 1;
232 }
233// --------------------------------------------------------------------------
234
235 if (!fSelectAll) {
236 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
237 if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
238 }
239
240// Put particle on the stack if selected
241
242 if (selected) {
243 nc++;
244 pSelected[i] = 1;
245 } // selected
246 } // particle loop final state
247
248// Write particles to stack
249
250 for (i = 0; i<np; i++) {
251 TParticle * iparticle = (TParticle *) fParticles->At(i);
252 Bool_t hasMother = (iparticle->GetFirstMother()>=0);
1a52e0ed 253 if (pSelected[i]) {
254
255 kf = iparticle->GetPdgCode();
1a52e0ed 256 ks = iparticle->GetStatusCode();
257
258 p[0] = iparticle->Px();
259 p[1] = iparticle->Py();
260 p[2] = iparticle->Pz();
e67e0fb7 261 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
262 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
263 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
264
1a52e0ed 265 tof = kconv*iparticle->T();
266
267 imo = -1;
268 TParticle* mother = 0;
269 if (hasMother) {
270 imo = iparticle->GetFirstMother();
271 mother = (TParticle *) fParticles->At(imo);
272 imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
273 } // if has mother
274
1a52e0ed 275 Bool_t tFlag = (fTrackIt && (ks == 1));
276
642f15cf 277 PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
1a52e0ed 278 KeepTrack(nt);
279 newPos[i] = nt;
280 } // if selected
281 } // particle loop
282 delete[] newPos;
283 delete[] pSelected;
1a52e0ed 284 if (nc>0) {
285 jev += nc;
286 if (jev >= fNpart || fNpart == -1) {
1a52e0ed 287 break;
288 }
1a52e0ed 289 }
290 } // event loop
291 MakeHeader();
292 SetHighWaterMark(nt);
293}
294
295
296//______________________________________________________________________________
297void AliGenDPMjet::KeepFullEvent()
298{
299 fKeep=1;
300}
301
302
303//______________________________________________________________________________
304/*void AliGenDPMjet::EvaluateCrossSections()
305{
306// Glauber Calculation of geometrical x-section
307//
308 Float_t xTot = 0.; // barn
309 Float_t xTotHard = 0.; // barn
310 Float_t xPart = 0.; // barn
311 Float_t xPartHard = 0.; // barn
312 Float_t sigmaHard = 0.1; // mbarn
313 Float_t bMin = 0.;
314 Float_t bMax = fDPMjet->GetProjRadius()+fDPMjet->GetTargRadius();
315 const Float_t kdib = 0.2;
316 Int_t kMax = Int_t((bMax-bMin)/kdib)+1;
317
318
319 printf("\n Projectile Radius (fm): %f \n",fDPMjet->GetProjRadius());
320 printf("\n Target Radius (fm): %f \n",fDPMjet->GetTargRadius());
321 Int_t i;
322 Float_t oldvalue= 0.;
323
324 Float_t* b = new Float_t[kMax];
325 Float_t* si1 = new Float_t[kMax];
326 Float_t* si2 = new Float_t[kMax];
327
328 for (i = 0; i < kMax; i++)
329 {
330 Float_t xb = bMin+i*kdib;
331 Float_t ov;
332 ov=fDPMjet->Profile(xb);
333 // ATT!->Manca la x-sec anel. nucleone-nucleone
334 Float_t gb = 2.*0.01*fDPMjet->TMath::Pi()*kdib*xb*(1.-TMath::Exp(-fDPMjet->GetXSFrac()*ov));
335 Float_t gbh = 2.*0.01*fDPMjet->TMath::Pi()*kdib*xb*sigmaHard*ov;
336 xTot+=gb;
337 xTotHard += gbh;
338 if (xb > fMinImpactParam && xb < fMaxImpactParam)
339 {
340 xPart += gb;
341 xPartHard += gbh;
342 }
343
344 if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
345 oldvalue = xTot;
346 printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
347 printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
348 if (i>0) {
349 si1[i] = gb/kdib;
350 si2[i] = gbh/gb;
351 b[i] = xb;
352 }
353 }
354
355 printf("\n Total cross section (barn): %f \n",xTot);
356 printf("\n Hard cross section (barn): %f \n \n",xTotHard);
357 printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
358 printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
359
360// Store result as a graph
361 b[0] = 0;
362 si1[0] = 0;
363 si2[0]=si2[1];
364
365 fDsigmaDb = new TGraph(i, b, si1);
366 fDnDb = new TGraph(i, b, si2);
367}*/
368
369
370//______________________________________________________________________________
371Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
372{
373//
374// Looks recursively if one of the daughters has been selected
375//
376// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
377 Int_t imin = -1;
378 Int_t imax = -1;
379 Int_t i;
380 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
381 Bool_t selected = kFALSE;
382 if (hasDaughters) {
383 imin = iparticle->GetFirstDaughter();
384 imax = iparticle->GetLastDaughter();
385 for (i = imin; i <= imax; i++){
386 TParticle * jparticle = (TParticle *) fParticles->At(i);
387 Int_t ip = jparticle->GetPdgCode();
388 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
389 selected=kTRUE; break;
390 }
391 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
392 }
393 } else {
394 return kFALSE;
395 }
396 return selected;
397}
398
399
400
401//______________________________________________________________________________
402Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
403{
404// Select flavor of particle
405// 0: all
406// 4: charm and beauty
407// 5: beauty
408 Bool_t res = 0;
409
410 if (fFlavor == 0) {
411 res = kTRUE;
412 } else {
413 Int_t ifl = TMath::Abs(pid/100);
414 if (ifl > 10) ifl/=10;
415 res = (fFlavor == ifl);
416 }
417//
418// This part if gamma writing is inhibited
419 if (fNoGammas)
420 res = res && (pid != kGamma && pid != kPi0);
421//
422 return res;
423}
424
425//______________________________________________________________________________
426Bool_t AliGenDPMjet::Stable(TParticle* particle)
427{
428// Return true for a stable particle
429//
430
431// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
432 if (particle->GetStatusCode() == 1) return kTRUE;
433 else return kFALSE;
434
435}
436
1a52e0ed 437//______________________________________________________________________________
438void AliGenDPMjet::MakeHeader()
439{
440// Builds the event header, to be called after each event
441 AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
442 ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
443 ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
444 ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
1a52e0ed 445 ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
446 fDPMjet->GetfIt());
1a52e0ed 447
448// Bookkeeping for kinematic bias
449 ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
450// Event Vertex
f2ba02a4 451 header->SetPrimaryVertex(fVertex);
1a52e0ed 452 gAlice->SetGenEventHeader(header);
453}
454
455
456
457//______________________________________________________________________________
ba758f5a 458AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /*rhs*/)
1a52e0ed 459{
460// Assignment operator
461 return *this;
462}
463
464
465
466//______________________________________________________________________________