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