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