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