git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Classed imported from EVGEN.
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
CommitLineData
8d2cd130 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$Log$
18Revision 1.69 2003/01/14 10:50:19 alibrary
19Cleanup of STEER coding conventions
20
21Revision 1.68 2002/12/11 09:16:16 morsch
22Use GetJets to fill header.
23
24Revision 1.67 2002/12/09 15:24:09 morsch
25Same trigger routine can use Pycell or Pyclus.
26
27Revision 1.66 2002/12/09 08:22:56 morsch
28UA1 jet finder (Pycell) for software triggering added.
29
30Revision 1.65 2002/11/19 08:57:10 morsch
31Configuration of pt-kick added.
32
33Revision 1.64 2002/11/15 00:43:06 morsch
34Changes for kPyJets
35- initial and final state g-radiation + pt-kick default
36- trigger based on parton clusters (using pyclus)
37- trigger jets are stored in header.
38
39Revision 1.63 2002/10/14 14:55:35 hristov
40Merging the VirtualMC branch to the main development branch (HEAD)
41
42Revision 1.52.4.4 2002/10/10 16:40:08 hristov
43Updating VirtualMC to v3-09-02
44
45Revision 1.62 2002/09/24 10:00:01 morsch
46CheckTrigger() corrected.
47
48Revision 1.61 2002/07/30 09:52:38 morsch
49Call SetGammaPhiRange() and SetGammaEtaRange() in the constructor.
50
51Revision 1.60 2002/07/19 14:49:03 morsch
52Typo corrected.
53
54Revision 1.59 2002/07/19 14:35:36 morsch
55Count total number of trials. Print mean Q, x1, x2.
56
57Revision 1.58 2002/07/17 10:04:09 morsch
58SetYHard method added.
59
60Revision 1.57 2002/05/22 13:22:53 morsch
61Process kPyMbNonDiffr added.
62
63Revision 1.56 2002/04/26 10:30:01 morsch
64Option kPyBeautyPbMNR added. (N. Carrer).
65
66Revision 1.55 2002/04/17 10:23:56 morsch
67Coding Rule violations corrected.
68
69Revision 1.54 2002/03/28 11:49:10 morsch
70Pass status code in SetTrack.
71
72Revision 1.53 2002/03/25 14:51:13 morsch
73New stack-fill and count options introduced (N. Carrer).
74
75Revision 1.51 2002/03/06 08:46:57 morsch
76- Loop until np-1
77- delete dyn. alloc. arrays (N. Carrer)
78
79Revision 1.50 2002/03/03 13:48:50 morsch
80Option kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
81NLO calculations (Nicola Carrer).
82
83Revision 1.49 2002/02/08 16:50:50 morsch
84Add name and title in constructor.
85
86Revision 1.48 2001/12/20 11:44:28 morsch
87Add kinematic bias for direct gamma production.
88
89Revision 1.47 2001/12/19 14:45:00 morsch
90Store number of trials in header.
91
92Revision 1.46 2001/12/19 10:36:19 morsch
93Add possibility if jet kinematic biasing.
94
95Revision 1.45 2001/11/28 08:06:52 morsch
96Use fMaxLifeTime parameter.
97
98Revision 1.44 2001/11/27 13:13:07 morsch
99Maximum lifetime for long-lived particles to be put on the stack is parameter.
100It can be set via SetMaximumLifetime(..).
101
102Revision 1.43 2001/10/21 18:35:56 hristov
103Several pointers were set to zero in the default constructors to avoid memory management problems
104
105Revision 1.42 2001/10/15 08:21:55 morsch
106Vertex truncation settings moved to AliGenMC.
107
108Revision 1.41 2001/10/08 08:45:42 morsch
109Possibility of vertex cut added.
110
111Revision 1.40 2001/09/25 11:30:23 morsch
112Pass event vertex to header.
113
114Revision 1.39 2001/07/27 17:09:36 morsch
115Use local SetTrack, KeepTrack and SetHighWaterMark methods
116to delegate either to local stack or to stack owned by AliRun.
117(Piotr Skowronski, A.M.)
118
119Revision 1.38 2001/07/13 10:58:54 morsch
120- Some coded moved to AliGenMC
121- Improved handling of secondary vertices.
122
123Revision 1.37 2001/06/28 11:17:28 morsch
124SetEventListRange setter added. Events in specified range are listed for
125debugging. (Yuri Kharlov)
126
127Revision 1.36 2001/03/30 07:05:49 morsch
128Final print-out in finish run.
129Write parton system for jet-production (preliminary solution).
130
131Revision 1.35 2001/03/09 13:03:40 morsch
132Process_t and Struc_Func_t moved to AliPythia.h
133
134Revision 1.34 2001/02/14 15:50:40 hristov
135The last particle in event marked using SetHighWaterMark
136
137Revision 1.33 2001/01/30 09:23:12 hristov
138Streamers removed (R.Brun)
139
140Revision 1.32 2001/01/26 19:55:51 hristov
141Major upgrade of AliRoot code
142
143Revision 1.31 2001/01/17 10:54:31 hristov
144Better protection against FPE
145
146Revision 1.30 2000/12/18 08:55:35 morsch
147Make AliPythia dependent generartors work with new scheme of random number generation
148
149Revision 1.29 2000/12/04 11:22:03 morsch
150Init of sRandom as in 1.15
151
152Revision 1.28 2000/12/02 11:41:39 morsch
153Use SetRandom() to initialize random number generator in constructor.
154
155Revision 1.27 2000/11/30 20:29:02 morsch
156Initialise static variable sRandom in constructor: sRandom = fRandom;
157
158Revision 1.26 2000/11/30 07:12:50 alibrary
159Introducing new Rndm and QA classes
160
161Revision 1.25 2000/10/18 19:11:27 hristov
162Division by zero fixed
163
164Revision 1.24 2000/09/18 10:41:35 morsch
165Add possibility to use nuclear structure functions from PDF library V8.
166
167Revision 1.23 2000/09/14 14:05:40 morsch
168dito
169
170Revision 1.22 2000/09/14 14:02:22 morsch
171- Correct conversion from mm to cm when passing particle vertex to MC.
172- Correct handling of fForceDecay == all.
173
174Revision 1.21 2000/09/12 14:14:55 morsch
175Call fDecayer->ForceDecay() at the beginning of Generate().
176
177Revision 1.20 2000/09/06 14:29:33 morsch
178Use AliPythia for event generation an AliDecayPythia for decays.
179Correct handling of "nodecay" option
180
181Revision 1.19 2000/07/11 18:24:56 fca
182Coding convention corrections + few minor bug fixes
183
184Revision 1.18 2000/06/30 12:40:34 morsch
185Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
186Geant units (cm).
187
188Revision 1.17 2000/06/09 20:34:07 morsch
189All coding rule violations except RS3 corrected
190
191Revision 1.16 2000/05/15 15:04:20 morsch
192The full event is written for fNtrack = -1
193Coding rule violations corrected.
194
195Revision 1.15 2000/04/26 10:14:24 morsch
196Particles array has one entry more than pythia particle list. Upper bound of
197particle loop changed to np-1 (R. Guernane, AM)
198
199Revision 1.14 2000/04/05 08:36:13 morsch
200Check status code of particles in Pythia event
201to avoid double counting as partonic state and final state particle.
202
203Revision 1.13 1999/11/09 07:38:48 fca
204Changes for compatibility with version 2.23 of ROOT
205
206Revision 1.12 1999/11/03 17:43:20 fca
207New version from G.Martinez & A.Morsch
208
209Revision 1.11 1999/09/29 09:24:14 fca
210Introduction of the Copyright and cvs Log
211*/
212
213//
214// Generator using the TPythia interface (via AliPythia)
215// to generate pp collisions.
216// Using SetNuclei() also nuclear modifications to the structure functions
217// can be taken into account. This makes, of course, only sense for the
218// generation of the products of hard processes (heavy flavor, jets ...)
219//
220// andreas.morsch@cern.ch
221//
222
223#include <TDatabasePDG.h>
224#include <TParticle.h>
225#include <TPDGCode.h>
226#include <TSystem.h>
227#include <TTree.h>
228
229#include "AliConst.h"
230#include "AliDecayerPythia.h"
231#include "AliGenPythia.h"
232#include "AliGenPythiaEventHeader.h"
233#include "AliPythia.h"
234#include "AliRun.h"
235
236 ClassImp(AliGenPythia)
237
238AliGenPythia::AliGenPythia()
239 :AliGenMC()
240{
241// Default Constructor
242 fParticles = 0;
243 fPythia = 0;
244 fDecayer = new AliDecayerPythia();
245 SetEventListRange();
246 SetJetPhiRange();
247 SetJetEtaRange();
248 SetJetEtRange();
249 SetGammaPhiRange();
250 SetGammaEtaRange();
251 SetPtKick();
252}
253
254AliGenPythia::AliGenPythia(Int_t npart)
255 :AliGenMC(npart)
256{
257// default charm production at 5. 5 TeV
258// semimuonic decay
259// structure function GRVHO
260//
261 fName = "Pythia";
262 fTitle= "Particle Generator using PYTHIA";
263 fXsection = 0.;
264 fNucA1=0;
265 fNucA2=0;
266 SetProcess();
267 SetStrucFunc();
268 SetForceDecay();
269 SetPtHard();
270 SetYHard();
271 SetEnergyCMS();
272 fDecayer = new AliDecayerPythia();
273 // Set random number generator
274 sRandom=fRandom;
275 fFlavorSelect = 0;
276 // Produced particles
277 fParticles = new TClonesArray("TParticle",1000);
278 fEventVertex.Set(3);
279 SetEventListRange();
280 SetJetPhiRange();
281 SetJetEtaRange();
282 SetJetEtRange();
283 SetGammaPhiRange();
284 SetGammaEtaRange();
285 SetJetReconstructionMode();
286 SetPtKick();
287 // Options determining what to keep in the stack (Heavy flavour generation)
288 fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
289 fFeedDownOpt = kTRUE; // allow feed down from higher family
290 // Fragmentation on/off
291 fFragmentation = kTRUE;
292 // Default counting mode
293 fCountMode = kCountAll;
294}
295
296AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
297{
298// copy constructor
299 Pythia.Copy(*this);
300}
301
302AliGenPythia::~AliGenPythia()
303{
304// Destructor
305}
306
307void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
308{
309 // Set a range of event numbers, for which a table
310 // of generated particle will be printed
311 fDebugEventFirst = eventFirst;
312 fDebugEventLast = eventLast;
313 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
314}
315
316void AliGenPythia::Init()
317{
318// Initialisation
319
320 SetMC(AliPythia::Instance());
321 fPythia=(AliPythia*) fgMCEvGen;
322//
323 fParentWeight=1./Float_t(fNpart);
324//
325// Forward Paramters to the AliPythia object
326 fDecayer->SetForceDecay(fForceDecay);
327 fDecayer->Init();
328
329
330 fPythia->SetCKIN(3,fPtHardMin);
331 fPythia->SetCKIN(4,fPtHardMax);
332 fPythia->SetCKIN(7,fYHardMin);
333 fPythia->SetCKIN(8,fYHardMax);
334
335 if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);
336 // Fragmentation?
337 if (fFragmentation) {
338 fPythia->SetMSTP(111,1);
339 } else {
340 fPythia->SetMSTP(111,0);
341 }
342
343
344// initial state radiation
345 fPythia->SetMSTP(61,fGinit);
346// final state radiation
347 fPythia->SetMSTP(71,fGfinal);
348// pt - kick
349 if (fPtKick > 0.) {
350 fPythia->SetMSTP(91,1);
351 fPythia->SetPARP(91,fPtKick);
352 } else {
353 fPythia->SetMSTP(91,0);
354 }
355
356 // fPythia->SetMSTJ(1,2);
357 //
358 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
359
360// Parent and Children Selection
361 switch (fProcess)
362 {
363 case kPyCharm:
364 case kPyCharmUnforced:
365 case kPyCharmPbMNR:
366 fParentSelect[0] = 411;
367 fParentSelect[1] = 421;
368 fParentSelect[2] = 431;
369 fParentSelect[3] = 4122;
370 fFlavorSelect = 4;
371 break;
372 case kPyD0PbMNR:
373 fParentSelect[0] = 421;
374 fFlavorSelect = 4;
375 break;
376 case kPyBeauty:
377 case kPyBeautyPbMNR:
378 fParentSelect[0]= 511;
379 fParentSelect[1]= 521;
380 fParentSelect[2]= 531;
381 fParentSelect[3]= 5122;
382 fParentSelect[4]= 5132;
383 fParentSelect[5]= 5232;
384 fParentSelect[6]= 5332;
385 fFlavorSelect = 5;
386 break;
387 case kPyBeautyUnforced:
388 fParentSelect[0] = 511;
389 fParentSelect[1] = 521;
390 fParentSelect[2] = 531;
391 fParentSelect[3] = 5122;
392 fParentSelect[4] = 5132;
393 fParentSelect[5] = 5232;
394 fParentSelect[6] = 5332;
395 fFlavorSelect = 5;
396 break;
397 case kPyJpsiChi:
398 case kPyJpsi:
399 fParentSelect[0] = 443;
400 break;
401 case kPyMb:
402 case kPyMbNonDiffr:
403 case kPyJets:
404 case kPyDirectGamma:
405 break;
406 }
407//
408// This counts the total number of calls to Pyevnt() per run.
409 fTrialsRun = 0;
410 fQ = 0.;
411 fX1 = 0.;
412 fX2 = 0.;
413 fNev = 0 ;
414//
415 AliGenMC::Init();
416}
417
418void AliGenPythia::Generate()
419{
420// Generate one event
421 fDecayer->ForceDecay();
422
423 Float_t polar[3] = {0,0,0};
424 Float_t origin[3] = {0,0,0};
425 Float_t p[3];
426// converts from mm/c to s
427 const Float_t kconv=0.001/2.999792458e8;
428//
429 Int_t nt=0;
430 Int_t jev=0;
431 Int_t j, kf;
432 fTrials=0;
433
434 // Set collision vertex position
435 if(fVertexSmear==kPerEvent) {
436 fPythia->SetMSTP(151,1);
437 for (j=0;j<3;j++) {
438 fPythia->SetPARP(151+j, fOsigma[j]*10.);
439 }
440 } else if (fVertexSmear==kPerTrack) {
441 fPythia->SetMSTP(151,0);
442 }
443// event loop
444 while(1)
445 {
446 fPythia->Pyevnt();
447 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
448 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
449 fTrials++;
450
451 fPythia->ImportParticles(fParticles,"All");
452
453//
454//
455//
456 Int_t i;
457
458 Int_t np = fParticles->GetEntriesFast();
459 if (np == 0 ) continue;
460// Get event vertex and discard the event if the Z coord. is too big
461 TParticle *iparticle = (TParticle *) fParticles->At(0);
462 Float_t distz = iparticle->Vz()/10.;
463 if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
464//
465 fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
466 fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
467 fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
468//
469 Int_t* pParent = new Int_t[np];
470 Int_t* pSelected = new Int_t[np];
471 Int_t* trackIt = new Int_t[np];
472 for (i=0; i< np; i++) {
473 pParent[i] = -1;
474 pSelected[i] = 0;
475 trackIt[i] = 0;
476 }
477
478 Int_t nc = 0; // Total n. of selected particles
479 Int_t nParents = 0; // Selected parents
480 Int_t nTkbles = 0; // Trackable particles
481 if (fProcess != kPyMb && fProcess != kPyJets &&
482 fProcess != kPyDirectGamma &&
483 fProcess != kPyMbNonDiffr) {
484
485 for (i = 0; i<np; i++) {
486 iparticle = (TParticle *) fParticles->At(i);
487 Int_t ks = iparticle->GetStatusCode();
488 kf = CheckPDGCode(iparticle->GetPdgCode());
489// No initial state partons
490 if (ks==21) continue;
491//
492// Heavy Flavor Selection
493//
494 // quark ?
495 kf = TMath::Abs(kf);
496 Int_t kfl = kf;
497 // meson ?
498 if (kfl > 10) kfl/=100;
499 // baryon
500 if (kfl > 10) kfl/=10;
501 if (kfl > 10) kfl/=10;
502
503 Int_t ipa = iparticle->GetFirstMother()-1;
504 Int_t kfMo = 0;
505
506 if (ipa > -1) {
507 TParticle * mother = (TParticle *) fParticles->At(ipa);
508 kfMo = TMath::Abs(mother->GetPdgCode());
509 }
510 // What to keep in Stack?
511 Bool_t flavorOK = kFALSE;
512 Bool_t selectOK = kFALSE;
513 if (fFeedDownOpt) {
514 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
515 } else {
516 if (kfl > fFlavorSelect) {
517 nc = -1;
518 break;
519 }
520 if (kfl == fFlavorSelect) flavorOK = kTRUE;
521 }
522 switch (fStackFillOpt) {
523 case kFlavorSelection:
524 selectOK = kTRUE;
525 break;
526 case kParentSelection:
527 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
528 break;
529 }
530 if (flavorOK && selectOK) {
531//
532// Heavy flavor hadron or quark
533//
534// Kinematic seletion on final state heavy flavor mesons
535 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
536 {
537 continue;
538 }
539 pSelected[i] = 1;
540 if (ParentSelected(kf)) ++nParents; // Update parent count
541// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
542 } else {
543// Kinematic seletion on decay products
544 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
545 && !KinematicSelection(iparticle, 1))
546 {
547 continue;
548 }
549//
550// Decay products
551// Select if mother was selected and is not tracked
552
553 if (pSelected[ipa] &&
554 !trackIt[ipa] && // mother will be tracked ?
555 kfMo != 5 && // mother is b-quark, don't store fragments
556 kfMo != 4 && // mother is c-quark, don't store fragments
557 kf != 92) // don't store string
558 {
559//
560// Semi-stable or de-selected: diselect decay products:
561//
562//
563 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
564 {
565 Int_t ipF = iparticle->GetFirstDaughter();
566 Int_t ipL = iparticle->GetLastDaughter();
567 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
568 }
569// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
570 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
571 }
572 }
573 if (pSelected[i] == -1) pSelected[i] = 0;
574 if (!pSelected[i]) continue;
575 // Count quarks only if you did not include fragmentation
576 if (fFragmentation && kf <= 10) continue;
577 nc++;
578// Decision on tracking
579 trackIt[i] = 0;
580//
581// Track final state particle
582 if (ks == 1) trackIt[i] = 1;
583// Track semi-stable particles
584 if ((ks ==1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
585// Track particles selected by process if undecayed.
586 if (fForceDecay == kNoDecay) {
587 if (ParentSelected(kf)) trackIt[i] = 1;
588 } else {
589 if (ParentSelected(kf)) trackIt[i] = 0;
590 }
591 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
592//
593//
594
595 } // particle selection loop
596 if (nc > 0) {
597 for (i = 0; i<np; i++) {
598 if (!pSelected[i]) continue;
599 TParticle * iparticle = (TParticle *) fParticles->At(i);
600 kf = CheckPDGCode(iparticle->GetPdgCode());
601 Int_t ks = iparticle->GetStatusCode();
602 p[0] = iparticle->Px();
603 p[1] = iparticle->Py();
604 p[2] = iparticle->Pz();
605 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
606 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
607 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
608 Float_t tof = kconv*iparticle->T();
609 Int_t ipa = iparticle->GetFirstMother()-1;
610 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
611 SetTrack(fTrackIt*trackIt[i] ,
612 iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
613 pParent[i] = nt;
614 KeepTrack(nt);
615 } // SetTrack loop
616 }
617 } else {
618 nc = GenerateMB();
619 } // mb ?
620
621 if (pParent) delete[] pParent;
622 if (pSelected) delete[] pSelected;
623 if (trackIt) delete[] trackIt;
624
625 if (nc > 0) {
626 switch (fCountMode) {
627 case kCountAll:
628 // printf(" Count all \n");
629 jev += nc;
630 break;
631 case kCountParents:
632 // printf(" Count parents \n");
633 jev += nParents;
634 break;
635 case kCountTrackables:
636 // printf(" Count trackable \n");
637 jev += nTkbles;
638 break;
639 }
640 if (jev >= fNpart || fNpart == -1) {
641 fKineBias=Float_t(fNpart)/Float_t(fTrials);
642 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
643
644 fQ += fPythia->GetVINT(51);
645 fX1 += fPythia->GetVINT(41);
646 fX2 += fPythia->GetVINT(42);
647 fTrialsRun += fTrials;
648 fNev++;
649 MakeHeader();
650 break;
651 }
652 }
653 } // event loop
654 SetHighWaterMark(nt);
655// adjust weight due to kinematic selection
656 AdjustWeights();
657// get cross-section
658 fXsection=fPythia->GetPARI(1);
659}
660
661Int_t AliGenPythia::GenerateMB()
662{
663//
664// Min Bias selection and other global selections
665//
666 Int_t i, kf, nt, iparent;
667 Int_t nc = 0;
668 Float_t p[3];
669 Float_t polar[3] = {0,0,0};
670 Float_t origin[3] = {0,0,0};
671// converts from mm/c to s
672 const Float_t kconv=0.001/2.999792458e8;
673
674 Int_t np = fParticles->GetEntriesFast();
675 Int_t* pParent = new Int_t[np];
676 for (i=0; i< np; i++) pParent[i] = -1;
677 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
678 TParticle* jet1 = (TParticle *) fParticles->At(6);
679 TParticle* jet2 = (TParticle *) fParticles->At(7);
680 if (!CheckTrigger(jet1, jet2)) return 0;
681 }
682
683 for (i = 0; i<np; i++) {
684 Int_t trackIt = 0;
685 TParticle * iparticle = (TParticle *) fParticles->At(i);
686 kf = CheckPDGCode(iparticle->GetPdgCode());
687 Int_t ks = iparticle->GetStatusCode();
688 Int_t km = iparticle->GetFirstMother();
689 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
690 (ks != 1) ||
691 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
692 nc++;
693 if (ks == 1) trackIt = 1;
694 Int_t ipa = iparticle->GetFirstMother()-1;
695
696 iparent = (ipa > -1) ? pParent[ipa] : -1;
697
698//
699// store track information
700 p[0] = iparticle->Px();
701 p[1] = iparticle->Py();
702 p[2] = iparticle->Pz();
703 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
704 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
705 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
706 Float_t tof=kconv*iparticle->T();
707 SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
708 tof, kPPrimary, nt, 1., ks);
709 KeepTrack(nt);
710 pParent[i] = nt;
711 } // select particle
712 } // particle loop
713
714 if (pParent) delete[] pParent;
715
716 printf("\n I've put %i particles on the stack \n",nc);
717 return nc;
718}
719
720
721void AliGenPythia::FinishRun()
722{
723// Print x-section summary
724 fPythia->Pystat(1);
725 fQ /= fNev;
726 fX1 /= fNev;
727 fX2 /= fNev;
728 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
729 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
730
731
732}
733
734void AliGenPythia::AdjustWeights()
735{
736// Adjust the weights after generation of all events
737//
738 TParticle *part;
739 Int_t ntrack=gAlice->GetNtrack();
740 for (Int_t i=0; i<ntrack; i++) {
741 part= gAlice->Particle(i);
742 part->SetWeight(part->GetWeight()*fKineBias);
743 }
744}
745
746void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
747{
748// Treat protons as inside nuclei with mass numbers a1 and a2
749 fNucA1 = a1;
750 fNucA2 = a2;
751}
752
753
754void AliGenPythia::MakeHeader()
755{
756// Builds the event header, to be called after each event
757 AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
758//
759// Event type
760 ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
761//
762// Number of trials
763 ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials);
764//
765// Event Vertex
766 header->SetPrimaryVertex(fEventVertex);
767//
768// Jets that have triggered
769 if (fProcess == kPyJets)
770 {
771 Int_t ntrig, njet;
772 Float_t jets[4][10];
773 GetJets(njet, ntrig, jets);
774
775 for (Int_t i = 0; i < ntrig; i++) {
776 ((AliGenPythiaEventHeader*) header)->AddJet(jets[0][i], jets[1][i], jets[2][i],
777 jets[3][i]);
778 }
779 }
780 gAlice->SetGenEventHeader(header);
781}
782
783
784Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
785{
786// Check the kinematic trigger condition
787//
788 Double_t eta[2];
789 eta[0] = jet1->Eta();
790 eta[1] = jet2->Eta();
791 Double_t phi[2];
792 phi[0] = jet1->Phi();
793 phi[1] = jet2->Phi();
794 Int_t pdg[2];
795 pdg[0] = jet1->GetPdgCode();
796 pdg[1] = jet2->GetPdgCode();
797 Bool_t triggered = kFALSE;
798
799 if (fProcess == kPyJets) {
800 Int_t njets = 0;
801 Int_t ntrig = 0;
802 Float_t jets[4][10];
803//
804// Use Pythia clustering on parton level to determine jet axis
805//
806 GetJets(njets, ntrig, jets);
807
808 if (ntrig) triggered = kTRUE;
809//
810 } else {
811 Int_t ij = 0;
812 Int_t ig = 1;
813 if (pdg[0] == kGamma) {
814 ij = 1;
815 ig = 0;
816 }
817 //Check eta range first...
818 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
819 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
820 {
821 //Eta is okay, now check phi range
822 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
823 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
824 {
825 triggered = kTRUE;
826 }
827 }
828 }
829 return triggered;
830}
831
832AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
833{
834// Assignment operator
835 return *this;
836}
837
838void AliGenPythia::LoadEvent()
839{
840//
841// Load event into Pythia Common Block
842//
843
844
845 Int_t npart = (Int_t) (gAlice->TreeK())->GetEntries();
846 (fPythia->GetPyjets())->N = npart;
847
848 for (Int_t part = 0; part < npart; part++) {
849 TParticle *MPart = gAlice->Particle(part);
850 Int_t kf = MPart->GetPdgCode();
851 Int_t ks = MPart->GetStatusCode();
852 Float_t px = MPart->Px();
853 Float_t py = MPart->Py();
854 Float_t pz = MPart->Pz();
855 Float_t e = MPart->Energy();
856 Float_t p = TMath::Sqrt(px * px + py * py + pz * pz);
857 Float_t m = TMath::Sqrt(e * e - p * p);
858
859
860 (fPythia->GetPyjets())->P[0][part] = px;
861 (fPythia->GetPyjets())->P[1][part] = py;
862 (fPythia->GetPyjets())->P[2][part] = pz;
863 (fPythia->GetPyjets())->P[3][part] = e;
864 (fPythia->GetPyjets())->P[4][part] = m;
865
866 (fPythia->GetPyjets())->K[1][part] = kf;
867 (fPythia->GetPyjets())->K[0][part] = ks;
868 }
869}
870
871void AliGenPythia::RecJetsUA1(Float_t eCellMin, Float_t eCellSeed, Float_t eMin, Float_t rMax,
872 Int_t& njets, Float_t jets [4][50])
873{
874//
875// Calls the Pythia jet finding algorithm to find jets in the current event
876//
877//
878// Configure detector (EMCAL like)
879//
880 fPythia->SetPARU(51,2.);
881 fPythia->SetMSTU(51,Int_t(96 * 2./0.7));
882 fPythia->SetMSTU(52,3 * 144);
883//
884// Configure Jet Finder
885//
886 fPythia->SetPARU(58, eCellMin);
887 fPythia->SetPARU(52, eCellSeed);
888 fPythia->SetPARU(53, eMin);
889 fPythia->SetPARU(54, rMax);
890 fPythia->SetMSTU(54, 2);
891//
892// Save jets
893 Int_t n = fPythia->GetN();
894
895//
896// Run Jet Finder
897 fPythia->Pycell(njets);
898 Int_t i;
899 for (i = 0; i < njets; i++) {
900 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
901 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
902 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
903 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
904
905 jets[0][i] = px;
906 jets[1][i] = py;
907 jets[2][i] = pz;
908 jets[3][i] = e;
909 }
910}
911
912
913
914void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
915{
916//
917// Calls the Pythia clustering algorithm to find jets in the current event
918//
919 Int_t n = fPythia->GetN();
920 nJets = 0;
921 nJetsTrig = 0;
922 if (fJetReconstruction == kCluster) {
923//
924// Configure cluster algorithm
925//
926 fPythia->SetPARU(43, 2.);
927 fPythia->SetMSTU(41, 1);
928//
929// Call cluster algorithm
930//
931 fPythia->Pyclus(nJets);
932//
933// Loading jets from common block
934//
935 } else {
936//
937// Configure detector (EMCAL like)
938//
939 fPythia->SetPARU(51,2.);
940 fPythia->SetMSTU(51,Int_t(96 * 2./0.7));
941 fPythia->SetMSTU(52,3 * 144);
942//
943// Configure Jet Finder
944//
945 fPythia->SetPARU(58, 0.0);
946 fPythia->SetPARU(52, 4.0);
947 fPythia->SetPARU(53, 10.0);
948 fPythia->SetPARU(54, 1.0);
949 fPythia->SetMSTU(54, 2);
950//
951// Run Jet Finder
952 fPythia->Pycell(nJets);
953 }
954
955 Int_t i;
956 for (i = 0; i < nJets; i++) {
957 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
958 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
959 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
960 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
961 Float_t pt = TMath::Sqrt(px * px + py * py);
962 Float_t phi = TMath::ATan2(py,px);
963 Float_t theta = TMath::ATan2(pt,pz);
964 Float_t et = e * TMath::Sin(theta);
965 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
966
967 if (
968 eta > fEtaMinJet && eta < fEtaMaxJet &&
969 phi > fPhiMinJet && eta < fPhiMaxJet &&
970 et > fEtMinJet && et < fEtMaxJet
971 )
972 {
973 jets[0][nJetsTrig] = px;
974 jets[1][nJetsTrig] = py;
975 jets[2][nJetsTrig] = pz;
976 jets[3][nJetsTrig] = e;
977 nJetsTrig++;
978
979 } else {
980// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
981 }
982 }
983}
984
985
986#ifdef never
987void AliGenPythia::Streamer(TBuffer &R__b)
988{
989 // Stream an object of class AliGenPythia.
990
991 if (R__b.IsReading()) {
992 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
993 AliGenerator::Streamer(R__b);
994 R__b >> (Int_t&)fProcess;
995 R__b >> (Int_t&)fStrucFunc;
996 R__b >> (Int_t&)fForceDecay;
997 R__b >> fEnergyCMS;
998 R__b >> fKineBias;
999 R__b >> fTrials;
1000 fParentSelect.Streamer(R__b);
1001 fChildSelect.Streamer(R__b);
1002 R__b >> fXsection;
1003// (AliPythia::Instance())->Streamer(R__b);
1004 R__b >> fPtHardMin;
1005 R__b >> fPtHardMax;
1006// if (fDecayer) fDecayer->Streamer(R__b);
1007 } else {
1008 R__b.WriteVersion(AliGenPythia::IsA());
1009 AliGenerator::Streamer(R__b);
1010 R__b << (Int_t)fProcess;
1011 R__b << (Int_t)fStrucFunc;
1012 R__b << (Int_t)fForceDecay;
1013 R__b << fEnergyCMS;
1014 R__b << fKineBias;
1015 R__b << fTrials;
1016 fParentSelect.Streamer(R__b);
1017 fChildSelect.Streamer(R__b);
1018 R__b << fXsection;
1019// R__b << fPythia;
1020 R__b << fPtHardMin;
1021 R__b << fPtHardMax;
1022 // fDecayer->Streamer(R__b);
1023 }
1024}
1025#endif
1026