Adding Generation via an external file.
[u/mrichter/AliRoot.git] / EVGEN / AliGenExtFile.cxx
CommitLineData
693caace 1#include "AliGenExtFile.h"
2#include "AliGenMUONlib.h"
3#include "AliMC.h"
4#include "AliRun.h"
5#include <TDirectory.h>
6#include <TFile.h>
7#include <TTree.h>
8#include <stdlib.h>
9 ClassImp(AliGenExtFile)
10 AliGenExtFile::AliGenExtFile()
11 :AliGenerator(-1)
12{
13 //
14 fName="ExtFile";
15 fTitle="Primaries from ext. File";
16 fFileName="dtujet93.root";
17 fTreeNtuple=0;
18 fNcurrent=0;
19//
20// Read all particles
21 fNpart=-1;
22}
23
24AliGenExtFile::AliGenExtFile(Int_t npart)
25 :AliGenerator(npart)
26{
27 //
28 fName="ExtFile";
29 fTitle="Primaries from ext. File";
30 fFileName="dtujet93.root";
31 fTreeNtuple=0;
32 fNcurrent=0;
33}
34
35//____________________________________________________________
36AliGenExtFile::~AliGenExtFile()
37{
38 delete fTreeNtuple;
39}
40
41//____________________________________________________________
42void AliGenExtFile::NtupleInit()
43{
44//
45// reset the existing file environment and open a new root file if
46// the pointer to the Fluka tree is null
47
48 TFile *File=0;
49 if (fTreeNtuple==0) {
50 if (!File) {
51 File = new TFile(fFileName);
52 File->cd();
53 cout<<"I have opened "<<fFileName<<" file "<<endl;
54 }
55// get the tree address in the Fluka boundary source file
56 fTreeNtuple = (TTree*)gDirectory->Get("h888");
57 } else {
58 File = fTreeNtuple->GetCurrentFile();
59 File->cd();
60 }
61
62 TTree *h2=fTreeNtuple;
63//Set branch addresses
64//Set branch addresses
65 h2->SetBranchAddress("Nihead",&Nihead);
66 h2->SetBranchAddress("Ihead",Ihead);
67 h2->SetBranchAddress("Nrhead",&Nrhead);
68 h2->SetBranchAddress("Rhead",Rhead);
69 h2->SetBranchAddress("Idpart",&Idpart);
70 h2->SetBranchAddress("Theta",&Theta);
71 h2->SetBranchAddress("Phi",&Phi);
72 h2->SetBranchAddress("P",&P);
73 h2->SetBranchAddress("E",&E);
74}
75
76
77//____________________________________________________________
78void AliGenExtFile::Generate()
79{
80
81 AliMC* pMC = AliMC::GetMC();
82
83 Float_t polar[3]= {0,0,0};
84 //
85 Float_t origin[3]={0,0,0};
86 Float_t p[3];
87 Float_t random[6];
88 Float_t prwn;
89 char name[100];
90 Float_t amass, charge, tlife;
91 Int_t itrtyp;
92 Int_t i, j, nt, Ntracks;
93 //
94 NtupleInit();
95 TTree *h2=fTreeNtuple;
96 Int_t nentries = (Int_t) h2->GetEntries();
97 // loop over number of particles
98 Int_t nb = (Int_t)h2->GetEvent(fNcurrent);
99 Int_t i5=Ihead[4];
100 Int_t i6=Ihead[5];
101
102 for (j=0;j<3;j++) origin[j]=fOrigin[j];
103 if(fVertexSmear==perEvent) {
104 pMC->Rndm(random,6);
105 for (j=0;j<3;j++) {
106 origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
107 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
108 }
109 }
110
111 if (fNcurrent >= nentries) {
112 printf("\n No more entries !!! !\n");
113 return;
114 }
115
116
117 if (i5==0) {
118 printf("\n This should never happen !\n");
119 } else {
120 printf("\n Next event contains %d tracks! \n", i6);
121 Ntracks=i6;
122 }
123 for (i=0; i<Ntracks; i++) {
124
125 pMC->Gfpart(Idpart, name, itrtyp,amass, charge, tlife);
126 prwn=sqrt((E+amass)*(E-amass));
127
128 Theta *= TMath::Pi()/180.;
129 Phi = (Phi-180)*TMath::Pi()/180.;
130 if(Theta<fThetaMin || Theta>fThetaMax ||
131 Phi<fPhiMin || Phi>fPhiMax ||
132 prwn<fPMin || prwn>fPMax)
133 {
134 ;
135 } else {
136 p[0]=prwn*TMath::Sin(Theta)*TMath::Cos(Phi);
137 p[1]=prwn*TMath::Sin(Theta)*TMath::Sin(Phi);
138 p[2]=prwn*TMath::Cos(Theta);
139
140 if(fVertexSmear==perTrack) {
141 pMC->Rndm(random,6);
142 for (j=0;j<3;j++) {
143 origin[j]=fOrigin[j]
144 +fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
145 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
146 }
147 }
148 gAlice->SetTrack(fTrackIt,-1,Idpart,p,origin,polar,0,"Primary",nt);
149 }
150 fNcurrent++;
151 nb = (Int_t)h2->GetEvent(fNcurrent);
152 }
153
154 TFile *File=0;
155// Get AliRun object or create it
156 if (!gAlice) {
157 gAlice = (AliRun*)File->Get("gAlice");
158 if (gAlice) printf("AliRun object found on file\n");
159 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
160 }
161 TTree *fAli=gAlice->TreeK();
162 if (fAli) File =fAli->GetCurrentFile();
163 File->cd();
164}
165
166
167
168
169
170
171
172
173