Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGenerator.cxx
Go to the documentation of this file.
1 // @(#)root/eg:$Id$
2 // Author: Ola Nordmann 21/09/95
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class TGenerator
13  \ingroup eg
14 
15 The interface to various event generators
16 
17 Is an base class, that defines the interface of ROOT to various
18 event generators. Every event generator should inherit from
19 TGenerator or its subclasses.
20 
21 Derived class can overload the member function GenerateEvent
22 to do the actual event generation (e.g., call PYEVNT or similar).
23 
24 The derived class should overload the member function
25 ImportParticles (both types) to read the internal storage of the
26 generated event into either the internal TObjArray or the passed
27 TClonesArray of TParticles.
28 
29 If the generator code stores event data in the /HEPEVT/ common block
30 Then the default implementation of ImportParticles should suffice.
31 The common block /HEPEVT/ is structed like
32 
33 \verbatim
34  // C
35  typedef struct {
36  Int_t nevhep; // Event number
37  Int_t nhep; // # of particles
38  Int_t isthep[4000]; // Status flag of i'th particle
39  Int_t idhep[4000]; // PDG # of particle
40  Int_t jmohep[4000][2]; // 1st & 2nd mother particle #
41  Int_t jdahep[4000][2]; // 1st & 2nd daughter particle #
42  Double_t phep[4000][5]; // 4-momentum and 1 word
43  Double_t vhep[4000][4]; // 4-position of production
44  } HEPEVT_DEF;
45 
46 
47  C Fortran
48  COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(4000),IDHEP(4000),
49  + JMOHEP(2,4000),JDAHEP(2,4000),PHEP(5,4000),VHEP(4,4000)
50  INTEGER NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
51  DOUBLE PRECISION PHEP,VHEP
52 \endverbatim
53 
54 The generic member functions SetParameter and GetParameter can be
55 overloaded to set and get parameters of the event generator.
56 
57 Note, if the derived class interfaces a (set of) Fortran common
58 blocks (like TPythia, TVenus does), one better make the derived
59 class a singleton. That is, something like
60 
61 \verbatim
62  class MyGenerator : public TGenerator
63  {
64  public:
65  static MyGenerator* Instance()
66  {
67  if (!fgInstance) fgInstance = new MyGenerator;
68  return fgInstance;
69  }
70  void GenerateEvent() { ... }
71  void ImportParticles(TClonesArray* a, Option_t opt="") {...}
72  Int_t ImportParticles(Option_t opt="") { ... }
73  Int_t SetParameter(const char* name, Double_t val) { ... }
74  Double_t GetParameter(const char* name) { ... }
75  virtual ~MyGenerator() { ... }
76  protected:
77  MyGenerator() { ... }
78  MyGenerator(const MyGenerator& o) { ... }
79  MyGenerator& operator=(const MyGenerator& o) { ... }
80  static MyGenerator* fgInstance;
81  ClassDef(MyGenerator,0);
82  };
83 \endverbatim
84 
85 Having multiple objects accessing the same common blocks is not
86 safe.
87 
88 Concrete TGenerator classes can be loaded in scripts and subseqent-
89 ly used in compiled code:
90 
91 \verbatim
92  // MyRun.h
93  class MyRun : public TObject
94  {
95  public:
96  static MyRun* Instance() { ... }
97  void SetGenerator(TGenerator* g) { fGenerator = g; }
98  void Run(Int_t n, Option_t* option="")
99  {
100  TFile* file = TFile::Open("file.root","RECREATE");
101  TTree* tree = new TTree("T","T");
102  TClonesArray* p = new TClonesArray("TParticles");
103  tree->Branch("particles", &p);
104  for (Int_t event = 0; event < n; event++) {
105  fGenerator->GenerateEvent();
106  fGenerator->ImportParticles(p,option);
107  tree->Fill();
108  }
109  file->Write();
110  file->Close();
111  }
112  ...
113  protected:
114  TGenerator* fGenerator;
115  ClassDef(MyRun,0);
116  };
117 
118  // Config.C
119  void Config()
120  {
121  MyRun* run = MyRun::Instance();
122  run->SetGenerator(MyGenerator::Instance());
123  }
124 
125  // main.cxx
126  int
127  main(int argc, char** argv)
128  {
129  TApplication app("", 0, 0);
130  gSystem->ProcessLine(".x Config.C");
131  MyRun::Instance()->Run(10);
132  return 0;
133  }
134 \endverbatim
135 
136 This is especially useful for example with TVirtualMC or similar.
137 */
138 
139 #include "TROOT.h"
140 #include "TGenerator.h"
141 #include "TDatabasePDG.h"
142 #include "TParticlePDG.h"
143 #include "TParticle.h"
144 #include "TObjArray.h"
145 #include "Hepevt.h"
146 #include "TVirtualPad.h"
147 #include "TView.h"
148 #include "TText.h"
149 #include "TPaveText.h"
150 #include "TClonesArray.h"
151 #include "Riostream.h"
152 
153 
154 ClassImp(TGenerator);
155 
156 ////////////////////////////////////////////////////////////////////////////////
157 /// Event generator default constructor
158 ///
159 
160 TGenerator::TGenerator(const char *name,const char *title): TNamed(name,title)
161 {
162  // Initialize particles table
163  TDatabasePDG::Instance();
164  //TDatabasePDG *pdg = TDatabasePDG::Instance();
165  //if (!pdg->ParticleList()) pdg->Init();
166 
167  fPtCut = 0;
168  fShowNeutrons = kTRUE;
169  fParticles = new TObjArray(10000);
170 }
171 
172 ////////////////////////////////////////////////////////////////////////////////
173 /// Event generator default destructor
174 ///
175 
176 TGenerator::~TGenerator()
177 {
178  //do nothing
179  if (fParticles) {
180  fParticles->Delete();
181  delete fParticles;
182  fParticles = 0;
183  }
184 }
185 
186 ////////////////////////////////////////////////////////////////////////////////
187 /// must be implemented in concrete class (see eg TPythia6)
188 
189 void TGenerator::GenerateEvent()
190 {
191 }
192 
193 ////////////////////////////////////////////////////////////////////////////////
194 ///
195 /// It reads the /HEPEVT/ common block which has been filled by the
196 /// GenerateEvent method. If the event generator does not use the
197 /// HEPEVT common block, This routine has to be overloaded by the
198 /// subclasses.
199 ///
200 /// The default action is to store only the stable particles (ISTHEP =
201 /// 1) This can be demanded explicitly by setting the option = "Final"
202 /// If the option = "All", all the particles are stored.
203 ///
204 
205 TObjArray* TGenerator::ImportParticles(Option_t *option)
206 {
207  fParticles->Clear();
208  Int_t numpart = HEPEVT.nhep;
209  if (!strcmp(option,"") || !strcmp(option,"Final")) {
210  for (Int_t i = 0; i<numpart; i++) {
211  if (HEPEVT.isthep[i] == 1) {
212 //
213 // Use the common block values for the TParticle constructor
214 //
215  TParticle *p = new TParticle(
216  HEPEVT.idhep[i],
217  HEPEVT.isthep[i],
218  HEPEVT.jmohep[i][0]-1,
219  HEPEVT.jmohep[i][1]-1,
220  HEPEVT.jdahep[i][0]-1,
221  HEPEVT.jdahep[i][1]-1,
222  HEPEVT.phep[i][0],
223  HEPEVT.phep[i][1],
224  HEPEVT.phep[i][2],
225  HEPEVT.phep[i][3],
226  HEPEVT.vhep[i][0],
227  HEPEVT.vhep[i][1],
228  HEPEVT.vhep[i][2],
229  HEPEVT.vhep[i][3]);
230  fParticles->Add(p);
231  }
232  }
233  } else if (!strcmp(option,"All")) {
234  for (Int_t i = 0; i<numpart; i++) {
235  TParticle *p = new TParticle(
236  HEPEVT.idhep[i],
237  HEPEVT.isthep[i],
238  HEPEVT.jmohep[i][0]-1,
239  HEPEVT.jmohep[i][1]-1,
240  HEPEVT.jdahep[i][0]-1,
241  HEPEVT.jdahep[i][1]-1,
242  HEPEVT.phep[i][0],
243  HEPEVT.phep[i][1],
244  HEPEVT.phep[i][2],
245  HEPEVT.phep[i][3],
246  HEPEVT.vhep[i][0],
247  HEPEVT.vhep[i][1],
248  HEPEVT.vhep[i][2],
249  HEPEVT.vhep[i][3]);
250  fParticles->Add(p);
251  }
252  }
253  return fParticles;
254 }
255 
256 ////////////////////////////////////////////////////////////////////////////////
257 ///
258 /// It reads the /HEPEVT/ common block which has been filled by the
259 /// GenerateEvent method. If the event generator does not use the
260 /// HEPEVT common block, This routine has to be overloaded by the
261 /// subclasses.
262 ///
263 /// The function loops on the generated particles and store them in
264 /// the TClonesArray pointed by the argument particles. The default
265 /// action is to store only the stable particles (ISTHEP = 1) This can
266 /// be demanded explicitly by setting the option = "Final" If the
267 /// option = "All", all the particles are stored.
268 ///
269 
270 Int_t TGenerator::ImportParticles(TClonesArray *particles, Option_t *option)
271 {
272  if (particles == 0) return 0;
273  TClonesArray &clonesParticles = *particles;
274  clonesParticles.Clear();
275  Int_t numpart = HEPEVT.nhep;
276  if (!strcmp(option,"") || !strcmp(option,"Final")) {
277  for (Int_t i = 0; i<numpart; i++) {
278  if (HEPEVT.isthep[i] == 1) {
279 //
280 // Use the common block values for the TParticle constructor
281 //
282  new(clonesParticles[i]) TParticle(
283  HEPEVT.idhep[i],
284  HEPEVT.isthep[i],
285  HEPEVT.jmohep[i][0]-1,
286  HEPEVT.jmohep[i][1]-1,
287  HEPEVT.jdahep[i][0]-1,
288  HEPEVT.jdahep[i][1]-1,
289  HEPEVT.phep[i][0],
290  HEPEVT.phep[i][1],
291  HEPEVT.phep[i][2],
292  HEPEVT.phep[i][3],
293  HEPEVT.vhep[i][0],
294  HEPEVT.vhep[i][1],
295  HEPEVT.vhep[i][2],
296  HEPEVT.vhep[i][3]);
297  }
298  }
299  } else if (!strcmp(option,"All")) {
300  for (Int_t i = 0; i<numpart; i++) {
301  new(clonesParticles[i]) TParticle(
302  HEPEVT.idhep[i],
303  HEPEVT.isthep[i],
304  HEPEVT.jmohep[i][0]-1,
305  HEPEVT.jmohep[i][1]-1,
306  HEPEVT.jdahep[i][0]-1,
307  HEPEVT.jdahep[i][1]-1,
308  HEPEVT.phep[i][0],
309  HEPEVT.phep[i][1],
310  HEPEVT.phep[i][2],
311  HEPEVT.phep[i][3],
312  HEPEVT.vhep[i][0],
313  HEPEVT.vhep[i][1],
314  HEPEVT.vhep[i][2],
315  HEPEVT.vhep[i][3]);
316  }
317  }
318  return numpart;
319 }
320 
321 ////////////////////////////////////////////////////////////////////////////////
322 ///browse generator
323 
324 void TGenerator::Browse(TBrowser *)
325 {
326  Draw();
327  gPad->Update();
328 }
329 
330 ////////////////////////////////////////////////////////////////////////////////
331 /// Compute distance from point px,py to objects in event
332 ///
333 
334 Int_t TGenerator::DistancetoPrimitive(Int_t px, Int_t py)
335 {
336  const Int_t big = 9999;
337  const Int_t inview = 0;
338  Int_t dist = big;
339  if (px > 50 && py > 50) dist = inview;
340  return dist;
341 }
342 
343 ////////////////////////////////////////////////////////////////////////////////
344 ///
345 /// Insert one event in the pad list
346 ///
347 
348 void TGenerator::Draw(Option_t *option)
349 {
350  // Create a default canvas if a canvas does not exist
351  if (!gPad) {
352  gROOT->MakeDefCanvas();
353  if (gPad->GetVirtCanvas())
354  gPad->GetVirtCanvas()->SetFillColor(13);
355  }
356 
357  static Float_t rbox = 1000;
358  Float_t rmin[3],rmax[3];
359  TView *view = gPad->GetView();
360  if (!strstr(option,"same")) {
361  if (view) { view->GetRange(rmin,rmax); rbox = rmax[2];}
362  gPad->Clear();
363  }
364 
365  AppendPad(option);
366 
367  view = gPad->GetView();
368  // compute 3D view
369  if (view) {
370  view->GetRange(rmin,rmax);
371  rbox = rmax[2];
372  } else {
373  view = TView::CreateView(1,0,0);
374  if (view) view->SetRange(-rbox,-rbox,-rbox, rbox,rbox,rbox );
375  }
376  const Int_t kColorProton = 4;
377  const Int_t kColorNeutron = 5;
378  const Int_t kColorAntiProton= 3;
379  const Int_t kColorPionPlus = 6;
380  const Int_t kColorPionMinus = 2;
381  const Int_t kColorKaons = 7;
382  const Int_t kColorElectrons = 0;
383  const Int_t kColorGamma = 18;
384 
385  Int_t nProtons = 0;
386  Int_t nNeutrons = 0;
387  Int_t nAntiProtons= 0;
388  Int_t nPionPlus = 0;
389  Int_t nPionMinus = 0;
390  Int_t nKaons = 0;
391  Int_t nElectrons = 0;
392  Int_t nGammas = 0;
393 
394  Int_t ntracks = fParticles->GetEntriesFast();
395  Int_t i,lwidth,color,lstyle;
396  TParticlePDG *ap;
397  TParticle *p;
398  const char *name;
399  Double_t etot,vx,vy,vz;
400  Int_t ninvol = 0;
401  for (i=0;i<ntracks;i++) {
402  p = (TParticle*)fParticles->UncheckedAt(i);
403  if(!p) continue;
404  ap = (TParticlePDG*)p->GetPDG();
405  vx = p->Vx();
406  vy = p->Vy();
407  vz = p->Vz();
408  if (vx*vx+vy*vy+vz*vz > rbox*rbox) continue;
409  Float_t pt = p->Pt();
410  if (pt < fPtCut) continue;
411  etot = p->Energy();
412  if (etot > 0.1) lwidth = Int_t(6*TMath::Log10(etot));
413  else lwidth = 1;
414  if (lwidth < 1) lwidth = 1;
415  lstyle = 1;
416  color = 0;
417  name = ap->GetName();
418  if (!strcmp(name,"n")) { if (!fShowNeutrons) continue;
419  color = kColorNeutron; nNeutrons++;}
420  if (!strcmp(name,"p")) { color = kColorProton; nProtons++;}
421  if (!strcmp(name,"p bar")) { color = kColorAntiProton; nAntiProtons++;}
422  if (!strcmp(name,"pi+")) { color = kColorPionPlus; nPionPlus++;}
423  if (!strcmp(name,"pi-")) { color = kColorPionMinus; nPionMinus++;}
424  if (!strcmp(name,"e+")) { color = kColorElectrons; nElectrons++;}
425  if (!strcmp(name,"e-")) { color = kColorElectrons; nElectrons++;}
426  if (!strcmp(name,"gamma")) { color = kColorGamma; nGammas++; lstyle = 3; }
427  if ( strstr(name,"K")) { color = kColorKaons; nKaons++;}
428  p->SetLineColor(color);
429  p->SetLineStyle(lstyle);
430  p->SetLineWidth(lwidth);
431  p->AppendPad();
432  ninvol++;
433  }
434 
435  // event title
436  TPaveText *pt = new TPaveText(-0.94,0.85,-0.25,0.98,"br");
437  pt->AddText((char*)GetName());
438  pt->AddText((char*)GetTitle());
439  pt->SetFillColor(42);
440  pt->Draw();
441 
442  // Annotate color codes
443  Int_t tcolor = 5;
444  if (gPad->GetFillColor() == 10) tcolor = 4;
445  TText *text = new TText(-0.95,-0.47,"Particles");
446  text->SetTextAlign(12);
447  text->SetTextSize(0.025);
448  text->SetTextColor(tcolor);
449  text->Draw();
450  text->SetTextColor(kColorGamma); text->DrawText(-0.95,-0.52,"(on screen)");
451  text->SetTextColor(kColorGamma); text->DrawText(-0.95,-0.57,"Gamma");
452  text->SetTextColor(kColorProton); text->DrawText(-0.95,-0.62,"Proton");
453  text->SetTextColor(kColorNeutron); text->DrawText(-0.95,-0.67,"Neutron");
454  text->SetTextColor(kColorAntiProton); text->DrawText(-0.95,-0.72,"AntiProton");
455  text->SetTextColor(kColorPionPlus); text->DrawText(-0.95,-0.77,"Pion +");
456  text->SetTextColor(kColorPionMinus); text->DrawText(-0.95,-0.82,"Pion -");
457  text->SetTextColor(kColorKaons); text->DrawText(-0.95,-0.87,"Kaons");
458  text->SetTextColor(kColorElectrons); text->DrawText(-0.95,-0.92,"Electrons,etc.");
459 
460  text->SetTextColor(tcolor);
461  text->SetTextAlign(32);
462  char tcount[32];
463  snprintf(tcount,12,"%d",ntracks); text->DrawText(-0.55,-0.47,tcount);
464  snprintf(tcount,12,"%d",ninvol); text->DrawText(-0.55,-0.52,tcount);
465  snprintf(tcount,12,"%d",nGammas); text->DrawText(-0.55,-0.57,tcount);
466  snprintf(tcount,12,"%d",nProtons); text->DrawText(-0.55,-0.62,tcount);
467  snprintf(tcount,12,"%d",nNeutrons); text->DrawText(-0.55,-0.67,tcount);
468  snprintf(tcount,12,"%d",nAntiProtons); text->DrawText(-0.55,-0.72,tcount);
469  snprintf(tcount,12,"%d",nPionPlus); text->DrawText(-0.55,-0.77,tcount);
470  snprintf(tcount,12,"%d",nPionMinus); text->DrawText(-0.55,-0.82,tcount);
471  snprintf(tcount,12,"%d",nKaons); text->DrawText(-0.55,-0.87,tcount);
472  snprintf(tcount,12,"%d",nElectrons); text->DrawText(-0.55,-0.92,tcount);
473 
474  text->SetTextAlign(12);
475  if (nPionPlus+nPionMinus) {
476  snprintf(tcount,31,"Protons/Pions= %4f",Float_t(nProtons)/Float_t(nPionPlus+nPionMinus));
477  } else {
478  strlcpy(tcount,"Protons/Pions= inf",31);
479  }
480  text->DrawText(-0.45,-0.92,tcount);
481 
482  if (nPionPlus+nPionMinus) {
483  snprintf(tcount,31,"Kaons/Pions= %4f",Float_t(nKaons)/Float_t(nPionPlus+nPionMinus));
484  } else {
485  strlcpy(tcount,"Kaons/Pions= inf",31);
486  }
487  text->DrawText(0.30,-0.92,tcount);
488 }
489 
490 
491 ////////////////////////////////////////////////////////////////////////////////
492 /// Execute action corresponding to one event
493 ///
494 
495 void TGenerator::ExecuteEvent(Int_t event, Int_t px, Int_t py)
496 {
497  if (gPad->GetView()) {
498  gPad->GetView()->ExecuteRotateView(event, px, py);
499  return;
500  }
501 }
502 
503 ////////////////////////////////////////////////////////////////////////////////
504 /// Return the number of particles in the stack
505 
506 Int_t TGenerator::GetNumberOfParticles() const
507 {
508  return fParticles->GetLast()+1;
509 }
510 
511 ////////////////////////////////////////////////////////////////////////////////
512 /// Returns pointer to primary number i;
513 ///
514 
515 TParticle *TGenerator::GetParticle(Int_t i) const
516 {
517  if (!fParticles) return 0;
518  Int_t n = fParticles->GetLast();
519  if (i < 0 || i > n) return 0;
520  return (TParticle*)fParticles->UncheckedAt(i);
521 }
522 
523 ////////////////////////////////////////////////////////////////////////////////
524 ///
525 /// Paint one event
526 ///
527 
528 void TGenerator::Paint(Option_t *)
529 {
530 }
531 
532 ////////////////////////////////////////////////////////////////////////////////
533 ///
534 /// Set Pt threshold below which primaries are not drawn
535 ///
536 
537 void TGenerator::SetPtCut(Float_t ptcut)
538 {
539  fPtCut = ptcut;
540  Draw();
541  gPad->Update();
542 }
543 
544 ////////////////////////////////////////////////////////////////////////////////
545 ///
546 /// Set lower and upper values of the view range
547 ///
548 
549 void TGenerator::SetViewRadius(Float_t rbox)
550 {
551  SetViewRange(-rbox,-rbox,-rbox,rbox,rbox,rbox);
552 }
553 
554 ////////////////////////////////////////////////////////////////////////////////
555 ///
556 /// Set lower and upper values of the view range
557 ///
558 
559 void TGenerator::SetViewRange(Float_t xmin, Float_t ymin, Float_t zmin, Float_t xmax, Float_t ymax, Float_t zmax)
560 {
561  TView *view = gPad->GetView();
562  if (!view) return;
563  view->SetRange(xmin,ymin,zmin,xmax,ymax,zmax);
564 
565  Draw();
566  gPad->Update();
567 }
568 
569 ////////////////////////////////////////////////////////////////////////////////
570 ///
571 /// Set flag to display or not neutrons
572 ///
573 
574 void TGenerator::ShowNeutrons(Bool_t show)
575 {
576  fShowNeutrons = show;
577  Draw();
578  gPad->Update();
579 }