Protection against div by 0 added.
[u/mrichter/AliRoot.git] / FASTSIM / fastGenEMCocktail.C
1 void fastGenEMCocktail(Int_t nev = 1, char* filename = "galice.root")
2 {
3 // load libraries
4
5   gSystem->Load("liblhapdf.so");      // Parton density functions
6   gSystem->Load("libpythia6.so");     // Pythia
7   gSystem->Load("libEG");
8   gSystem->Load("libEGPythia6");
9   gSystem->Load("libAliPythia6.so");  // ALICE specific implementations
10
11 // Runloader
12     
13   AliRunLoader* rl = AliRunLoader::Open("galice.root","FASTRUN","recreate");
14     
15   rl->SetCompressionLevel(2);
16   rl->SetNumberOfEventsPerFile(nev);
17   rl->LoadKinematics("RECREATE");
18   rl->MakeTree("E");
19   gAlice->SetRunLoader(rl);
20
21 // Create stack
22   rl->MakeStack();
23   AliStack* stack      = rl->Stack();
24  
25 // Header
26   AliHeader* header = rl->GetHeader();
27
28 // Create and Initialize Generator
29
30   AliGenEMCocktail *gener = new AliGenEMCocktail();
31
32 //=======================================================================
33 // Set External decayer
34   TVirtualMCDecayer *decayer = new AliDecayerPythia();
35
36   gener->SetNPart(1000);               // source multiplicity per event
37   gener->SetPtRange(0.,20.);
38   gener->SetYRange(-1.,1.);
39   gener->SetPhiRange(0., 360.);
40   gener->SetOrigin(0.,0.,0.); 
41   gener->SetSigma(0.,0.,0.);
42   gener->SetVertexSmear(kPerEvent);
43   gener->SetTrackingFlag(0);
44   gener->SetDecayMode(kElectronEM);    // select cocktail:
45                                        // kElectronEM   => single electron
46                                        // kDiElectronEM => electron-positron
47                                        // kGammaEM      => single photon
48   gener->SetDecayer(decayer);
49   gener->SetWeightingMode(kNonAnalog); // select weighting:
50                                        // kNonAnalog => weight ~ dN/dp_T
51                                        // kAnalog    => weight ~ 1
52   gener->CreateCocktail();
53   gener->Init();
54   gener->SetStack(stack);
55
56 //
57 //                        Event Loop
58 //
59   Int_t iev;
60      
61   for (iev = 0; iev < nev; iev++) {
62
63     printf("\n \n Event number %d \n \n", iev);
64         
65 // Initialize event
66     header->Reset(0,iev);
67     rl->SetEventNumber(iev);
68     stack->Reset();
69     rl->MakeTree("K");
70 //  stack->ConnectTree();
71     
72 // Generate event
73     gener->Generate();
74 // Analysis
75     Int_t npart = stack->GetNprimary();
76     printf("Analyse %d Particles\n", npart);
77     for (Int_t part=0; part<npart; part++) {
78       TParticle *MPart = stack->Particle(part);
79       Int_t mpart  = MPart->GetPdgCode();
80 //      printf("Particle %d\n", mpart);
81     }
82         
83 // Finish event
84     header->SetNprimary(stack->GetNprimary());
85     header->SetNtrack(stack->GetNtrack());  
86
87 // I/O
88 //      
89     stack->FinishEvent();
90     header->SetStack(stack);
91     rl->TreeE()->Fill();
92     rl->WriteKinematics("OVERWRITE");
93
94   } // event loop
95
96 //
97 //                         Termination
98 // Generator
99   gener->FinishRun();
100 // Write file
101   rl->WriteHeader("OVERWRITE");
102   gener->Write();
103   rl->Write();
104 }
105
106
107
108
109
110
111
112
113
114
115
116