Introduction of the Copyright and cvs Log
[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 */
19
20 #include "AliGenFLUKAsource.h"
21 #include "AliGenMUONlib.h"
22 #include "AliMC.h"
23 #include "AliRun.h"
24 #include "AliPDG.h"
25
26 #include <TDirectory.h>
27 #include <TFile.h>
28 #include <TTree.h>
29 #include <stdlib.h>
30  ClassImp(AliGenFLUKAsource)
31      AliGenFLUKAsource::AliGenFLUKAsource()
32          :AliGenerator(-1)
33 {
34     //
35     fName="FLUKA";
36     fTitle="FLUKA Boundary Source";
37     // Read in all particle types by default
38     fIkine=6;
39     // Set maximum admitted age of particles to 1.e-05 by default 
40     fAgeMax=1.e-05;
41     // Do not add weight
42     fAddWeight=1.;
43     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads 
44     // from  specific boundary source to the chamber (fZshift=4.5;),else there 
45     // is no need to shift as it reads boundary source for the whole volume of 
46     // the Muon Arm; the default value corresponds to boundary source for the
47     // whole volume of the MUON Arm 
48     fZshift=0;
49     // Set the default file 
50     fFileName="flukasource.root";
51
52     fTreeFluka=0;
53 //
54 //  Read all particles
55     fNpart=-1;
56     
57 }
58
59 AliGenFLUKAsource::AliGenFLUKAsource(Int_t npart)
60     :AliGenerator(npart)
61 {
62     //
63     fName="FLUKA";
64     fTitle="FLUKA Boundary Source";
65     // Read in all particle types by default
66     fIkine=6;
67     // Set maximum admitted age of particles to 1.e-05 by default 
68     fAgeMax=1.e-05;
69     // Do not add weight
70     fAddWeight=1.;
71     // Shift the z-coordinate of the impact point by 4.5 cm only if it reads 
72     // from  specific boundary source to the chamber (fZshift=4.5;),else there 
73     // is no need to shift as it reads boundary source for the whole volume of 
74     // the Muon Arm; the default value corresponds to boundary source for the
75     // whole volume of the MUON Arm 
76     fZshift=0;
77     // Set the default file 
78     fFileName="flukasource.root";
79
80     fTreeFluka=0;
81
82 }
83
84 //____________________________________________________________
85 AliGenFLUKAsource::~AliGenFLUKAsource()
86 {
87     delete fTreeFluka;
88 }
89
90 //____________________________________________________________
91 void AliGenFLUKAsource::FlukaInit() 
92 {
93 //
94 // reset the existing file environment and open a new root file if
95 // the pointer to the Fluka tree is null
96     
97     TFile *File=0;
98     if (fTreeFluka==0) {
99         if (!File) {
100             File = new TFile(fFileName);
101             File->cd();
102             cout<<"I have opened "<<fFileName<<" file "<<endl;
103         }
104 // get the tree address in the Fluka boundary source file
105         fTreeFluka = (TTree*)gDirectory->Get("h1");
106     } else {
107         File = fTreeFluka->GetCurrentFile();
108         File->cd();
109     }
110
111     TTree *h2=fTreeFluka;
112 //Set branch addresses
113     h2->SetBranchAddress("Ip",&Ip);
114     h2->SetBranchAddress("Ipp",&Ipp);
115     h2->SetBranchAddress("Xi",&Xi);
116     h2->SetBranchAddress("Yi",&Yi);
117     h2->SetBranchAddress("Zi",&Zi);
118     h2->SetBranchAddress("Px",&Px);
119     h2->SetBranchAddress("Py",&Py);
120     h2->SetBranchAddress("Pz",&Pz);
121     h2->SetBranchAddress("Ekin",&Ekin);
122     h2->SetBranchAddress("Zv",&Zv);
123     h2->SetBranchAddress("Rv",&Rv);
124     h2->SetBranchAddress("Itra",&Itra);
125     h2->SetBranchAddress("Igas",&Igas);
126     h2->SetBranchAddress("Wgt",&Wgt);
127     h2->SetBranchAddress("Etag",&Etag);
128     h2->SetBranchAddress("Ptg",&Ptg);
129     h2->SetBranchAddress("Age",&Age);
130 }
131
132 //____________________________________________________________
133 void AliGenFLUKAsource::Generate()
134 {
135
136   const Int_t ifluge[28]={kProton, kProtonBar, kElectron, kPositron,
137                           kNuE, kNuEBar, kGamma, kNeutron, kNeutronBar,
138                           kMuonPlus, kMuonMinus, kK0Long , kPiPlus, kPiMinus,
139                           kKPlus, kKMinus, kLambda0, kLambda0Bar, kK0Short,
140                           kSigmaMinus, kSigmaPlus, kSigma0, kPi0, kK0, kK0Bar,
141                           0,kNuMu,kNuMuBar};
142   Float_t polar[3]= {0,0,0};
143   //
144   Float_t origin[3];
145   Float_t p[3];
146   Float_t prwn;
147   Float_t wgt, fwgt;
148   Float_t phi;
149   char name[100];
150   Float_t amass, charge, tlife;
151   Int_t itrtyp;
152   Int_t iwgt;
153   Int_t i, j, part, nt;
154   static Int_t irwn;
155   //
156   Float_t random[2];
157   //
158   FlukaInit();
159   TTree *h2=fTreeFluka;
160   Int_t nentries = (Int_t) h2->GetEntries();
161   if (fNpart == -1) fNpart=Int_t(nentries*fFrac);
162   
163   // loop over number of particles
164   Int_t nb=0;
165   for (i=0; i<fNpart;i++) {
166     Int_t ev=gMC->CurrentEvent();
167     Int_t entry=fNpart*(ev-1)+i; 
168     nb = (Int_t)h2->GetEvent(entry); 
169     if (irwn > nentries) {
170       printf("No more entries in the FLUKA boundary source file\n");
171       TFile *File=0;
172       // Get AliRun object or create it 
173       if (!gAlice) {
174         gAlice = (AliRun*)File->Get("gAlice");
175         if (gAlice) printf("AliRun object found on file\n");
176         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
177       }
178       TTree *fAli=gAlice->TreeK();
179       if (fAli) File =fAli->GetCurrentFile();
180       File->cd();
181       printf("Generate - I'm out \n");
182       return;
183     }   
184     if (Ip > 28 || Ip < 0) {
185       irwn++;
186       continue;
187     }
188  
189     if ((Ip != fIkine && fIkine != 6 && fIkine != 9) || Age > fAgeMax){
190        irwn++;
191        continue;
192     } else if (fIkine == 9) {
193        if (Ip == 7 || 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=sqrt(Ekin*Ekin + 2.*0.93956563);
208        part=13;
209     } else {
210        part=ifluge[int(Ip)-1];
211        gMC->Gfpart(part, name, itrtyp,  
212                    amass, charge, tlife); 
213        prwn=sqrt(Ekin*Ekin + 2.*amass);
214     }
215     origin[0]=Xi;
216     origin[1]=Yi;
217     origin[2]=Zi;
218
219     p[0]=Px*prwn;
220     p[1]=Py*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(1,-1,part,p,origin,polar,0,"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