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