]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PYTHIA6/AliGenPythia.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
kPyLhwgMb added.
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
... / ...
CommitLineData
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/* $Id$ */
17
18//
19// Generator using the TPythia interface (via AliPythia)
20// to generate pp collisions.
21// Using SetNuclei() also nuclear modifications to the structure functions
22// can be taken into account. This makes, of course, only sense for the
23// generation of the products of hard processes (heavy flavor, jets ...)
24//
25// andreas.morsch@cern.ch
26//
27
28#include <TDatabasePDG.h>
29#include <TParticle.h>
30#include <TPDGCode.h>
31#include <TSystem.h>
32#include <TTree.h>
33#include "AliConst.h"
34#include "AliDecayerPythia.h"
35#include "AliGenPythia.h"
36#include "AliHeader.h"
37#include "AliGenPythiaEventHeader.h"
38#include "AliPythia.h"
39#include "AliPythiaRndm.h"
40#include "AliRun.h"
41#include "AliStack.h"
42#include "AliRunLoader.h"
43#include "AliMC.h"
44#include "pyquenCommon.h"
45
46ClassImp(AliGenPythia)
47
48
49AliGenPythia::AliGenPythia():
50 AliGenMC(),
51 fProcess(kPyCharm),
52 fStrucFunc(kCTEQ5L),
53 fEnergyCMS(5500.),
54 fKineBias(0.),
55 fTrials(0),
56 fTrialsRun(0),
57 fQ(0.),
58 fX1(0.),
59 fX2(0.),
60 fEventTime(0.),
61 fInteractionRate(0.),
62 fTimeWindow(0.),
63 fCurSubEvent(0),
64 fEventsTime(0),
65 fNev(0),
66 fFlavorSelect(0),
67 fXsection(0.),
68 fPythia(0),
69 fPtHardMin(0.),
70 fPtHardMax(1.e4),
71 fYHardMin(-1.e10),
72 fYHardMax(1.e10),
73 fGinit(1),
74 fGfinal(1),
75 fHadronisation(1),
76 fNpartons(0),
77 fReadFromFile(0),
78 fQuench(0),
79 fPtKick(1.),
80 fFullEvent(kTRUE),
81 fDecayer(new AliDecayerPythia()),
82 fDebugEventFirst(-1),
83 fDebugEventLast(-1),
84 fEtMinJet(0.),
85 fEtMaxJet(1.e4),
86 fEtaMinJet(-20.),
87 fEtaMaxJet(20.),
88 fPhiMinJet(0.),
89 fPhiMaxJet(2.* TMath::Pi()),
90 fJetReconstruction(kCell),
91 fEtaMinGamma(-20.),
92 fEtaMaxGamma(20.),
93 fPhiMinGamma(0.),
94 fPhiMaxGamma(2. * TMath::Pi()),
95 fPycellEtaMax(2.),
96 fPycellNEta(274),
97 fPycellNPhi(432),
98 fPycellThreshold(0.),
99 fPycellEtSeed(4.),
100 fPycellMinEtJet(10.),
101 fPycellMaxRadius(1.),
102 fStackFillOpt(kFlavorSelection),
103 fFeedDownOpt(kTRUE),
104 fFragmentation(kTRUE),
105 fSetNuclei(kFALSE),
106 fNewMIS(kFALSE),
107 fHFoff(kFALSE),
108 fTriggerParticle(0),
109 fTriggerEta(0.9),
110 fCountMode(kCountAll),
111 fHeader(0),
112 fRL(0),
113 fFileName(0)
114{
115// Default Constructor
116 SetNuclei(0,0);
117 if (!AliPythiaRndm::GetPythiaRandom())
118 AliPythiaRndm::SetPythiaRandom(GetRandom());
119}
120
121AliGenPythia::AliGenPythia(Int_t npart)
122 :AliGenMC(npart),
123 fProcess(kPyCharm),
124 fStrucFunc(kCTEQ5L),
125 fEnergyCMS(5500.),
126 fKineBias(0.),
127 fTrials(0),
128 fTrialsRun(0),
129 fQ(0.),
130 fX1(0.),
131 fX2(0.),
132 fEventTime(0.),
133 fInteractionRate(0.),
134 fTimeWindow(0.),
135 fCurSubEvent(0),
136 fEventsTime(0),
137 fNev(0),
138 fFlavorSelect(0),
139 fXsection(0.),
140 fPythia(0),
141 fPtHardMin(0.),
142 fPtHardMax(1.e4),
143 fYHardMin(-1.e10),
144 fYHardMax(1.e10),
145 fGinit(kTRUE),
146 fGfinal(kTRUE),
147 fHadronisation(kTRUE),
148 fNpartons(0),
149 fReadFromFile(kFALSE),
150 fQuench(kFALSE),
151 fPtKick(1.),
152 fFullEvent(kTRUE),
153 fDecayer(new AliDecayerPythia()),
154 fDebugEventFirst(-1),
155 fDebugEventLast(-1),
156 fEtMinJet(0.),
157 fEtMaxJet(1.e4),
158 fEtaMinJet(-20.),
159 fEtaMaxJet(20.),
160 fPhiMinJet(0.),
161 fPhiMaxJet(2.* TMath::Pi()),
162 fJetReconstruction(kCell),
163 fEtaMinGamma(-20.),
164 fEtaMaxGamma(20.),
165 fPhiMinGamma(0.),
166 fPhiMaxGamma(2. * TMath::Pi()),
167 fPycellEtaMax(2.),
168 fPycellNEta(274),
169 fPycellNPhi(432),
170 fPycellThreshold(0.),
171 fPycellEtSeed(4.),
172 fPycellMinEtJet(10.),
173 fPycellMaxRadius(1.),
174 fStackFillOpt(kFlavorSelection),
175 fFeedDownOpt(kTRUE),
176 fFragmentation(kTRUE),
177 fSetNuclei(kFALSE),
178 fNewMIS(kFALSE),
179 fHFoff(kFALSE),
180 fTriggerParticle(0),
181 fTriggerEta(0.9),
182 fCountMode(kCountAll),
183 fHeader(0),
184 fRL(0),
185 fFileName(0)
186{
187// default charm production at 5. 5 TeV
188// semimuonic decay
189// structure function GRVHO
190//
191 fName = "Pythia";
192 fTitle= "Particle Generator using PYTHIA";
193 SetForceDecay();
194 // Set random number generator
195 if (!AliPythiaRndm::GetPythiaRandom())
196 AliPythiaRndm::SetPythiaRandom(GetRandom());
197 fParticles = new TClonesArray("TParticle",1000);
198 SetNuclei(0,0);
199 }
200
201AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
202 :AliGenMC(Pythia),
203 fProcess(kPyCharm),
204 fStrucFunc(kCTEQ5L),
205 fEnergyCMS(5500.),
206 fKineBias(0.),
207 fTrials(0),
208 fTrialsRun(0),
209 fQ(0.),
210 fX1(0.),
211 fX2(0.),
212 fEventTime(0.),
213 fInteractionRate(0.),
214 fTimeWindow(0.),
215 fCurSubEvent(0),
216 fEventsTime(0),
217 fNev(0),
218 fFlavorSelect(0),
219 fXsection(0.),
220 fPythia(0),
221 fPtHardMin(0.),
222 fPtHardMax(1.e4),
223 fYHardMin(-1.e10),
224 fYHardMax(1.e10),
225 fGinit(kTRUE),
226 fGfinal(kTRUE),
227 fHadronisation(kTRUE),
228 fNpartons(0),
229 fReadFromFile(kFALSE),
230 fQuench(kFALSE),
231 fPtKick(1.),
232 fFullEvent(kTRUE),
233 fDecayer(new AliDecayerPythia()),
234 fDebugEventFirst(-1),
235 fDebugEventLast(-1),
236 fEtMinJet(0.),
237 fEtMaxJet(1.e4),
238 fEtaMinJet(-20.),
239 fEtaMaxJet(20.),
240 fPhiMinJet(0.),
241 fPhiMaxJet(2.* TMath::Pi()),
242 fJetReconstruction(kCell),
243 fEtaMinGamma(-20.),
244 fEtaMaxGamma(20.),
245 fPhiMinGamma(0.),
246 fPhiMaxGamma(2. * TMath::Pi()),
247 fPycellEtaMax(2.),
248 fPycellNEta(274),
249 fPycellNPhi(432),
250 fPycellThreshold(0.),
251 fPycellEtSeed(4.),
252 fPycellMinEtJet(10.),
253 fPycellMaxRadius(1.),
254 fStackFillOpt(kFlavorSelection),
255 fFeedDownOpt(kTRUE),
256 fFragmentation(kTRUE),
257 fSetNuclei(kFALSE),
258 fNewMIS(kFALSE),
259 fHFoff(kFALSE),
260 fTriggerParticle(0),
261 fTriggerEta(0.9),
262 fCountMode(kCountAll),
263 fHeader(0),
264 fRL(0),
265 fFileName(0)
266{
267// copy constructor
268 Pythia.Copy(*this);
269}
270
271AliGenPythia::~AliGenPythia()
272{
273// Destructor
274 if(fEventsTime) delete fEventsTime;
275}
276
277void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
278{
279// Generate pileup using user specified rate
280 fInteractionRate = rate;
281 fTimeWindow = timewindow;
282 GeneratePileup();
283}
284
285void AliGenPythia::GeneratePileup()
286{
287// Generate sub events time for pileup
288 fEventsTime = 0;
289 if(fInteractionRate == 0.) {
290 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
291 return;
292 }
293
294 Int_t npart = NumberParticles();
295 if(npart < 0) {
296 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
297 return;
298 }
299
300 if(fEventsTime) delete fEventsTime;
301 fEventsTime = new TArrayF(npart);
302 TArrayF &array = *fEventsTime;
303 for(Int_t ipart = 0; ipart < npart; ipart++)
304 array[ipart] = 0.;
305
306 Float_t eventtime = 0.;
307 while(1)
308 {
309 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
310 if(eventtime > fTimeWindow) break;
311 array.Set(array.GetSize()+1);
312 array[array.GetSize()-1] = eventtime;
313 }
314
315 eventtime = 0.;
316 while(1)
317 {
318 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
319 if(TMath::Abs(eventtime) > fTimeWindow) break;
320 array.Set(array.GetSize()+1);
321 array[array.GetSize()-1] = eventtime;
322 }
323
324 SetNumberParticles(fEventsTime->GetSize());
325}
326
327void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
328 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
329{
330// Set pycell parameters
331 fPycellEtaMax = etamax;
332 fPycellNEta = neta;
333 fPycellNPhi = nphi;
334 fPycellThreshold = thresh;
335 fPycellEtSeed = etseed;
336 fPycellMinEtJet = minet;
337 fPycellMaxRadius = r;
338}
339
340
341
342void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
343{
344 // Set a range of event numbers, for which a table
345 // of generated particle will be printed
346 fDebugEventFirst = eventFirst;
347 fDebugEventLast = eventLast;
348 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
349}
350
351void AliGenPythia::Init()
352{
353// Initialisation
354
355 SetMC(AliPythia::Instance());
356 fPythia=(AliPythia*) fMCEvGen;
357
358//
359 fParentWeight=1./Float_t(fNpart);
360//
361// Forward Paramters to the AliPythia object
362 fDecayer->SetForceDecay(fForceDecay);
363 fDecayer->Init();
364
365
366 fPythia->SetCKIN(3,fPtHardMin);
367 fPythia->SetCKIN(4,fPtHardMax);
368 fPythia->SetCKIN(7,fYHardMin);
369 fPythia->SetCKIN(8,fYHardMax);
370
371 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
372 // Fragmentation?
373 if (fFragmentation) {
374 fPythia->SetMSTP(111,1);
375 } else {
376 fPythia->SetMSTP(111,0);
377 }
378
379
380// initial state radiation
381 fPythia->SetMSTP(61,fGinit);
382// final state radiation
383 fPythia->SetMSTP(71,fGfinal);
384// pt - kick
385 if (fPtKick > 0.) {
386 fPythia->SetMSTP(91,1);
387 fPythia->SetPARP(91,fPtKick);
388 } else {
389 fPythia->SetMSTP(91,0);
390 }
391
392
393 if (fReadFromFile) {
394 fRL = AliRunLoader::Open(fFileName, "Partons");
395 fRL->LoadKinematics();
396 fRL->LoadHeader();
397 } else {
398 fRL = 0x0;
399 }
400// Switch off Heavy Flavors on request
401 if (fHFoff) {
402 fPythia->SetMSTP(58, 3);
403 fPythia->SetMSTJ(45, 3);
404 for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
405 }
406 //
407 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
408
409// Parent and Children Selection
410 switch (fProcess)
411 {
412 case kPyOldUEQ2ordered:
413 case kPyOldUEQ2ordered2:
414 case kPyOldPopcorn:
415 break;
416 case kPyCharm:
417 case kPyCharmUnforced:
418 case kPyCharmPbPbMNR:
419 case kPyCharmpPbMNR:
420 case kPyCharmppMNR:
421 case kPyCharmppMNRwmi:
422 fParentSelect[0] = 411;
423 fParentSelect[1] = 421;
424 fParentSelect[2] = 431;
425 fParentSelect[3] = 4122;
426 fFlavorSelect = 4;
427 break;
428 case kPyD0PbPbMNR:
429 case kPyD0pPbMNR:
430 case kPyD0ppMNR:
431 fParentSelect[0] = 421;
432 fFlavorSelect = 4;
433 break;
434 case kPyDPlusPbPbMNR:
435 case kPyDPluspPbMNR:
436 case kPyDPlusppMNR:
437 fParentSelect[0] = 411;
438 fFlavorSelect = 4;
439 break;
440 case kPyDPlusStrangePbPbMNR:
441 case kPyDPlusStrangepPbMNR:
442 case kPyDPlusStrangeppMNR:
443 fParentSelect[0] = 431;
444 fFlavorSelect = 4;
445 break;
446 case kPyBeauty:
447 case kPyBeautyPbPbMNR:
448 case kPyBeautypPbMNR:
449 case kPyBeautyppMNR:
450 case kPyBeautyppMNRwmi:
451 fParentSelect[0]= 511;
452 fParentSelect[1]= 521;
453 fParentSelect[2]= 531;
454 fParentSelect[3]= 5122;
455 fParentSelect[4]= 5132;
456 fParentSelect[5]= 5232;
457 fParentSelect[6]= 5332;
458 fFlavorSelect = 5;
459 break;
460 case kPyBeautyUnforced:
461 fParentSelect[0] = 511;
462 fParentSelect[1] = 521;
463 fParentSelect[2] = 531;
464 fParentSelect[3] = 5122;
465 fParentSelect[4] = 5132;
466 fParentSelect[5] = 5232;
467 fParentSelect[6] = 5332;
468 fFlavorSelect = 5;
469 break;
470 case kPyJpsiChi:
471 case kPyJpsi:
472 fParentSelect[0] = 443;
473 break;
474 case kPyMb:
475 case kPyMbNonDiffr:
476 case kPyMbMSEL1:
477 case kPyJets:
478 case kPyDirectGamma:
479 case kPyLhwgMb:
480 break;
481 case kPyW:
482 case kPyZ:
483 break;
484 }
485//
486//
487// JetFinder for Trigger
488//
489// Configure detector (EMCAL like)
490//
491 fPythia->SetPARU(51, fPycellEtaMax);
492 fPythia->SetMSTU(51, fPycellNEta);
493 fPythia->SetMSTU(52, fPycellNPhi);
494//
495// Configure Jet Finder
496//
497 fPythia->SetPARU(58, fPycellThreshold);
498 fPythia->SetPARU(52, fPycellEtSeed);
499 fPythia->SetPARU(53, fPycellMinEtJet);
500 fPythia->SetPARU(54, fPycellMaxRadius);
501 fPythia->SetMSTU(54, 2);
502//
503// This counts the total number of calls to Pyevnt() per run.
504 fTrialsRun = 0;
505 fQ = 0.;
506 fX1 = 0.;
507 fX2 = 0.;
508 fNev = 0 ;
509//
510//
511//
512 AliGenMC::Init();
513//
514//
515//
516 if (fSetNuclei) {
517 fDyBoost = 0;
518 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
519 }
520
521 if (fQuench) {
522 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
523 }
524 fPythia->SetPARJ(200, 0.0);
525
526 if (fQuench == 3) {
527 // Nestor's change of the splittings
528 fPythia->SetPARJ(200, 0.8);
529 fPythia->SetMSTJ(41, 1); // QCD radiation only
530 fPythia->SetMSTJ(42, 2); // angular ordering
531 fPythia->SetMSTJ(44, 2); // option to run alpha_s
532 fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element
533 fPythia->SetMSTJ(50, 0); // No coherence in first branching
534 fPythia->SetPARJ(82, 1.); // Cut off for parton showers
535 }
536}
537
538void AliGenPythia::Generate()
539{
540// Generate one event
541
542 fDecayer->ForceDecay();
543
544 Float_t polar[3] = {0,0,0};
545 Float_t origin[3] = {0,0,0};
546 Float_t p[4];
547// converts from mm/c to s
548 const Float_t kconv=0.001/2.999792458e8;
549//
550 Int_t nt=0;
551 Int_t jev=0;
552 Int_t j, kf;
553 fTrials=0;
554 fEventTime = 0.;
555
556
557
558 // Set collision vertex position
559 if (fVertexSmear == kPerEvent) Vertex();
560
561// event loop
562 while(1)
563 {
564//
565// Produce event
566//
567//
568// Switch hadronisation off
569//
570 fPythia->SetMSTJ(1, 0);
571//
572// Either produce new event or read partons from file
573//
574 if (!fReadFromFile) {
575 if (!fNewMIS) {
576 fPythia->Pyevnt();
577 } else {
578 fPythia->Pyevnw();
579 }
580 fNpartons = fPythia->GetN();
581 } else {
582 printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
583 fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
584 fPythia->SetN(0);
585 LoadEvent(fRL->Stack(), 0 , 1);
586 fPythia->Pyedit(21);
587 }
588
589//
590// Run quenching routine
591//
592 if (fQuench == 1) {
593 fPythia->Quench();
594 } else if (fQuench == 2){
595 fPythia->Pyquen(208., 0, 0.);
596 } else if (fQuench == 3) {
597 // Quenching is via multiplicative correction of the splittings
598 }
599
600//
601// Switch hadronisation on
602//
603 fPythia->SetMSTJ(1, 1);
604//
605// .. and perform hadronisation
606// printf("Calling hadronisation %d\n", fPythia->GetN());
607 fPythia->Pyexec();
608 fTrials++;
609 fPythia->ImportParticles(fParticles,"All");
610 Boost();
611//
612//
613//
614 Int_t i;
615
616
617 Int_t np = fParticles->GetEntriesFast();
618
619 if (np == 0) continue;
620//
621
622//
623 Int_t* pParent = new Int_t[np];
624 Int_t* pSelected = new Int_t[np];
625 Int_t* trackIt = new Int_t[np];
626 for (i = 0; i < np; i++) {
627 pParent[i] = -1;
628 pSelected[i] = 0;
629 trackIt[i] = 0;
630 }
631
632 Int_t nc = 0; // Total n. of selected particles
633 Int_t nParents = 0; // Selected parents
634 Int_t nTkbles = 0; // Trackable particles
635 if (fProcess != kPyMb && fProcess != kPyJets &&
636 fProcess != kPyDirectGamma &&
637 fProcess != kPyMbNonDiffr &&
638 fProcess != kPyMbMSEL1 &&
639 fProcess != kPyW && fProcess != kPyZ &&
640 fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
641
642 for (i = 0; i < np; i++) {
643 TParticle* iparticle = (TParticle *) fParticles->At(i);
644 Int_t ks = iparticle->GetStatusCode();
645 kf = CheckPDGCode(iparticle->GetPdgCode());
646// No initial state partons
647 if (ks==21) continue;
648//
649// Heavy Flavor Selection
650//
651 // quark ?
652 kf = TMath::Abs(kf);
653 Int_t kfl = kf;
654 // Resonance
655
656 if (kfl > 100000) kfl %= 100000;
657 if (kfl > 10000) kfl %= 10000;
658 // meson ?
659 if (kfl > 10) kfl/=100;
660 // baryon
661 if (kfl > 10) kfl/=10;
662 Int_t ipa = iparticle->GetFirstMother()-1;
663 Int_t kfMo = 0;
664//
665// Establish mother daughter relation between heavy quarks and mesons
666//
667 if (kf >= fFlavorSelect && kf <= 6) {
668 Int_t idau = iparticle->GetFirstDaughter() - 1;
669 if (idau > -1) {
670 TParticle* daughter = (TParticle *) fParticles->At(idau);
671 Int_t pdgD = daughter->GetPdgCode();
672 if (pdgD == 91 || pdgD == 92) {
673 Int_t jmin = daughter->GetFirstDaughter() - 1;
674 Int_t jmax = daughter->GetLastDaughter() - 1;
675 for (Int_t j = jmin; j <= jmax; j++)
676 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
677 } // is string or cluster
678 } // has daughter
679 } // heavy quark
680
681
682 if (ipa > -1) {
683 TParticle * mother = (TParticle *) fParticles->At(ipa);
684 kfMo = TMath::Abs(mother->GetPdgCode());
685 }
686
687 // What to keep in Stack?
688 Bool_t flavorOK = kFALSE;
689 Bool_t selectOK = kFALSE;
690 if (fFeedDownOpt) {
691 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
692 } else {
693 if (kfl > fFlavorSelect) {
694 nc = -1;
695 break;
696 }
697 if (kfl == fFlavorSelect) flavorOK = kTRUE;
698 }
699 switch (fStackFillOpt) {
700 case kFlavorSelection:
701 selectOK = kTRUE;
702 break;
703 case kParentSelection:
704 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
705 break;
706 }
707 if (flavorOK && selectOK) {
708//
709// Heavy flavor hadron or quark
710//
711// Kinematic seletion on final state heavy flavor mesons
712 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
713 {
714 continue;
715 }
716 pSelected[i] = 1;
717 if (ParentSelected(kf)) ++nParents; // Update parent count
718// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
719 } else {
720// Kinematic seletion on decay products
721 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
722 && !KinematicSelection(iparticle, 1))
723 {
724 continue;
725 }
726//
727// Decay products
728// Select if mother was selected and is not tracked
729
730 if (pSelected[ipa] &&
731 !trackIt[ipa] && // mother will be tracked ?
732 kfMo != 5 && // mother is b-quark, don't store fragments
733 kfMo != 4 && // mother is c-quark, don't store fragments
734 kf != 92) // don't store string
735 {
736//
737// Semi-stable or de-selected: diselect decay products:
738//
739//
740 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
741 {
742 Int_t ipF = iparticle->GetFirstDaughter();
743 Int_t ipL = iparticle->GetLastDaughter();
744 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
745 }
746// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
747 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
748 }
749 }
750 if (pSelected[i] == -1) pSelected[i] = 0;
751 if (!pSelected[i]) continue;
752 // Count quarks only if you did not include fragmentation
753 if (fFragmentation && kf <= 10) continue;
754
755 nc++;
756// Decision on tracking
757 trackIt[i] = 0;
758//
759// Track final state particle
760 if (ks == 1) trackIt[i] = 1;
761// Track semi-stable particles
762 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
763// Track particles selected by process if undecayed.
764 if (fForceDecay == kNoDecay) {
765 if (ParentSelected(kf)) trackIt[i] = 1;
766 } else {
767 if (ParentSelected(kf)) trackIt[i] = 0;
768 }
769 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
770//
771//
772
773 } // particle selection loop
774 if (nc > 0) {
775 for (i = 0; i<np; i++) {
776 if (!pSelected[i]) continue;
777 TParticle * iparticle = (TParticle *) fParticles->At(i);
778 kf = CheckPDGCode(iparticle->GetPdgCode());
779 Int_t ks = iparticle->GetStatusCode();
780 p[0] = iparticle->Px();
781 p[1] = iparticle->Py();
782 p[2] = iparticle->Pz();
783 p[3] = iparticle->Energy();
784
785 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
786 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
787 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
788
789 Float_t tof = kconv*iparticle->T();
790 Int_t ipa = iparticle->GetFirstMother()-1;
791 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
792
793 PushTrack(fTrackIt*trackIt[i], iparent, kf,
794 p[0], p[1], p[2], p[3],
795 origin[0], origin[1], origin[2], tof,
796 polar[0], polar[1], polar[2],
797 kPPrimary, nt, 1., ks);
798 pParent[i] = nt;
799 KeepTrack(nt);
800 } // PushTrack loop
801 }
802 } else {
803 nc = GenerateMB();
804 } // mb ?
805
806 GetSubEventTime();
807
808 delete[] pParent;
809 delete[] pSelected;
810 delete[] trackIt;
811
812 if (nc > 0) {
813 switch (fCountMode) {
814 case kCountAll:
815 // printf(" Count all \n");
816 jev += nc;
817 break;
818 case kCountParents:
819 // printf(" Count parents \n");
820 jev += nParents;
821 break;
822 case kCountTrackables:
823 // printf(" Count trackable \n");
824 jev += nTkbles;
825 break;
826 }
827 if (jev >= fNpart || fNpart == -1) {
828 fKineBias=Float_t(fNpart)/Float_t(fTrials);
829
830 fQ += fPythia->GetVINT(51);
831 fX1 += fPythia->GetVINT(41);
832 fX2 += fPythia->GetVINT(42);
833 fTrialsRun += fTrials;
834 fNev++;
835 MakeHeader();
836 break;
837 }
838 }
839 } // event loop
840 SetHighWaterMark(nt);
841// adjust weight due to kinematic selection
842// AdjustWeights();
843// get cross-section
844 fXsection=fPythia->GetPARI(1);
845}
846
847Int_t AliGenPythia::GenerateMB()
848{
849//
850// Min Bias selection and other global selections
851//
852 Int_t i, kf, nt, iparent;
853 Int_t nc = 0;
854 Float_t p[4];
855 Float_t polar[3] = {0,0,0};
856 Float_t origin[3] = {0,0,0};
857// converts from mm/c to s
858 const Float_t kconv=0.001/2.999792458e8;
859
860
861
862 Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
863
864
865
866 Int_t* pParent = new Int_t[np];
867 for (i=0; i< np; i++) pParent[i] = -1;
868 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
869 TParticle* jet1 = (TParticle *) fParticles->At(6);
870 TParticle* jet2 = (TParticle *) fParticles->At(7);
871 if (!CheckTrigger(jet1, jet2)) {
872 delete [] pParent;
873 return 0;
874 }
875 }
876
877 if (fTriggerParticle) {
878 Bool_t triggered = kFALSE;
879 for (i = 0; i < np; i++) {
880 TParticle * iparticle = (TParticle *) fParticles->At(i);
881 kf = CheckPDGCode(iparticle->GetPdgCode());
882 if (kf != fTriggerParticle) continue;
883 if (iparticle->Pt() == 0.) continue;
884 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
885 triggered = kTRUE;
886 break;
887 }
888 if (!triggered) {
889 delete [] pParent;
890 return 0;
891 }
892 }
893
894
895 // Check if there is a ccbar or bbbar pair with at least one of the two
896 // in fYMin < y < fYMax
897 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
898 TParticle *hvq;
899 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
900 Float_t yQ;
901 Int_t pdgQ;
902 for(i=0; i<np; i++) {
903 hvq = (TParticle*)fParticles->At(i);
904 pdgQ = hvq->GetPdgCode();
905 if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
906 if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
907 yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
908 (hvq->Energy()-hvq->Pz()+1.e-13));
909 if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
910 }
911 if (!theQ || !theQbar || !inYcut) {
912 delete[] pParent;
913 return 0;
914 }
915 }
916
917 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
918 if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
919 && (fCutOnChild == 1) ) {
920 if ( !CheckKinematicsOnChild() ) {
921 delete[] pParent;
922 return 0;
923 }
924 }
925
926
927 for (i = 0; i < np; i++) {
928 Int_t trackIt = 0;
929 TParticle * iparticle = (TParticle *) fParticles->At(i);
930 kf = CheckPDGCode(iparticle->GetPdgCode());
931 Int_t ks = iparticle->GetStatusCode();
932 Int_t km = iparticle->GetFirstMother();
933 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
934 (ks != 1) ||
935 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
936 nc++;
937 if (ks == 1) trackIt = 1;
938 Int_t ipa = iparticle->GetFirstMother()-1;
939
940 iparent = (ipa > -1) ? pParent[ipa] : -1;
941
942//
943// store track information
944 p[0] = iparticle->Px();
945 p[1] = iparticle->Py();
946 p[2] = iparticle->Pz();
947 p[3] = iparticle->Energy();
948
949
950 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
951 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
952 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
953
954 Float_t tof = fEventTime + kconv * iparticle->T();
955
956 PushTrack(fTrackIt*trackIt, iparent, kf,
957 p[0], p[1], p[2], p[3],
958 origin[0], origin[1], origin[2], tof,
959 polar[0], polar[1], polar[2],
960 kPPrimary, nt, 1., ks);
961 //
962 // Special Treatment to store color-flow
963 //
964 if (ks == 3 || ks == 13 || ks == 14) {
965 TParticle* particle = 0;
966 if (fStack) {
967 particle = fStack->Particle(nt);
968 } else {
969 particle = gAlice->Stack()->Particle(nt);
970 }
971 particle->SetFirstDaughter(fPythia->GetK(2, i));
972 particle->SetLastDaughter(fPythia->GetK(3, i));
973 }
974
975 KeepTrack(nt);
976 pParent[i] = nt;
977 SetHighWaterMark(nt);
978
979 } // select particle
980 } // particle loop
981
982 delete[] pParent;
983
984 return 1;
985}
986
987
988void AliGenPythia::FinishRun()
989{
990// Print x-section summary
991 fPythia->Pystat(1);
992
993 if (fNev > 0.) {
994 fQ /= fNev;
995 fX1 /= fNev;
996 fX2 /= fNev;
997 }
998
999 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1000 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1001}
1002
1003void AliGenPythia::AdjustWeights() const
1004{
1005// Adjust the weights after generation of all events
1006//
1007 if (gAlice) {
1008 TParticle *part;
1009 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1010 for (Int_t i=0; i<ntrack; i++) {
1011 part= gAlice->GetMCApp()->Particle(i);
1012 part->SetWeight(part->GetWeight()*fKineBias);
1013 }
1014 }
1015}
1016
1017void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
1018{
1019// Treat protons as inside nuclei with mass numbers a1 and a2
1020
1021 fAProjectile = a1;
1022 fATarget = a2;
1023 fSetNuclei = kTRUE;
1024}
1025
1026
1027void AliGenPythia::MakeHeader()
1028{
1029//
1030// Make header for the simulated event
1031//
1032 if (gAlice) {
1033 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1034 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1035 }
1036
1037// Builds the event header, to be called after each event
1038 if (fHeader) delete fHeader;
1039 fHeader = new AliGenPythiaEventHeader("Pythia");
1040//
1041// Event type
1042 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1043//
1044// Number of trials
1045 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1046//
1047// Event Vertex
1048 fHeader->SetPrimaryVertex(fVertex);
1049//
1050// Jets that have triggered
1051
1052 if (fProcess == kPyJets)
1053 {
1054 Int_t ntrig, njet;
1055 Float_t jets[4][10];
1056 GetJets(njet, ntrig, jets);
1057
1058
1059 for (Int_t i = 0; i < ntrig; i++) {
1060 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1061 jets[3][i]);
1062 }
1063 }
1064//
1065// Copy relevant information from external header, if present.
1066//
1067 Float_t uqJet[4];
1068
1069 if (fRL) {
1070 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1071 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1072 {
1073 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
1074
1075
1076 exHeader->TriggerJet(i, uqJet);
1077 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1078 }
1079 }
1080//
1081// Store quenching parameters
1082//
1083 if (fQuench){
1084 Double_t z[4];
1085 Double_t xp, yp;
1086 if (fQuench == 1) {
1087 // Pythia::Quench()
1088 fPythia->GetQuenchingParameters(xp, yp, z);
1089 } else {
1090 // Pyquen
1091 Double_t r1 = PARIMP.rb1;
1092 Double_t r2 = PARIMP.rb2;
1093 Double_t b = PARIMP.b1;
1094 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1095 Double_t phi = PARIMP.psib1;
1096 xp = r * TMath::Cos(phi);
1097 yp = r * TMath::Sin(phi);
1098
1099 }
1100 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1101 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1102 }
1103//
1104// Store pt^hard
1105 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1106//
1107// Pass header
1108//
1109 AddHeader(fHeader);
1110}
1111
1112void AliGenPythia::AddHeader(AliGenEventHeader* header)
1113{
1114 // Add header to container or runloader
1115 if (fContainer) {
1116 fContainer->AddHeader(header);
1117 } else {
1118 AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
1119 }
1120}
1121
1122
1123Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1124{
1125// Check the kinematic trigger condition
1126//
1127 Double_t eta[2];
1128 eta[0] = jet1->Eta();
1129 eta[1] = jet2->Eta();
1130 Double_t phi[2];
1131 phi[0] = jet1->Phi();
1132 phi[1] = jet2->Phi();
1133 Int_t pdg[2];
1134 pdg[0] = jet1->GetPdgCode();
1135 pdg[1] = jet2->GetPdgCode();
1136 Bool_t triggered = kFALSE;
1137
1138 if (fProcess == kPyJets) {
1139 Int_t njets = 0;
1140 Int_t ntrig = 0;
1141 Float_t jets[4][10];
1142//
1143// Use Pythia clustering on parton level to determine jet axis
1144//
1145 GetJets(njets, ntrig, jets);
1146
1147 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1148//
1149 } else {
1150 Int_t ij = 0;
1151 Int_t ig = 1;
1152 if (pdg[0] == kGamma) {
1153 ij = 1;
1154 ig = 0;
1155 }
1156 //Check eta range first...
1157 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1158 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1159 {
1160 //Eta is okay, now check phi range
1161 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1162 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1163 {
1164 triggered = kTRUE;
1165 }
1166 }
1167 }
1168 return triggered;
1169}
1170
1171
1172
1173Bool_t AliGenPythia::CheckKinematicsOnChild(){
1174//
1175//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1176//
1177 Bool_t checking = kFALSE;
1178 Int_t j, kcode, ks, km;
1179 Int_t nPartAcc = 0; //number of particles in the acceptance range
1180 Int_t numberOfAcceptedParticles = 1;
1181 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1182 Int_t npart = fParticles->GetEntriesFast();
1183
1184 for (j = 0; j<npart; j++) {
1185 TParticle * jparticle = (TParticle *) fParticles->At(j);
1186 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1187 ks = jparticle->GetStatusCode();
1188 km = jparticle->GetFirstMother();
1189
1190 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1191 nPartAcc++;
1192 }
1193 if( numberOfAcceptedParticles <= nPartAcc){
1194 checking = kTRUE;
1195 break;
1196 }
1197 }
1198
1199 return checking;
1200}
1201
1202
1203AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
1204{
1205// Assignment operator
1206 rhs.Copy(*this);
1207 return *this;
1208}
1209
1210void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1211{
1212//
1213// Load event into Pythia Common Block
1214//
1215
1216 Int_t npart = stack -> GetNprimary();
1217 Int_t n0 = 0;
1218
1219 if (!flag) {
1220 (fPythia->GetPyjets())->N = npart;
1221 } else {
1222 n0 = (fPythia->GetPyjets())->N;
1223 (fPythia->GetPyjets())->N = n0 + npart;
1224 }
1225
1226
1227 for (Int_t part = 0; part < npart; part++) {
1228 TParticle *mPart = stack->Particle(part);
1229
1230 Int_t kf = mPart->GetPdgCode();
1231 Int_t ks = mPart->GetStatusCode();
1232 Int_t idf = mPart->GetFirstDaughter();
1233 Int_t idl = mPart->GetLastDaughter();
1234
1235 if (reHadr) {
1236 if (ks == 11 || ks == 12) {
1237 ks -= 10;
1238 idf = -1;
1239 idl = -1;
1240 }
1241 }
1242
1243 Float_t px = mPart->Px();
1244 Float_t py = mPart->Py();
1245 Float_t pz = mPart->Pz();
1246 Float_t e = mPart->Energy();
1247 Float_t m = mPart->GetCalcMass();
1248
1249
1250 (fPythia->GetPyjets())->P[0][part+n0] = px;
1251 (fPythia->GetPyjets())->P[1][part+n0] = py;
1252 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1253 (fPythia->GetPyjets())->P[3][part+n0] = e;
1254 (fPythia->GetPyjets())->P[4][part+n0] = m;
1255
1256 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1257 (fPythia->GetPyjets())->K[0][part+n0] = ks;
1258 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1259 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1260 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1261 }
1262}
1263
1264
1265void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1266{
1267//
1268// Calls the Pythia jet finding algorithm to find jets in the current event
1269//
1270//
1271//
1272// Save jets
1273 Int_t n = fPythia->GetN();
1274
1275//
1276// Run Jet Finder
1277 fPythia->Pycell(njets);
1278 Int_t i;
1279 for (i = 0; i < njets; i++) {
1280 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1281 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1282 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1283 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1284
1285 jets[0][i] = px;
1286 jets[1][i] = py;
1287 jets[2][i] = pz;
1288 jets[3][i] = e;
1289 }
1290}
1291
1292
1293
1294void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1295{
1296//
1297// Calls the Pythia clustering algorithm to find jets in the current event
1298//
1299 Int_t n = fPythia->GetN();
1300 nJets = 0;
1301 nJetsTrig = 0;
1302 if (fJetReconstruction == kCluster) {
1303//
1304// Configure cluster algorithm
1305//
1306 fPythia->SetPARU(43, 2.);
1307 fPythia->SetMSTU(41, 1);
1308//
1309// Call cluster algorithm
1310//
1311 fPythia->Pyclus(nJets);
1312//
1313// Loading jets from common block
1314//
1315 } else {
1316
1317//
1318// Run Jet Finder
1319 fPythia->Pycell(nJets);
1320 }
1321
1322 Int_t i;
1323 for (i = 0; i < nJets; i++) {
1324 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1325 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1326 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1327 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1328 Float_t pt = TMath::Sqrt(px * px + py * py);
1329 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
1330 Float_t theta = TMath::ATan2(pt,pz);
1331 Float_t et = e * TMath::Sin(theta);
1332 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
1333 if (
1334 eta > fEtaMinJet && eta < fEtaMaxJet &&
1335 phi > fPhiMinJet && phi < fPhiMaxJet &&
1336 et > fEtMinJet && et < fEtMaxJet
1337 )
1338 {
1339 jets[0][nJetsTrig] = px;
1340 jets[1][nJetsTrig] = py;
1341 jets[2][nJetsTrig] = pz;
1342 jets[3][nJetsTrig] = e;
1343 nJetsTrig++;
1344// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1345 } else {
1346// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1347 }
1348 }
1349}
1350
1351void AliGenPythia::GetSubEventTime()
1352{
1353 // Calculates time of the next subevent
1354 fEventTime = 0.;
1355 if (fEventsTime) {
1356 TArrayF &array = *fEventsTime;
1357 fEventTime = array[fCurSubEvent++];
1358 }
1359 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1360 return;
1361}
1362
1363#ifdef never
1364void AliGenPythia::Streamer(TBuffer &R__b)
1365{
1366 // Stream an object of class AliGenPythia.
1367
1368 if (R__b.IsReading()) {
1369 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1370 AliGenerator::Streamer(R__b);
1371 R__b >> (Int_t&)fProcess;
1372 R__b >> (Int_t&)fStrucFunc;
1373 R__b >> (Int_t&)fForceDecay;
1374 R__b >> fEnergyCMS;
1375 R__b >> fKineBias;
1376 R__b >> fTrials;
1377 fParentSelect.Streamer(R__b);
1378 fChildSelect.Streamer(R__b);
1379 R__b >> fXsection;
1380// (AliPythia::Instance())->Streamer(R__b);
1381 R__b >> fPtHardMin;
1382 R__b >> fPtHardMax;
1383// if (fDecayer) fDecayer->Streamer(R__b);
1384 } else {
1385 R__b.WriteVersion(AliGenPythia::IsA());
1386 AliGenerator::Streamer(R__b);
1387 R__b << (Int_t)fProcess;
1388 R__b << (Int_t)fStrucFunc;
1389 R__b << (Int_t)fForceDecay;
1390 R__b << fEnergyCMS;
1391 R__b << fKineBias;
1392 R__b << fTrials;
1393 fParentSelect.Streamer(R__b);
1394 fChildSelect.Streamer(R__b);
1395 R__b << fXsection;
1396// R__b << fPythia;
1397 R__b << fPtHardMin;
1398 R__b << fPtHardMax;
1399 // fDecayer->Streamer(R__b);
1400 }
1401}
1402#endif
1403
1404
1405