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