5a5ade01edca76f3e6f478d37e3e0e3bc8fb9008
[u/mrichter/AliRoot.git] / EVGEN / AliGenFLUKAsource.cxx
1 #include "AliGenFLUKAsource.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(AliGenFLUKAsource)
10      AliGenFLUKAsource::AliGenFLUKAsource()
11          :AliGenerator(-1)
12 {
13     //
14     fName="FLUKA";
15     fTitle="FLUKA Boundary Source";
16     // Read in all particle types by default
17     fIkine=6;
18     // Set maximum admitted age of particles to 1.e-05 by default 
19     fAgeMax=1.e-05;
20     // Do not add weight
21     fAddWeight=1.;
22     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads 
23     // from  specific boundary source to the chamber (fZshift=4.5;),else there 
24     // is no need to shift as it reads boundary source for the whole volume of 
25     // the Muon Arm; the default value corresponds to boundary source for the
26     // whole volume of the MUON Arm 
27     fZshift=0;
28     // Set the default file 
29     fFileName="flukasource.root";
30
31     fTreeFluka=0;
32 //
33 //  Read all particles
34     fNpart=-1;
35     
36 }
37
38 AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
39     :AliGenerator(npart)
40 {
41     //
42     fName="FLUKA";
43     fTitle="FLUKA Boundary Source";
44     // Read in all particle types by default
45     fIkine=6;
46     // Set maximum admitted age of particles to 1.e-05 by default 
47     fAgeMax=1.e-05;
48     // Do not add weight
49     fAddWeight=1.;
50     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads 
51     // from  specific boundary source to the chamber (fZshift=4.5;),else there 
52     // is no need to shift as it reads boundary source for the whole volume of 
53     // the Muon Arm; the default value corresponds to boundary source for the
54     // whole volume of the MUON Arm 
55     fZshift=0;
56     // Set the default file 
57     fFileName="flukasource.root";
58
59     fTreeFluka=0;
60
61 }
62
63 //____________________________________________________________
64 AliGenFLUKAsource::~AliGenFLUKAsource()
65 {
66     delete fTreeFluka;
67 }
68
69 //____________________________________________________________
70 void AliGenFLUKAsource::FlukaInit() 
71 {
72 //
73 // reset the existing file environment and open a new root file if
74 // the pointer to the Fluka tree is null
75     
76     TFile *File=0;
77     if (fTreeFluka==0) {
78         if (!File) {
79             File = new TFile(fFileName);
80             File->cd();
81             cout<<"I have opened "<<fFileName<<" file "<<endl;
82         }
83 // get the tree address in the Fluka boundary source file
84         fTreeFluka = (TTree*)gDirectory->Get("h1");
85     } else {
86         File = fTreeFluka->GetCurrentFile();
87         File->cd();
88     }
89
90     TTree *h2=fTreeFluka;
91 //Set branch addresses
92     h2->SetBranchAddress("Ip",&Ip);
93     h2->SetBranchAddress("Ipp",&Ipp);
94     h2->SetBranchAddress("Xi",&Xi);
95     h2->SetBranchAddress("Yi",&Yi);
96     h2->SetBranchAddress("Zi",&Zi);
97     h2->SetBranchAddress("Px",&Px);
98     h2->SetBranchAddress("Py",&Py);
99     h2->SetBranchAddress("Pz",&Pz);
100     h2->SetBranchAddress("Ekin",&Ekin);
101     h2->SetBranchAddress("Zv",&Zv);
102     h2->SetBranchAddress("Rv",&Rv);
103     h2->SetBranchAddress("Itra",&Itra);
104     h2->SetBranchAddress("Igas",&Igas);
105     h2->SetBranchAddress("Wgt",&Wgt);
106     h2->SetBranchAddress("Etag",&Etag);
107     h2->SetBranchAddress("Ptg",&Ptg);
108     h2->SetBranchAddress("Age",&Age);
109 }
110
111 //____________________________________________________________
112 void AliGenFLUKAsource::Generate()
113 {
114
115   AliMC* pMC = AliMC::GetMC();
116
117   const Int_t ifluge[28]={14, 15, 3, 2, 4, 4, 1, 13, 25, 5, 6, 10, 8, 9,
118                          11, 12,18, 26, 16, 21, 19, 20, 7, 16,16,0,0,0};
119   Float_t polar[3]= {0,0,0};
120   //
121   Float_t origin[3];
122   Float_t p[3];
123   Float_t prwn;
124   Float_t wgt, fwgt;
125   Float_t phi;
126   char name[100];
127   Float_t amass, charge, tlife;
128   Int_t itrtyp;
129   Int_t iwgt;
130   Int_t i, j, part, nt;
131   static Int_t irwn;
132   //
133   Float_t random[2];
134   //
135   FlukaInit();
136   TTree *h2=fTreeFluka;
137   Int_t nentries = (Int_t) h2->GetEntries();
138   if (fNpart == -1) fNpart=nentries;
139   
140   // loop over number of particles
141   Int_t nb=0;
142   for (i=0; i<fNpart;i++) {
143     Int_t ev=pMC->CurrentEvent();
144     Int_t entry=fNpart*(ev-1)+i; 
145     nb = (Int_t)h2->GetEvent(entry); 
146     if (irwn > nentries) {
147       printf("No more entries in the FLUKA boundary source file\n");
148       return;
149     }   
150     if (Ip > 28 || Ip < 0) {
151       irwn++;
152       continue;
153     }
154  
155     if ((Ip != fIkine && fIkine != 6 && fIkine != 9) || Age > fAgeMax){
156        irwn++;
157        continue;
158     } else if (fIkine == 9) {
159        if (Ip == 7 || Ip == 8 || Age > fAgeMax) { 
160           irwn++;
161           continue;
162        }
163     }
164     
165
166     irwn++;
167     printf("\n Particle type: %f \n \n ", Ip);
168     
169     if (Ip ==7){
170        prwn=Ekin;
171        part=1;
172     } else if (Ip == 8) {
173        prwn=sqrt(Ekin*Ekin + 2.*0.93956563);
174        part=13;
175     } else {
176        part=ifluge[int(Ip)-1];
177        pMC->Gfpart(part, name, itrtyp,  
178                    amass, charge, tlife); 
179        prwn=sqrt(Ekin*Ekin + 2.*amass);
180     }
181     origin[0]=Xi;
182     origin[1]=Yi;
183     origin[2]=Zi;
184
185     p[0]=Px*prwn;
186     p[1]=Py*prwn;
187     p[2]=Pz*prwn;
188
189     //handle particle weight correctly
190     wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
191     iwgt=Int_t(wgt);
192     fwgt=wgt-Float_t(iwgt);
193     pMC->Rndm(random,2);
194     if (random[0] < fwgt) iwgt++;
195     if (part==1 && iwgt>100) iwgt=100;
196     Int_t nstack=0;
197     for (j=0; j<iwgt; j++) {
198         gAlice->SetTrack(1,-1,part,p,origin,polar,0,"Primary",nt);
199         pMC->Rndm(random,2);
200         phi=2*random[1]*TMath::Pi();
201         Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
202         Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
203         p[0]=pn1;
204         p[1]=pn2;
205         Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
206         Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
207         origin[0]=on1;
208         origin[1]=on2;
209         nstack++;
210     }
211     if (nstack == 0) continue;
212   }
213  
214     TFile *File=0;
215 // Get AliRun object or create it 
216     if (!gAlice) {
217         gAlice = (AliRun*)File->Get("gAlice");
218         if (gAlice) printf("AliRun object found on file\n");
219         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
220     }
221     TTree *fAli=gAlice->TreeK();
222     if (fAli) File =fAli->GetCurrentFile();
223     File->cd();
224 }
225
226
227
228
229
230
231
232
233