New version from G.Martinez & A.Morsch
[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 /*
17 $Log$
18 Revision 1.6  1999/09/29 09:24:12  fca
19 Introduction of the Copyright and cvs Log
20
21 */
22
23 #include "AliGenFLUKAsource.h"
24 #include "AliGenMUONlib.h"
25 #include "AliMC.h"
26 #include "AliRun.h"
27 #include "AliPDG.h"
28 #include <TDirectory.h>
29 #include <TFile.h>
30 #include <TTree.h>
31 #include <stdlib.h>
32  ClassImp(AliGenFLUKAsource)
33      AliGenFLUKAsource::AliGenFLUKAsource()
34          :AliGenerator(-1)
35 {
36     //
37     fName="FLUKA";
38     fTitle="FLUKA Boundary Source";
39     // Read in all particle types by default
40     fIkine=6;
41     // Set maximum admitted age of particles to 1.e-05 by default 
42     fAgeMax=1.e-05;
43     // Do not add weight
44     fAddWeight=1.;
45     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads 
46     // from  specific boundary source to the chamber (fZshift=4.5;),else there 
47     // is no need to shift as it reads boundary source for the whole volume of 
48     // the Muon Arm; the default value corresponds to boundary source for the
49     // whole volume of the MUON Arm 
50     fZshift=0;
51     // Set the default file 
52     fFileName="flukasource.root";
53
54     fTreeFluka=0;
55     fTreeChain = new TChain("h1");
56 //
57 //  Read all particles
58     fNpart=-1;
59     
60 }
61
62 AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
63     :AliGenerator(npart)
64 {
65     //
66     fName="FLUKA";
67     fTitle="FLUKA Boundary Source";
68     // Read in all particle types by default
69     fIkine=6;
70     // Set maximum admitted age of particles to 1.e-05 by default 
71     fAgeMax=1.e-05;
72     // Do not add weight
73     fAddWeight=1.;
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     fZshift=0;
80     // Set the default file 
81     fFileName="flukasource.root";
82
83     fTreeFluka=0;
84     fTreeChain = new TChain("h1"); 
85 }
86
87 //____________________________________________________________
88 AliGenFLUKAsource::~AliGenFLUKAsource()
89 {
90  if (fTreeFluka) delete fTreeFluka;
91  if (fTreeChain) delete fTreeChain;
92  // if (fFileName)  delete fFileName;
93 }
94
95 void AliGenFLUKAsource::AddFile(const Text_t *filename)
96 {
97     fTreeChain->Add(filename);
98     
99 }
100
101
102 //____________________________________________________________
103 void AliGenFLUKAsource::FlukaInit() 
104 {
105     TChain *h2=fTreeChain;
106 //Set branch addresses
107     h2->SetBranchAddress("Ip",&Ip);
108     h2->SetBranchAddress("Ipp",&Ipp);
109     h2->SetBranchAddress("Xi",&Xi);
110     h2->SetBranchAddress("Yi",&Yi);
111     h2->SetBranchAddress("Zi",&Zi);
112     h2->SetBranchAddress("Px",&Px);
113     h2->SetBranchAddress("Py",&Py);
114     h2->SetBranchAddress("Pz",&Pz);
115     h2->SetBranchAddress("Ekin",&Ekin);
116     h2->SetBranchAddress("Zv",&Zv);
117     h2->SetBranchAddress("Rv",&Rv);
118     h2->SetBranchAddress("Itra",&Itra);
119     h2->SetBranchAddress("Igas",&Igas);
120     h2->SetBranchAddress("Wgt",&Wgt);
121     h2->SetBranchAddress("Etag",&Etag);
122     h2->SetBranchAddress("Ptg",&Ptg);
123     h2->SetBranchAddress("Age",&Age);
124 }
125
126 //____________________________________________________________
127 void AliGenFLUKAsource::Generate()
128
129     AliMC* gMC = AliMC::GetMC();
130
131     const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
132                           kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
133                           kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
134                           kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
135                           kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
136                           0,kNuMu,kNuMuBar};
137     Float_t polar[3]= {0,0,0};
138   //
139     Float_t origin[3];
140     Float_t p[3];
141     Float_t prwn;
142     Float_t wgt, fwgt;
143     Float_t phi;
144     char name[100];
145     Float_t amass, charge, tlife;
146     Int_t itrtyp;
147     Int_t iwgt;
148     Int_t i, j, part, nt;
149     static Int_t irwn;
150     //
151     Float_t random[2];
152   //
153     FlukaInit();
154     TChain *h2=fTreeChain;
155     Int_t nentries = (Int_t) h2->GetEntries();
156     if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
157   
158   // loop over number of particles
159     Int_t nb=0;
160     Int_t ev=gMC->CurrentEvent();
161     for (i=0; i<fNpart;i++) {
162         Int_t entry=fNpart*(ev-1)+i; 
163         nb = (Int_t)h2->GetEvent(entry); 
164         if (irwn > nentries) {
165             printf("No more entries in the FLUKA boundary source file\n");
166             TFile *File=0;
167             // Get AliRun object or create it 
168             if (!gAlice) {
169                 gAlice = (AliRun*)File->Get("gAlice");
170                 if (gAlice) printf("AliRun object found on file\n");
171                 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
172             }
173             TTree *fAli=gAlice->TreeK();
174             if (fAli) File =fAli->GetCurrentFile();
175             File->cd();
176             printf("Generate - I'm out \n");
177             return;
178         }   
179         if (Ip > 28 || Ip < 0) {
180             irwn++;
181             continue;
182         }
183         
184         if ((Ip != fIkine && fIkine != 6 && fIkine != 9 && fIkine != 10) || Age > fAgeMax){
185             irwn++;
186             continue;
187         } else if (fIkine == 9) {
188             if (Ip == 7 || Ip == 8 || Age > fAgeMax) { 
189                 irwn++;
190                 continue;
191             }
192         } else if (fIkine == 10) {
193             if (Ip == 8 || Age > fAgeMax) { 
194                 irwn++;
195                 continue;
196             }
197         }
198     
199
200         irwn++;
201         printf("\n Particle type: %f \n \n ", Ip);
202         
203         if (Ip ==7){
204             prwn=Ekin;
205             part=1;
206         } else if (Ip == 8) {
207             prwn=Ekin*sqrt(1. + 2.*0.93956563/Ekin);
208             part=13;
209         } else {
210             part=ifluge[int(Ip)-1];
211             gMC->Gfpart(part, name, itrtyp,  
212                         amass, charge, tlife); 
213             prwn=Ekin*sqrt(1. + 2.*amass/Ekin);
214         }
215         origin[0]=Yi;
216         origin[1]=Xi;
217         origin[2]=Zi;
218         
219         p[0]=Py*prwn;
220         p[1]=Px*prwn;
221         p[2]=Pz*prwn;
222
223         //handle particle weight correctly
224         wgt = (part == 13) ? Wgt*fAddWeight : Wgt;
225         iwgt=Int_t(wgt);
226         fwgt=wgt-Float_t(iwgt);
227         gMC->Rndm(random,2);
228         if (random[0] < fwgt) iwgt++;
229         if (part==1 && iwgt>100) iwgt=100;
230         Int_t nstack=0;
231         for (j=0; j<iwgt; j++) {
232             gAlice->SetTrack(fTrackIt,-1,part,p,origin,polar,Age,"Primary",nt);
233             gMC->Rndm(random,2);
234             phi=2*random[1]*TMath::Pi();
235             Float_t pn1=p[0]*TMath::Sin(phi) - p[1]*TMath::Cos(phi);
236             Float_t pn2=p[0]*TMath::Cos(phi) + p[1]*TMath::Sin(phi);
237             p[0]=pn1;
238             p[1]=pn2;
239             Float_t on1=origin[0]*TMath::Sin(phi)-origin[1]*TMath::Cos(phi);
240             Float_t on2=origin[0]*TMath::Cos(phi)+origin[1]*TMath::Sin(phi);
241             origin[0]=on1;
242             origin[1]=on2;
243             nstack++;
244         }
245         if (nstack == 0) continue;
246     }
247     
248     TFile *File=0;
249 // Get AliRun object or create it 
250     if (!gAlice) {
251         gAlice = (AliRun*)File->Get("gAlice");
252         if (gAlice) printf("AliRun object found on file\n");
253         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
254     }
255     TTree *fAli=gAlice->TreeK();
256     if (fAli) File =fAli->GetCurrentFile();
257     File->cd();
258 }
259
260
261
262
263
264
265
266
267