Add kinematic bias for direct gamma production.
[u/mrichter/AliRoot.git] / EVGEN / AliGenPythia.cxx
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$
18 Revision 1.47  2001/12/19 14:45:00  morsch
19 Store number of trials in header.
20
21 Revision 1.46  2001/12/19 10:36:19  morsch
22 Add possibility if jet kinematic biasing.
23
24 Revision 1.45  2001/11/28 08:06:52  morsch
25 Use fMaxLifeTime parameter.
26
27 Revision 1.44  2001/11/27 13:13:07  morsch
28 Maximum lifetime for long-lived particles to be put on the stack is parameter.
29 It can be set via SetMaximumLifetime(..).
30
31 Revision 1.43  2001/10/21 18:35:56  hristov
32 Several pointers were set to zero in the default constructors to avoid memory management problems
33
34 Revision 1.42  2001/10/15 08:21:55  morsch
35 Vertex truncation settings moved to AliGenMC.
36
37 Revision 1.41  2001/10/08 08:45:42  morsch
38 Possibility of vertex cut added.
39
40 Revision 1.40  2001/09/25 11:30:23  morsch
41 Pass event vertex to header.
42
43 Revision 1.39  2001/07/27 17:09:36  morsch
44 Use local SetTrack, KeepTrack and SetHighWaterMark methods
45 to delegate either to local stack or to stack owned by AliRun.
46 (Piotr Skowronski, A.M.)
47
48 Revision 1.38  2001/07/13 10:58:54  morsch
49 - Some coded moved to AliGenMC
50 - Improved handling of secondary vertices.
51
52 Revision 1.37  2001/06/28 11:17:28  morsch
53 SetEventListRange setter added. Events in specified range are listed for
54 debugging. (Yuri Kharlov)
55
56 Revision 1.36  2001/03/30 07:05:49  morsch
57 Final print-out in finish run.
58 Write parton system for jet-production (preliminary solution).
59
60 Revision 1.35  2001/03/09 13:03:40  morsch
61 Process_t and Struc_Func_t moved to AliPythia.h
62
63 Revision 1.34  2001/02/14 15:50:40  hristov
64 The last particle in event marked using SetHighWaterMark
65
66 Revision 1.33  2001/01/30 09:23:12  hristov
67 Streamers removed (R.Brun)
68
69 Revision 1.32  2001/01/26 19:55:51  hristov
70 Major upgrade of AliRoot code
71
72 Revision 1.31  2001/01/17 10:54:31  hristov
73 Better protection against FPE
74
75 Revision 1.30  2000/12/18 08:55:35  morsch
76 Make AliPythia dependent generartors work with new scheme of random number generation
77
78 Revision 1.29  2000/12/04 11:22:03  morsch
79 Init of sRandom as in 1.15
80
81 Revision 1.28  2000/12/02 11:41:39  morsch
82 Use SetRandom() to initialize random number generator in constructor.
83
84 Revision 1.27  2000/11/30 20:29:02  morsch
85 Initialise static variable sRandom in constructor: sRandom = fRandom;
86
87 Revision 1.26  2000/11/30 07:12:50  alibrary
88 Introducing new Rndm and QA classes
89
90 Revision 1.25  2000/10/18 19:11:27  hristov
91 Division by zero fixed
92
93 Revision 1.24  2000/09/18 10:41:35  morsch
94 Add possibility to use nuclear structure functions from PDF library V8.
95
96 Revision 1.23  2000/09/14 14:05:40  morsch
97 dito
98
99 Revision 1.22  2000/09/14 14:02:22  morsch
100 - Correct conversion from mm to cm when passing particle vertex to MC.
101 - Correct handling of fForceDecay == all.
102
103 Revision 1.21  2000/09/12 14:14:55  morsch
104 Call fDecayer->ForceDecay() at the beginning of Generate().
105
106 Revision 1.20  2000/09/06 14:29:33  morsch
107 Use AliPythia for event generation an AliDecayPythia for decays.
108 Correct handling of "nodecay" option
109
110 Revision 1.19  2000/07/11 18:24:56  fca
111 Coding convention corrections + few minor bug fixes
112
113 Revision 1.18  2000/06/30 12:40:34  morsch
114 Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
115 Geant units (cm).
116
117 Revision 1.17  2000/06/09 20:34:07  morsch
118 All coding rule violations except RS3 corrected
119
120 Revision 1.16  2000/05/15 15:04:20  morsch
121 The full event is written for fNtrack = -1
122 Coding rule violations corrected.
123
124 Revision 1.15  2000/04/26 10:14:24  morsch
125 Particles array has one entry more than pythia particle list. Upper bound of
126 particle loop changed to np-1 (R. Guernane, AM)
127
128 Revision 1.14  2000/04/05 08:36:13  morsch
129 Check status code of particles in Pythia event
130 to avoid double counting as partonic state and final state particle.
131
132 Revision 1.13  1999/11/09 07:38:48  fca
133 Changes for compatibility with version 2.23 of ROOT
134
135 Revision 1.12  1999/11/03 17:43:20  fca
136 New version from G.Martinez & A.Morsch
137
138 Revision 1.11  1999/09/29 09:24:14  fca
139 Introduction of the Copyright and cvs Log
140 */
141
142 #include "AliGenPythia.h"
143 #include "AliGenPythiaEventHeader.h"
144 #include "AliDecayerPythia.h"
145 #include "AliRun.h"
146 #include "AliPythia.h"
147 #include "AliPDG.h"
148 #include <TParticle.h>
149 #include <TSystem.h>
150
151  ClassImp(AliGenPythia)
152
153 AliGenPythia::AliGenPythia()
154                  :AliGenMC()
155 {
156 // Default Constructor
157   fParticles = 0;
158   fPythia    = 0;
159   fDecayer   = new AliDecayerPythia();
160   SetEventListRange();
161   SetJetPhiRange();
162   SetJetEtaRange();
163 }
164
165 AliGenPythia::AliGenPythia(Int_t npart)
166                  :AliGenMC(npart)
167 {
168 // default charm production at 5. 5 TeV
169 // semimuonic decay
170 // structure function GRVHO
171 //
172     fXsection  = 0.;
173     fNucA1=0;
174     fNucA2=0;
175     SetProcess();
176     SetStrucFunc();
177     SetForceDecay();
178     SetPtHard();
179     SetEnergyCMS();
180     fDecayer = new AliDecayerPythia();
181     // Set random number generator 
182     sRandom=fRandom;
183     fFlavorSelect   = 0;
184     // Produced particles  
185     fParticles = new TClonesArray("TParticle",1000);
186     fEventVertex.Set(3);
187     SetEventListRange();
188     SetJetPhiRange();
189     SetJetEtaRange();
190 }
191
192 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
193 {
194 // copy constructor
195 }
196
197 AliGenPythia::~AliGenPythia()
198 {
199 // Destructor
200 }
201
202 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
203 {
204   // Set a range of event numbers, for which a table
205   // of generated particle will be printed
206   fDebugEventFirst = eventFirst;
207   fDebugEventLast  = eventLast;
208   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
209 }
210
211 void AliGenPythia::Init()
212 {
213 // Initialisation
214   SetMC(AliPythia::Instance());
215     fPythia=(AliPythia*) fgMCEvGen;
216 //
217     fParentWeight=1./Float_t(fNpart);
218 //
219 //  Forward Paramters to the AliPythia object
220     fDecayer->SetForceDecay(fForceDecay);    
221     fDecayer->Init();
222
223
224     fPythia->SetCKIN(3,fPtHardMin);
225     fPythia->SetCKIN(4,fPtHardMax);    
226     if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);  
227     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
228
229     //    fPythia->Pylist(0);
230     //    fPythia->Pystat(2);
231 //  Parent and Children Selection
232     switch (fProcess) 
233     {
234     case kPyCharm:
235         fParentSelect[0] =  411;
236         fParentSelect[1] =  421;
237         fParentSelect[2] =  431;
238         fParentSelect[3] = 4122;        
239         fFlavorSelect    = 4;
240         break;
241     case kPyCharmUnforced:
242         fParentSelect[0] =   411;
243         fParentSelect[1] =   421;
244         fParentSelect[2] =   431;
245         fParentSelect[3] =  4122;
246         fFlavorSelect    =     4;       
247         break;
248     case kPyBeauty:
249         fParentSelect[0]=  511;
250         fParentSelect[1]=  521;
251         fParentSelect[2]=  531;
252         fParentSelect[3]= 5122;
253         fParentSelect[4]= 5132;
254         fParentSelect[5]= 5232;
255         fParentSelect[6]= 5332;
256         fFlavorSelect   = 5;    
257         break;
258     case kPyBeautyUnforced:
259         fParentSelect[0] =  511;
260         fParentSelect[1] =  521;
261         fParentSelect[2] =  531;
262         fParentSelect[3] = 5122;
263         fParentSelect[4] = 5132;
264         fParentSelect[5] = 5232;
265         fParentSelect[6] = 5332;
266         fFlavorSelect    = 5;   
267         break;
268     case kPyJpsiChi:
269     case kPyJpsi:
270         fParentSelect[0] = 443;
271         break;
272     case kPyMb:
273     case kPyJets:
274     case kPyDirectGamma:
275         break;
276     }
277     AliGenMC::Init();
278 }
279
280 void AliGenPythia::Generate()
281 {
282 // Generate one event
283     fDecayer->ForceDecay();
284
285     Float_t polar[3]   =   {0,0,0};
286     Float_t origin[3]  =   {0,0,0};
287     Float_t p[3];
288 //  converts from mm/c to s
289     const Float_t kconv=0.001/2.999792458e8;
290 //
291     Int_t nt=0;
292     Int_t jev=0;
293     Int_t j, kf;
294     fTrials=0;
295
296     //  Set collision vertex position 
297     if(fVertexSmear==kPerEvent) {
298         fPythia->SetMSTP(151,1);
299         for (j=0;j<3;j++) {
300             fPythia->SetPARP(151+j, fOsigma[j]*10.);
301         }
302     } else if (fVertexSmear==kPerTrack) {
303         fPythia->SetMSTP(151,0);
304     }
305 //  event loop    
306     while(1)
307     {
308         fPythia->Pyevnt();
309         if (gAlice->GetEvNumber()>=fDebugEventFirst &&
310             gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
311         fTrials++;
312         
313         fPythia->ImportParticles(fParticles,"All");
314
315 //
316 //
317 //
318         Int_t i;
319         
320         Int_t np = fParticles->GetEntriesFast();
321         if (np == 0 ) continue;
322 // Get event vertex and discard the event if the Z coord. is too big    
323         TParticle *iparticle = (TParticle *) fParticles->At(0);
324         Float_t distz = iparticle->Vz()/10.;
325         if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
326 //
327         fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
328         fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
329         fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
330 //
331         Int_t* pParent   = new Int_t[np];
332         Int_t* pSelected = new Int_t[np];
333         Int_t* trackIt   = new Int_t[np];
334         for (i=0; i< np-1; i++) {
335             pParent[i]   = -1;
336             pSelected[i] =  0;
337         }
338         printf("\n **************************************************%d\n",np);
339         Int_t nc = 0;
340         if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) {
341             
342             for (i = 0; i<np-1; i++) {
343                 iparticle = (TParticle *) fParticles->At(i);
344                 Int_t ks = iparticle->GetStatusCode();
345                 kf = CheckPDGCode(iparticle->GetPdgCode());
346 // No initial state partons
347                 if (ks==21) continue;
348 //
349 // Heavy Flavor Selection
350 //
351                 // quark ?
352                 kf = TMath::Abs(kf);
353                 Int_t kfl = kf;
354                 // meson ?
355                 if  (kfl > 10) kfl/=100;
356                 // baryon
357                 if (kfl > 10) kfl/=10;
358                 if (kfl > 10) kfl/=10;
359
360                 Int_t ipa = iparticle->GetFirstMother()-1;
361                 Int_t kfMo = 0;
362                 
363                 if (ipa > -1) {
364                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
365                     kfMo = TMath::Abs(mother->GetPdgCode());
366                 }
367 //              printf("\n particle (all)  %d %d %d", i, pSelected[i], kf);
368                 if (kfl >= fFlavorSelect) { 
369 //
370 // Heavy flavor hadron or quark
371 //
372 // Kinematic seletion on final state heavy flavor mesons
373                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) 
374                     {
375                         nc = -1;
376                         break;
377                     }
378                     pSelected[i] = 1;
379 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
380                 } else {
381 // Kinematic seletion on decay products
382                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) 
383                         && !KinematicSelection(iparticle, 1))
384                     {
385                         nc = -1;
386                         break;
387                     }
388 //
389 // Decay products 
390 // Select if mother was selected and is not tracked
391
392                     if (pSelected[ipa] && 
393                         !trackIt[ipa]  &&     // mother will be  tracked ?
394                         kfMo !=  5 &&         // mother is b-quark, don't store fragments          
395                         kfMo !=  4 &&         // mother is c-quark, don't store fragments 
396                         kf   != 92)           // don't store string
397                     {
398 //
399 // Semi-stable or de-selected: diselect decay products:
400 // 
401 //
402                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
403                         {
404                             Int_t ipF = iparticle->GetFirstDaughter();
405                             Int_t ipL = iparticle->GetLastDaughter();   
406                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
407                         }
408 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
409                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
410                     }
411                 }
412                 if (pSelected[i] == -1) pSelected[i] = 0;
413                 if (!pSelected[i]) continue;
414                 nc++;
415 // Decision on tracking
416                 trackIt[i] = 0;
417 //
418 // Track final state particle
419                 if (ks == 1) trackIt[i] = 1;
420 // Track semi-stable particles
421                 if ((ks ==1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
422 // Track particles selected by process if undecayed. 
423                 if (fForceDecay == kNoDecay) {
424                     if (ParentSelected(kf)) trackIt[i] = 1;
425                 } else {
426                     if (ParentSelected(kf)) trackIt[i] = 0;
427                 }
428 //
429 //
430
431             } // particle selection loop
432             if (nc > -1) {
433                 for (i = 0; i<np-1; i++) {
434                     if (!pSelected[i]) continue;
435                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
436                     kf = CheckPDGCode(iparticle->GetPdgCode());
437                     p[0] = iparticle->Px();
438                     p[1] = iparticle->Py();
439                     p[2] = iparticle->Pz();
440                     origin[0] = fOrigin[0]+iparticle->Vx()/10.;
441                     origin[1] = fOrigin[1]+iparticle->Vy()/10.;
442                     origin[2] = fOrigin[2]+iparticle->Vz()/10.;
443                     Float_t tof   = kconv*iparticle->T();
444                     Int_t ipa     = iparticle->GetFirstMother()-1;
445                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
446                     SetTrack(fTrackIt*trackIt[i] ,
447                                      iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.);
448                     pParent[i] = nt;
449                     KeepTrack(nt); 
450                 } //  SetTrack loop
451             }
452         } else {
453             nc = GenerateMB();
454         } // mb ?
455
456         if (nc > 0) {
457             jev+=nc;
458             if (jev >= fNpart || fNpart == -1) {
459                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
460                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
461                 MakeHeader();
462                 break;
463             }
464         }
465     } // event loop
466     SetHighWaterMark(nt);
467 //  adjust weight due to kinematic selection
468     AdjustWeights();
469 //  get cross-section
470     fXsection=fPythia->GetPARI(1);
471 }
472
473 Int_t  AliGenPythia::GenerateMB()
474 {
475     Int_t i, kf, nt, iparent;
476     Int_t nc = 0;
477     Float_t p[3];
478     Float_t polar[3]   =   {0,0,0};
479     Float_t origin[3]  =   {0,0,0};
480 //  converts from mm/c to s
481     const Float_t kconv=0.001/2.999792458e8;
482     
483     Int_t np = fParticles->GetEntriesFast();
484     Int_t* pParent = new Int_t[np];
485     for (i=0; i< np-1; i++) pParent[i] = -1;
486     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
487         TParticle* jet1 = (TParticle *) fParticles->At(6);
488         TParticle* jet2 = (TParticle *) fParticles->At(7);
489         if (!CheckTrigger(jet1, jet2)) return 0;
490     }
491     
492     for (i = 0; i<np-1; i++) {
493         Int_t trackIt = 0;
494         TParticle *  iparticle = (TParticle *) fParticles->At(i);
495         kf = CheckPDGCode(iparticle->GetPdgCode());
496         Int_t ks = iparticle->GetStatusCode();
497         Int_t km = iparticle->GetFirstMother();
498 //      printf("\n Particle: %d %d %d", i, kf, ks);
499         
500         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
501             (ks != 1) ||
502             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
503             nc++;
504             if (ks == 1) trackIt = 1;
505             Int_t ipa = iparticle->GetFirstMother()-1;
506
507             iparent = (ipa > -1) ? pParent[ipa] : -1;
508
509 //
510 // store track information
511             p[0] = iparticle->Px();
512             p[1] = iparticle->Py();
513             p[2] = iparticle->Pz();
514             origin[0] = fOrigin[0]+iparticle->Vx()/10.;
515             origin[1] = fOrigin[1]+iparticle->Vy()/10.;
516             origin[2] = fOrigin[2]+iparticle->Vz()/10.;
517             Float_t tof=kconv*iparticle->T();
518             SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
519                          tof, kPPrimary, nt);
520             KeepTrack(nt);
521             pParent[i] = nt;
522         } // select particle
523     } // particle loop 
524
525     if (pParent) delete[] pParent;
526     
527     printf("\n I've put %i particles on the stack \n",nc);
528     return nc;
529 }
530
531
532 void AliGenPythia::FinishRun()
533 {
534 // Print x-section summary
535     fPythia->Pystat(1);
536 }
537
538 void AliGenPythia::AdjustWeights()
539 {
540 // Adjust the weights after generation of all events
541 //
542     TParticle *part;
543     Int_t ntrack=gAlice->GetNtrack();
544     for (Int_t i=0; i<ntrack; i++) {
545         part= gAlice->Particle(i);
546         part->SetWeight(part->GetWeight()*fKineBias);
547     }
548 }
549     
550 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
551 {
552 // Treat protons as inside nuclei with mass numbers a1 and a2  
553     fNucA1 = a1;
554     fNucA2 = a2;
555 }
556
557
558 void AliGenPythia::MakeHeader()
559 {
560 // Builds the event header, to be called after each event
561     AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
562     ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
563     ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials);
564     header->SetPrimaryVertex(fEventVertex);
565     gAlice->SetGenEventHeader(header);
566 }
567         
568
569 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
570 {
571 // Check the kinematic trigger condition
572 //
573     Double_t eta[2];
574     eta[0] = jet1->Eta();
575     eta[1] = jet2->Eta();
576     Double_t phi[2];
577     phi[0] = jet1->Phi();
578     phi[1] = jet2->Phi();
579     Int_t    pdg[2]; 
580     pdg[0] = jet1->GetPdgCode();
581     pdg[1] = jet2->GetPdgCode();    
582     Bool_t   triggered = kFALSE;
583
584     if (fProcess == kPyJets) {
585         //Check eta range first...
586         if ((eta[0] < fEtaMaxJet && eta[0] > fEtaMinJet) ||
587             (eta[1] < fEtaMaxJet && eta[1] > fEtaMinJet))
588         {
589             //Eta is okay, now check phi range
590             if ((phi[0] < fPhiMaxJet && phi[0] > fPhiMinJet) ||
591                 (phi[1] < fPhiMaxJet && phi[1] > fPhiMinJet))
592             {
593                 triggered = kTRUE;
594             }
595         }
596     } else {
597         Int_t ij = 0;
598         Int_t ig = 1;
599         if (pdg[0] == kGamma) {
600             ij = 1;
601             ig = 0;
602         }
603         //Check eta range first...
604         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
605             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
606         {
607             //Eta is okay, now check phi range
608             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
609                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
610             {
611                 triggered = kTRUE;
612             }
613         }
614     }
615     
616     
617     return triggered;
618 }
619           
620 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
621 {
622 // Assignment operator
623     return *this;
624 }
625
626
627
628 #ifdef never
629 void AliGenPythia::Streamer(TBuffer &R__b)
630 {
631    // Stream an object of class AliGenPythia.
632
633    if (R__b.IsReading()) {
634       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
635       AliGenerator::Streamer(R__b);
636       R__b >> (Int_t&)fProcess;
637       R__b >> (Int_t&)fStrucFunc;
638       R__b >> (Int_t&)fForceDecay;
639       R__b >> fEnergyCMS;
640       R__b >> fKineBias;
641       R__b >> fTrials;
642       fParentSelect.Streamer(R__b);
643       fChildSelect.Streamer(R__b);
644       R__b >> fXsection;
645 //      (AliPythia::Instance())->Streamer(R__b);
646       R__b >> fPtHardMin;
647       R__b >> fPtHardMax;
648 //      if (fDecayer) fDecayer->Streamer(R__b);
649    } else {
650       R__b.WriteVersion(AliGenPythia::IsA());
651       AliGenerator::Streamer(R__b);
652       R__b << (Int_t)fProcess;
653       R__b << (Int_t)fStrucFunc;
654       R__b << (Int_t)fForceDecay;
655       R__b << fEnergyCMS;
656       R__b << fKineBias;
657       R__b << fTrials;
658       fParentSelect.Streamer(R__b);
659       fChildSelect.Streamer(R__b);
660       R__b << fXsection;
661 //      R__b << fPythia;
662       R__b << fPtHardMin;
663       R__b << fPtHardMax;
664       //     fDecayer->Streamer(R__b);
665    }
666 }
667 #endif
668