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