Line data Source code
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 : // Event generator that using an instance of type AliGenReader
19 : // reads particles from a file and applies cuts.
20 : // Example: In your Config.C you can include the following lines
21 : // AliGenExtFile *gener = new AliGenExtFile(-1);
22 : // gener->SetMomentumRange(0,999);
23 : // gener->SetPhiRange(-180.,180.);
24 : // gener->SetThetaRange(0,180);
25 : // gener->SetYRange(-999,999);
26 : // AliGenReaderTreeK * reader = new AliGenReaderTreeK();
27 : // reader->SetFileName("myFileWithTreeK.root");
28 : // gener->SetReader(reader);
29 : // gener->Init();
30 :
31 :
32 : #include <Riostream.h>
33 :
34 : #include "AliLog.h"
35 : #include "AliGenExtFile.h"
36 : #include "AliRunLoader.h"
37 : #include "AliHeader.h"
38 : #include "AliStack.h"
39 : #include "AliGenEventHeader.h"
40 : #include "AliGenReader.h"
41 :
42 : #include <TParticle.h>
43 : #include <TFile.h>
44 : #include <TTree.h>
45 :
46 :
47 : using std::cout;
48 : using std::endl;
49 : using std::map;
50 :
51 6 : ClassImp(AliGenExtFile)
52 :
53 : AliGenExtFile::AliGenExtFile()
54 0 : :AliGenMC(),
55 0 : fFileName(0),
56 0 : fReader(0),
57 0 : fStartEvent(0)
58 0 : {
59 : // Constructor
60 : //
61 : // Read all particles
62 0 : fNpart = -1;
63 0 : }
64 :
65 : AliGenExtFile::AliGenExtFile(Int_t npart)
66 0 : :AliGenMC(npart),
67 0 : fFileName(0),
68 0 : fReader(0),
69 0 : fStartEvent(0)
70 0 : {
71 : // Constructor
72 0 : fName = "ExtFile";
73 0 : fTitle = "Primaries from ext. File";
74 0 : }
75 :
76 : //____________________________________________________________
77 0 : AliGenExtFile::~AliGenExtFile()
78 0 : {
79 : // Destructor
80 0 : delete fReader;
81 0 : }
82 :
83 : //___________________________________________________________
84 : void AliGenExtFile::Init()
85 : {
86 : // Initialize
87 0 : if (fReader) fReader->Init();
88 0 : }
89 :
90 : //___________________________________________________________
91 : void AliGenExtFile::Generate()
92 : {
93 : // Generate particles
94 :
95 : Double_t polar[3] = {0,0,0};
96 : //
97 0 : Double_t origin[3] = {0,0,0};
98 : Double_t time = 0.;
99 : Double_t p[4];
100 0 : Float_t random[6];
101 0 : Int_t i=0, j, nt;
102 : //
103 : //
104 : // Fast forward up to start Event
105 0 : for (Int_t ie=0; ie<fStartEvent; ++ie ) {
106 0 : Int_t nTracks = fReader->NextEvent();
107 0 : if (nTracks == 0) {
108 : // printf("\n No more events !!! !\n");
109 0 : Warning("AliGenExtFile::Generate","\nNo more events in external file!!!\nLast event may be empty or incomplete.\n");
110 0 : return;
111 : }
112 0 : for (i = 0; i < nTracks; ++i) {
113 0 : if (NULL == fReader->NextParticle())
114 0 : AliFatal("Error while skipping tracks");
115 : }
116 0 : cout << "Skipping event " << ie << endl;
117 0 : }
118 0 : fStartEvent = 0; // do not skip events the second time
119 :
120 0 : while(1) {
121 0 : if (fVertexSmear == kPerEvent) Vertex();
122 0 : Int_t nTracks = fReader->NextEvent();
123 0 : if (nTracks == 0) {
124 : // printf("\n No more events !!! !\n");
125 0 : Warning("AliGenExtFile::Generate","\nNo more events in external file!!!\nLast event may be empty or incomplete.\n");
126 0 : return;
127 : }
128 :
129 : //
130 : // Particle selection loop
131 : //
132 : // The selection criterium for the external file generator is as follows:
133 : //
134 : // 1) All tracks are subject to the cuts defined by AliGenerator, i.e.
135 : // fThetaMin, fThetaMax, fPhiMin, fPhiMax, fPMin, fPMax, fPtMin, fPtMax,
136 : // fYMin, fYMax.
137 : // If the particle does not satisfy these cuts, it is not put on the
138 : // stack.
139 : // 2) If fCutOnChild and some specific child is selected (e.g. if
140 : // fForceDecay==kSemiElectronic) the event is rejected if NOT EVEN ONE
141 : // child falls into the child-cuts.
142 : TParticle* iparticle = 0x0;
143 :
144 0 : if(fCutOnChild) {
145 : // Count the selected children
146 : Int_t nSelected = 0;
147 0 : while ((iparticle=fReader->NextParticle()) ) {
148 0 : Int_t kf = CheckPDGCode(iparticle->GetPdgCode());
149 0 : kf = TMath::Abs(kf);
150 0 : if (ChildSelected(kf) && KinematicSelection(iparticle, 1)) {
151 0 : nSelected++;
152 0 : }
153 : }
154 0 : if (!nSelected) continue; // No particle selected: Go to next event
155 0 : fReader->RewindEvent();
156 0 : }
157 :
158 : //
159 : // Stack selection loop
160 : //
161 0 : class SelectorLogic { // need to do recursive back tracking, requires a "nested" function
162 : private:
163 : Int_t idCount;
164 : map<Int_t,Int_t> firstMotherMap;
165 : map<Int_t,Int_t> secondMotherMap;
166 : map<Int_t,Bool_t> selectedIdMap;
167 : map<Int_t,Int_t> newIdMap;
168 : void selectMothersToo(Int_t particleId) {
169 0 : Int_t mum1 = firstMotherMap[particleId];
170 0 : if (mum1 > -1 && !selectedIdMap[mum1]) {
171 0 : selectedIdMap[mum1] = true;
172 0 : selectMothersToo(mum1);
173 0 : }
174 0 : Int_t mum2 = secondMotherMap[particleId];
175 0 : if (mum2 > -1 && !selectedIdMap[mum2]) {
176 0 : selectedIdMap[mum2] = true;
177 0 : selectMothersToo(mum2);
178 0 : }
179 0 : }
180 : public:
181 0 : SelectorLogic():idCount(0), firstMotherMap(), secondMotherMap(), selectedIdMap(), newIdMap() {}
182 : void init() {
183 0 : idCount = 0;
184 0 : }
185 : void setData(Int_t id, Int_t mum1, Int_t mum2, Bool_t selected) {
186 0 : idCount++; // we know that this function is called in succession of ids, so counting is fine to determine max id
187 0 : firstMotherMap[id] = mum1;
188 0 : secondMotherMap[id] = mum2;
189 0 : selectedIdMap[id] = selected;
190 0 : }
191 : void reselectCuttedMothersAndRemapIDs() {
192 0 : for (Int_t id = 0; id < idCount; ++id) {
193 0 : if (selectedIdMap[id]) {
194 0 : selectMothersToo(id);
195 0 : }
196 : }
197 : Int_t newId0 = 0;
198 0 : for (Int_t id = 0; id < idCount; id++) {
199 0 : if (selectedIdMap[id]) {
200 0 : newIdMap[id] = newId0; ++newId0;
201 0 : } else {
202 0 : newIdMap[id] = -1;
203 : }
204 : }
205 0 : }
206 : Bool_t isSelected(Int_t id) {
207 0 : return selectedIdMap[id];
208 : }
209 : Int_t newId(Int_t id) {
210 0 : if (id == -1) return -1;
211 0 : return newIdMap[id];
212 0 : }
213 : };
214 0 : SelectorLogic selector;
215 0 : selector.init();
216 0 : for (i = 0; i < nTracks; i++) {
217 0 : TParticle* jparticle = fReader->NextParticle();
218 0 : selector.setData(i,
219 0 : jparticle->GetFirstMother(),
220 0 : jparticle->GetSecondMother(),
221 0 : KinematicSelection(jparticle,0));
222 : }
223 0 : selector.reselectCuttedMothersAndRemapIDs();
224 0 : fReader->RewindEvent();
225 :
226 : //
227 : // Stack filling loop
228 : //
229 0 : fNprimaries = 0;
230 0 : for (i = 0; i < nTracks; i++) {
231 0 : TParticle* jparticle = fReader->NextParticle();
232 0 : Bool_t selected = selector.isSelected(i);
233 0 : if (!selected) {
234 0 : continue;
235 : }
236 0 : Int_t parent = selector.newId(jparticle->GetFirstMother());
237 : // printf("particle %d -> %d, with mother %d -> %d\n", i, selector.newId(i), jparticle->GetFirstMother(), parent);
238 :
239 0 : p[0] = jparticle->Px();
240 0 : p[1] = jparticle->Py();
241 0 : p[2] = jparticle->Pz();
242 0 : p[3] = jparticle->Energy();
243 :
244 0 : Int_t idpart = jparticle->GetPdgCode();
245 0 : if(fVertexSmear==kPerTrack)
246 : {
247 0 : Rndm(random,6);
248 0 : for (j = 0; j < 3; j++) {
249 0 : origin[j]=fOrigin[j]+
250 0 : fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
251 0 : TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
252 : }
253 0 : Rndm(random,2);
254 0 : time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
255 0 : TMath::Cos(2*random[0]*TMath::Pi())*
256 0 : TMath::Sqrt(-2*TMath::Log(random[1]));
257 0 : } else {
258 0 : origin[0] = fVertex[0] + jparticle->Vx();
259 0 : origin[1] = fVertex[1] + jparticle->Vy();
260 0 : origin[2] = fVertex[2] + jparticle->Vz();
261 0 : time = fTime + jparticle->T();
262 : }
263 0 : Int_t doTracking = fTrackIt && selected && (jparticle->TestBit(kTransportBit));
264 :
265 0 : PushTrack(doTracking, parent, idpart,
266 0 : p[0], p[1], p[2], p[3], origin[0], origin[1], origin[2], time,
267 : polar[0], polar[1], polar[2],
268 0 : kPPrimary, nt, 1., jparticle->GetStatusCode());
269 :
270 0 : KeepTrack(nt);
271 0 : fNprimaries++;
272 0 : } // track loop
273 :
274 : // Generated event header
275 0 : AliGenEventHeader * header = fReader->GetGenEventHeader();
276 0 : if (!header) header = new AliGenEventHeader();
277 0 : header->SetNProduced(fNprimaries);
278 0 : header->SetPrimaryVertex(fVertex);
279 0 : header->SetInteractionTime(fTime);
280 0 : AddHeader(header);
281 : break;
282 :
283 0 : } // event loop
284 :
285 0 : SetHighWaterMark(nt);
286 0 : CdEventFile();
287 0 : }
288 :
289 : //___________________________________________________________
290 : void AliGenExtFile::CdEventFile()
291 : {
292 : // CD back to the event file
293 0 : AliRunLoader::Instance()->CdGAFile();
294 0 : }
295 :
296 :
297 :
298 :
299 :
|
