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