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