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