added stuff
[u/mrichter/AliRoot.git] / EVGEN / AliGenFLUKAsource.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 /* $Id$ */
17
18 // Read background particles from a FLUKA boundary source file
19 // This is a very special generator that works for background studies for the muon-spectrometer.
20 // The input files come from FLUKA simulations.
21 // Files can be chained. 
22 // Author: andreas.morsch@cern.ch
23
24 #include "AliGenFLUKAsource.h"
25 #include <stdlib.h>
26
27 #include <TDatabasePDG.h>
28 #include <TPDGCode.h>
29 #include <RVersion.h>
30 #include <TChain.h>
31 #include <TFile.h>
32 #include <TTree.h>
33 #include <TVirtualMC.h>
34
35 #include "AliRun.h"
36
37 ClassImp(AliGenFLUKAsource)
38
39 AliGenFLUKAsource::AliGenFLUKAsource()
40     :AliGenerator(-1), 
41      fIkine(6),
42      fAgeMax(1.e-5), 
43      fAddWeight(1.),
44      fZshift(0.),
45      fFrac(0.),
46      fSourceId(-1),
47      fFileName(0),
48      fTreeChain(0),
49      fTreeFluka(0),
50      fIp(0.),
51      fIpp(0.),
52      fXi(0.),
53      fYi(0.),
54      fZi(0.),
55      fPx(0.),
56      fPy(0.),
57      fPz(0.),
58      fEkin(0.),
59      fZv(0.),
60      fRv(0.),
61      fItra(0.),
62      fIgas(0.),
63      fWgt(0.),
64      fEtag(0.),
65      fPtg(0.),
66      fAge(0.)
67 {
68     // Constructor
69     fName="FLUKA";
70     fTitle="FLUKA Boundary Source";
71     // Read in all particle types by default
72     // Set maximum admitted age of particles to 1.e-05 by default 
73     // Do not add weight
74     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads 
75     // from  specific boundary source to the chamber (fZshift=4.5;),else there 
76     // is no need to shift as it reads boundary source for the whole volume of 
77     // the Muon Arm; the default value corresponds to boundary source for the
78     // whole volume of the MUON Arm 
79     // Set the default file 
80     fTreeChain = new TChain("h1");
81 //
82 //  Read all particles
83     fNpart=-1;
84 }
85
86 AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
87     :AliGenerator(npart), 
88      fIkine(6),
89      fAgeMax(1.e-5), 
90      fAddWeight(1.),
91      fZshift(0.),
92      fFrac(0.),
93      fSourceId(-1),
94      fFileName(""),
95      fTreeChain(new TChain("h1")),
96      fTreeFluka(0),
97      fIp(0.),
98      fIpp(0.),
99      fXi(0.),
100      fYi(0.),
101      fZi(0.),
102      fPx(0.),
103      fPy(0.),
104      fPz(0.),
105      fEkin(0.),
106      fZv(0.),
107      fRv(0.),
108      fItra(0.),
109      fIgas(0.),
110      fWgt(0.),
111      fEtag(0.),
112      fPtg(0.),
113      fAge(0.)
114 {
115     // Constructor
116     fName  = "FLUKA";
117     fTitle = "FLUKA Boundary Source";
118 }
119
120 AliGenFLUKAsource::AliGenFLUKAsource(const AliGenFLUKAsource & FLUKAsource):
121     AliGenerator(FLUKAsource), 
122     fIkine(6),
123     fAgeMax(1.e-5), 
124     fAddWeight(1.),
125     fZshift(0.),
126     fFrac(0.),
127     fSourceId(-1),
128     fFileName(0),
129     fTreeChain(0),
130     fTreeFluka(0),
131     fIp(0.),
132     fIpp(0.),
133     fXi(0.),
134     fYi(0.),
135     fZi(0.),
136     fPx(0.),
137     fPy(0.),
138     fPz(0.),
139     fEkin(0.),
140     fZv(0.),
141     fRv(0.),
142     fItra(0.),
143     fIgas(0.),
144     fWgt(0.),
145     fEtag(0.),
146     fPtg(0.),
147     fAge(0.)
148 {
149 // Copy constructor
150     FLUKAsource.Copy(*this);
151 }
152
153
154 //____________________________________________________________
155 AliGenFLUKAsource::~AliGenFLUKAsource()
156 {
157 // Destructor
158  if (fTreeFluka) delete fTreeFluka;
159  if (fTreeChain) delete fTreeChain;
160 }
161
162 void AliGenFLUKAsource::AddFile(const Text_t *filename)
163 {
164 // Add a file to the chain
165     fTreeChain->Add(filename);
166     
167 }
168
169
170 //____________________________________________________________
171 void AliGenFLUKAsource::FlukaInit() 
172 {
173 // Set branch addresses of data entries
174     TChain *h2=fTreeChain;
175 //Set branch addresses
176     h2->SetBranchAddress("Ip",&fIp);
177     h2->SetBranchAddress("Ipp",&fIpp);
178     h2->SetBranchAddress("Xi",&fXi);
179     h2->SetBranchAddress("Yi",&fYi);
180     h2->SetBranchAddress("Zi",&fZi);
181     h2->SetBranchAddress("Px",&fPx);
182     h2->SetBranchAddress("Py",&fPy);
183     h2->SetBranchAddress("Pz",&fPz);
184     h2->SetBranchAddress("Ekin",&fEkin);
185     h2->SetBranchAddress("Zv",&fZv);
186     h2->SetBranchAddress("Rv",&fRv);
187     h2->SetBranchAddress("Itra",&fItra);
188     h2->SetBranchAddress("Igas",&fIgas);
189     h2->SetBranchAddress("Wgt",&fWgt);
190     h2->SetBranchAddress("Etag",&fEtag);
191     h2->SetBranchAddress("Ptg",&fPtg);
192     h2->SetBranchAddress("Age",&fAge);
193 }
194
195 //____________________________________________________________
196 void AliGenFLUKAsource::Generate()
197 {
198 // Generate one event 
199
200     const Int_t kIfluge[28]={kProton, kProtonBar, kElectron, kPositron,
201                           kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
202                           kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
203                           kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
204                           kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
205                           0,kNuMu,kNuMuBar};
206     Float_t polar[3]= {0,0,0};
207   //
208     Float_t origin[3];
209     Float_t p[3];
210     Float_t prwn;
211     Float_t wgt, fwgt;
212     Float_t phi;
213     Float_t amass;
214     Int_t iwgt;
215     Int_t i, j, part, nt;
216     static Int_t irwn=0;
217     //
218     Float_t random[2];
219   //
220     FlukaInit();
221     TChain *h2=fTreeChain;
222     Int_t nentries = (Int_t) h2->GetEntries();
223     if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
224     
225
226   // loop over number of particles
227     Int_t nb=0;
228     Int_t ev=gMC->CurrentEvent();
229     for (i=0; i<fNpart;i++) {
230         Int_t entry=fNpart*(ev-1)+i; 
231         nb = (Int_t)h2->GetEvent(entry); 
232         if (irwn > nentries) {
233             printf("No more entries in the FLUKA boundary source file\n");
234             TFile *pFile=0;
235             // Get AliRun object or create it 
236             if (!gAlice) {
237                 gAlice = (AliRun*)pFile->Get("gAlice");
238                 if (gAlice) printf("AliRun object found on file\n");
239                 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
240             }
241             TTree *fAli=gAlice->TreeK();
242             if (fAli) pFile =fAli->GetCurrentFile();
243             pFile->cd();
244             printf("Generate - I'm out \n");
245             return;
246         }   
247         
248         Int_t ifip = Int_t(fIp);
249         
250
251         if (fSourceId != -1 && fIgas !=fSourceId) {
252             irwn++;
253             continue;
254         }
255         
256         if (ifip > 28 || ifip < 0) {
257             irwn++;
258             continue;
259         }
260         
261         if ((ifip != fIkine && fIkine != kAll && fIkine != kCharged 
262              && fIkine != 10) || fAge > fAgeMax){
263             irwn++;
264             continue;
265         } else if (fIkine == kCharged) {
266             if (ifip == 7 || ifip == 8 || fAge > fAgeMax) { 
267                 irwn++;
268                 continue;
269             }
270         } else if (fIkine == kNoNeutron) {
271             if (ifip == 8 || fAge > fAgeMax) { 
272                 irwn++;
273                 continue;
274             }
275         }
276     
277
278         irwn++;
279 //
280 // PDG code from FLUKA particle type (ifip)
281         part=kIfluge[int(ifip)-1];      
282 //
283 // Calculate momentum from kinetic energy and mass of the particle
284 #if ROOT_VERSION_CODE > 197895
285         amass = gMC->ParticleMass(part);
286 #else
287         amass = (TDatabasePDG::Instance())->GetParticle(part)->Mass();
288 #endif
289         prwn=fEkin*sqrt(1. + 2.*amass/fEkin);
290
291
292         origin[0]=fYi;
293         origin[1]=fXi;
294         origin[2]=fZi;
295         
296         p[0]=fPy*prwn;
297         p[1]=fPx*prwn;
298         p[2]=fPz*prwn;
299
300         //handle particle weight correctly
301         wgt = (part == 13) ? fWgt*fAddWeight : fWgt;
302         iwgt=Int_t(wgt);
303         fwgt=wgt-Float_t(iwgt);
304         Rndm(random,2);
305         if (random[0] < fwgt) iwgt++;
306         if (part==1 && iwgt>100) iwgt=100;
307         Int_t nstack=0;
308         for (j=0; j<iwgt; j++) {
309             PushTrack(fTrackIt,-1,part,p,origin,polar,fAge,kPPrimary,nt);
310             Rndm(random,2);
311             phi=2*random[1]*TMath::Pi();
312             Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
313             Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
314             p[0]=pn1;
315             p[1]=pn2;
316             Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
317             Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
318             origin[0]=on1;
319             origin[1]=on2;
320             nstack++;
321         }
322         if (nstack == 0) continue;
323     }
324     
325     TFile *pFile=0;
326 // Get AliRun object or create it 
327     if (!gAlice) {
328         gAlice = (AliRun*)pFile->Get("gAlice");
329         if (gAlice) printf("AliRun object found on file\n");
330         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
331     }
332     TTree *fAli=gAlice->TreeK();
333     if (fAli) pFile =fAli->GetCurrentFile();
334     pFile->cd();
335 }
336
337
338 AliGenFLUKAsource& AliGenFLUKAsource::operator=(const  AliGenFLUKAsource& rhs)
339 {
340 // Assignment operator
341     rhs.Copy(*this);
342     return (*this);
343 }
344
345
346 void AliGenFLUKAsource::Copy(TObject &) const
347 {
348     Fatal("Copy","Not implemented!\n");
349 }
350
351
352
353
354