AliGenExtFile derives from AliGenMC. Generate() uses methods from AliGenMC (N. Carrer)
[u/mrichter/AliRoot.git] / EVGEN / AliGenExtFile.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 /*
17 $Log$
18 Revision 1.20  2002/03/22 09:43:28  morsch
19 Don't delete the TParticle.
20
21 Revision 1.19  2002/02/08 16:50:50  morsch
22 Add name and title in constructor.
23
24 Revision 1.18  2001/11/12 14:31:00  morsch
25 Memory leaks fixed. (M. Bondila)
26
27 Revision 1.17  2001/11/09 09:12:58  morsch
28 Generalization by using AliGenReader object to read particles from file.
29
30 Revision 1.16  2001/07/27 17:09:35  morsch
31 Use local SetTrack, KeepTrack and SetHighWaterMark methods
32 to delegate either to local stack or to stack owned by AliRun.
33 (Piotr Skowronski, A.M.)
34
35 Revision 1.15  2001/01/23 13:29:37  morsch
36 Add method SetParticleCode and enum type Code_t to handle both PDG (new ntuples)
37 and GEANT3 codes (old ntuples) in input file.
38
39 Revision 1.14  2000/12/21 16:24:06  morsch
40 Coding convention clean-up
41
42 Revision 1.13  2000/11/30 07:12:50  alibrary
43 Introducing new Rndm and QA classes
44
45 Revision 1.12  2000/10/27 13:54:45  morsch
46 Remove explicite reference to input file from constuctor.
47
48 Revision 1.11  2000/10/02 21:28:06  fca
49 Removal of useless dependecies via forward declarations
50
51 Revision 1.10  2000/07/11 18:24:55  fca
52 Coding convention corrections + few minor bug fixes
53
54 Revision 1.9  2000/06/14 15:20:09  morsch
55 Include clean-up (IH)
56
57 Revision 1.8  2000/06/09 20:36:44  morsch
58 All coding rule violations except RS3 corrected
59
60 Revision 1.7  2000/02/16 14:56:27  morsch
61 Convert geant particle code into pdg code before putting particle on the stack.
62
63 Revision 1.6  1999/11/09 07:38:48  fca
64 Changes for compatibility with version 2.23 of ROOT
65
66 Revision 1.5  1999/09/29 09:24:12  fca
67 Introduction of the Copyright and cvs Log
68
69 */
70
71
72 // Event generator that using an instance of type AliGenReader
73 // reads particles from a file and applies cuts. 
74
75 #include <iostream.h>
76
77 #include "AliGenExtFile.h"
78 #include "AliRun.h"
79
80 #include <TParticle.h>
81 #include <TFile.h>
82 #include <TTree.h>
83
84
85  ClassImp(AliGenExtFile)
86
87 AliGenExtFile::AliGenExtFile()
88   :AliGenMC()
89 {
90 //  Constructor
91 //
92 //  Read all particles
93     fNpart  =- 1;
94     fReader =  0;
95 }
96
97 AliGenExtFile::AliGenExtFile(Int_t npart)
98     :AliGenMC(npart)
99 {
100 //  Constructor
101     fName   = "ExtFile";
102     fTitle  = "Primaries from ext. File";
103     fReader = 0;
104 }
105
106 AliGenExtFile::AliGenExtFile(const AliGenExtFile & ExtFile)
107 {
108 // copy constructor
109 }
110 //____________________________________________________________
111 AliGenExtFile::~AliGenExtFile()
112 {
113 // Destructor
114     delete fReader;
115 }
116
117 //___________________________________________________________
118 void AliGenExtFile::Init()
119 {
120 // Initialize
121     if (fReader) fReader->Init();
122 }
123
124     
125 void AliGenExtFile::Generate()
126 {
127 // Generate particles
128
129   Float_t polar[3]  = {0,0,0};
130   //
131   Float_t origin[3] = {0,0,0};
132   Float_t p[3];
133   Float_t random[6];
134   Int_t i, j, nt;
135   //
136   for (j=0;j<3;j++) origin[j]=fOrigin[j];
137   if(fVertexSmear == kPerTrack) {
138     Rndm(random,6);
139     for (j = 0; j < 3; j++) {
140         origin[j] += fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
141             TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
142     }
143   }
144
145   while(1) {
146     Int_t nTracks = fReader->NextEvent();       
147     if (nTracks == 0) {
148       // printf("\n No more events !!! !\n");
149       Warning("AliGenExtFile::Generate","\nNo more events in external file!!!\nLast event may be empty or incomplete.\n");
150       return;
151     }
152
153     //
154     // Particle selection loop
155     //
156     // The selction criterium for the external file generator is as follows:
157     //
158     // 1) All tracs are subjects to the cuts defined by AliGenerator, i.e.
159     //    fThetaMin, fThetaMax, fPhiMin, fPhiMax, fPMin, fPMax, fPtMin, fPtMax,
160     //    fYMin, fYMax.
161     //    If the particle does not satisfy these cuts, it is not put on the
162     //    stack.
163     // 2) If fCutOnChild and some specific child is selected (e.g. if
164     //    fForceDecay==kSemiElectronic) the event is rejected if NOT EVEN ONE
165     //    child falls into the child-cuts.
166     if(fCutOnChild) {
167       // Count the selected children
168       Int_t nSelected = 0;
169
170       for (i = 0; i < nTracks; i++) {
171         TParticle* iparticle = fReader->NextParticle();
172         Int_t kf = CheckPDGCode(iparticle->GetPdgCode());
173         kf = TMath::Abs(kf);
174         if (ChildSelected(kf) && KinematicSelection(iparticle, 1)) {
175           nSelected++;
176         }
177       }
178       if (!nSelected) continue;    // No particle selected:  Go to next event
179       fReader->RewindEvent();
180     }
181
182     //
183     // Stack filling loop
184     //
185     for (i = 0; i < nTracks; i++) {
186
187       TParticle* iparticle = fReader->NextParticle();
188       if (!KinematicSelection(iparticle,0)) {
189         Double_t  pz   = iparticle->Pz();
190         Double_t  e    = iparticle->Energy();
191         Double_t  y;
192         if ((e-pz) == 0) {
193           y = 20.;
194         } else if ((e+pz) == 0.) {
195           y = -20.;
196         } else {
197           y = 0.5*TMath::Log((e+pz)/(e-pz));      
198         }
199         printf("\n Not selected %d %f %f %f %f %f", i,
200                iparticle->Theta(),
201                iparticle->Phi(),
202                iparticle->P(),
203                iparticle->Pt(),
204                y);
205         delete iparticle;
206         continue;
207       }
208       p[0] = iparticle->Px();
209       p[1] = iparticle->Py();
210       p[2] = iparticle->Pz();
211       Int_t idpart = iparticle->GetPdgCode();
212       if(fVertexSmear==kPerTrack) {
213           Rndm(random,6);
214           for (j = 0; j < 3; j++) {
215               origin[j]=fOrigin[j]
216                   +fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
217                   TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
218           }
219       }
220       Int_t decayed = iparticle->GetFirstDaughter();
221       Int_t doTracking = fTrackIt && (decayed < 0) &&
222                          (TMath::Abs(idpart) > 10);
223       // printf("*** pdg, first daughter, trk = %d, %d, %d\n",
224       //   idpart,decayed, doTracking);
225       SetTrack(doTracking,-1,idpart,p,origin,polar,0,kPPrimary,nt);
226       KeepTrack(nt);
227     } // track loop
228
229     break;
230
231   } // event loop
232
233   SetHighWaterMark(nt);
234
235   TFile *pFile=0;
236 // Get AliRun object or create it 
237   if (!gAlice) {
238       gAlice = (AliRun*)pFile->Get("gAlice");
239       if (gAlice) printf("AliRun object found on file\n");
240       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
241   }
242   TTree *fAli=gAlice->TreeK();
243   if (fAli) pFile =fAli->GetCurrentFile();
244   pFile->cd();
245 }
246
247
248 AliGenExtFile& AliGenExtFile::operator=(const  AliGenExtFile& rhs)
249 {
250 // Assignment operator
251     return *this;
252 }
253
254
255
256
257
258
259
260