PreReading of MC information on demand.
[u/mrichter/AliRoot.git] / STEER / AliMCEvent.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, 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 //-------------------------------------------------------------------------
19 //     Class for Kinematic Events
20 //     Author: Andreas Morsch, CERN
21 //-------------------------------------------------------------------------
22 #include <TArrow.h>
23 #include <TMarker.h>
24 #include <TH2F.h>
25 #include <TTree.h>
26 #include <TFile.h>
27 #include <TParticle.h>
28 #include <TClonesArray.h>
29 #include <TRefArray.h>
30 #include <TList.h>
31
32 #include "AliLog.h"
33 #include "AliMCEvent.h"
34 #include "AliStack.h"
35 #include "AliTrackReference.h"
36 #include "AliHeader.h"
37 #include "AliGenEventHeader.h"
38
39
40 Int_t AliMCEvent::fgkBgLabelOffset(10000000);
41
42
43 AliMCEvent::AliMCEvent():
44     AliVEvent(),
45     fStack(0),
46     fMCParticles(0),
47     fMCParticleMap(0),
48     fHeader(new AliHeader()),
49     fTRBuffer(0),
50     fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
51     fTreeTR(0),
52     fTmpTreeTR(0),
53     fTmpFileTR(0),
54     fNprimaries(-1),
55     fNparticles(-1),
56     fSubsidiaryEvents(0),
57     fPrimaryOffset(0),
58     fSecondaryOffset(0),
59     fExternal(0)
60 {
61     // Default constructor
62 }
63
64 AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
65     AliVEvent(mcEvnt),
66     fStack(mcEvnt.fStack),
67     fMCParticles(mcEvnt.fMCParticles),
68     fMCParticleMap(mcEvnt.fMCParticleMap),
69     fHeader(mcEvnt.fHeader),
70     fTRBuffer(mcEvnt.fTRBuffer),
71     fTrackReferences(mcEvnt.fTrackReferences),
72     fTreeTR(mcEvnt.fTreeTR),
73     fTmpTreeTR(mcEvnt.fTmpTreeTR),
74     fTmpFileTR(mcEvnt.fTmpFileTR),
75     fNprimaries(mcEvnt.fNprimaries),
76     fNparticles(mcEvnt.fNparticles),
77     fSubsidiaryEvents(0),
78     fPrimaryOffset(0),
79     fSecondaryOffset(0),
80     fExternal(0)
81
82 // Copy constructor
83 }
84
85
86 AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
87 {
88     // assignment operator
89     if (this!=&mcEvnt) { 
90         AliVEvent::operator=(mcEvnt); 
91     }
92   
93     return *this; 
94 }
95
96 void AliMCEvent::ConnectTreeE (TTree* tree)
97 {
98     // Connect the event header tree
99     tree->SetBranchAddress("Header", &fHeader);
100 }
101
102 void AliMCEvent::ConnectTreeK (TTree* tree)
103 {
104     // Connect the kinematics tree to the stack
105     if (!fMCParticles) fMCParticles = new TClonesArray("AliMCParticle",1000);
106     //
107     fStack = fHeader->Stack();
108     fStack->ConnectTree(tree);
109     //
110     // Load the event
111     fStack->GetEvent();
112     fNparticles = fStack->GetNtrack();
113     fNprimaries = fStack->GetNprimary();
114     Int_t iev  = fHeader->GetEvent();
115     Int_t ievr = fHeader->GetEventNrInRun();
116     
117     AliInfo(Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n", 
118                  iev, ievr, fNparticles, fNprimaries));
119  
120     // This is a cache for the TParticles converted to MCParticles on user request
121     if (fMCParticleMap) {
122         fMCParticleMap->Clear();
123         fMCParticles->Delete();
124         if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
125     }
126     else
127         fMCParticleMap = new TRefArray(fNparticles);
128 }
129
130 void AliMCEvent::ConnectTreeTR (TTree* tree)
131 {
132     // Connect the track reference tree
133     fTreeTR = tree;
134     
135     if (fTreeTR->GetBranch("AliRun")) {
136         if (fTmpFileTR) {
137             fTmpFileTR->Close();
138             delete fTmpFileTR;
139         }
140         // This is an old format with one branch per detector not in synch with TreeK
141         ReorderAndExpandTreeTR();
142     } else {
143         // New format 
144         fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
145     }
146 }
147
148 Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
149 {
150     // Retrieve entry i
151     if (i < 0 || i >= fNparticles) {
152         AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
153         particle = 0;
154         trefs    = 0;
155         return (-1);
156     }
157     particle = fStack->Particle(i);
158     if (fTreeTR) {
159         fTreeTR->GetEntry(fStack->TreeKEntry(i));
160         trefs    = fTRBuffer;
161         return trefs->GetEntries();
162     } else {
163         trefs = 0;
164         return -1;
165     }
166 }
167
168
169 void AliMCEvent::Clean()
170 {
171     // Clean-up before new trees are connected
172 //    if (fHeader) {
173 //      delete fHeader;
174 //      fHeader = 0;
175 //    }
176
177     delete fStack; fStack = 0;
178
179     // Clear TR
180     if (fTRBuffer) {
181         fTRBuffer->Delete();
182         delete fTRBuffer;
183         fTRBuffer = 0;
184     }
185 }
186
187 void AliMCEvent::FinishEvent()
188 {
189   // Clean-up after event
190   //    
191     fStack->Reset(0);
192     fMCParticles->Delete();
193     fMCParticleMap->Clear();
194     if (fTRBuffer)
195       fTRBuffer->Delete();
196     fTrackReferences->Delete();
197     fNparticles = -1;
198     fNprimaries = -1;    
199     fStack      =  0;
200 //    fSubsidiaryEvents->Clear();
201     fSubsidiaryEvents = 0;
202 }
203
204
205
206 void AliMCEvent::DrawCheck(Int_t i, Int_t search)
207 {
208     //
209     // Simple event display for debugging
210     if (!fTreeTR) {
211         AliWarning("No Track Reference information available");
212         return;
213     } 
214     
215     if (i > -1 && i < fNparticles) {
216         fTreeTR->GetEntry(fStack->TreeKEntry(i));
217     } else {
218         AliWarning("AliMCEvent::GetEntry: Index out of range");
219     }
220     
221     Int_t nh = fTRBuffer->GetEntries();
222     
223     
224     if (search) {
225         while(nh <= search && i < fNparticles - 1) {
226             i++;
227             fTreeTR->GetEntry(fStack->TreeKEntry(i));
228             nh =  fTRBuffer->GetEntries();
229         }
230         printf("Found Hits at %5d\n", i);
231     }
232     TParticle* particle = fStack->Particle(i);
233     
234     TH2F*    h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
235     Float_t x0 = particle->Vx();
236     Float_t y0 = particle->Vy();
237
238     Float_t x1 = particle->Vx() + particle->Px() * 50.;
239     Float_t y1 = particle->Vy() + particle->Py() * 50.;
240     
241     TArrow*  a = new TArrow(x0, y0, x1, y1, 0.01);
242     h->Draw();
243     a->SetLineColor(2);
244     
245     a->Draw();
246     
247     for (Int_t ih = 0; ih < nh; ih++) {
248         AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
249         TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
250         m->Draw();
251         m->SetMarkerSize(0.4);
252         
253     }
254 }
255
256
257 void AliMCEvent::ReorderAndExpandTreeTR()
258 {
259 //
260 //  Reorder and expand the track reference tree in order to match the kinematics tree.
261 //  Copy the information from different branches into one
262 //
263 //  TreeTR
264
265     fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
266     fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
267     if (!fTRBuffer)  fTRBuffer = new TClonesArray("AliTrackReference", 100);
268     fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 32000, 0);
269     
270
271 //
272 //  Activate the used branches only. Otherwisw we get a bad memory leak.
273     fTreeTR->SetBranchStatus("*",        0);
274     fTreeTR->SetBranchStatus("AliRun.*", 1);
275     fTreeTR->SetBranchStatus("ITS.*",    1);
276     fTreeTR->SetBranchStatus("TPC.*",    1);
277     fTreeTR->SetBranchStatus("TRD.*",    1);
278     fTreeTR->SetBranchStatus("TOF.*",    1);
279     fTreeTR->SetBranchStatus("FRAME.*",  1);
280     fTreeTR->SetBranchStatus("MUON.*",   1);
281 //
282 //  Connect the active branches
283     TClonesArray* trefs[7];
284     for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
285     if (fTreeTR){
286         // make branch for central track references
287         if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
288         if (fTreeTR->GetBranch("ITS"))    fTreeTR->SetBranchAddress("ITS",    &trefs[1]);
289         if (fTreeTR->GetBranch("TPC"))    fTreeTR->SetBranchAddress("TPC",    &trefs[2]);
290         if (fTreeTR->GetBranch("TRD"))    fTreeTR->SetBranchAddress("TRD",    &trefs[3]);
291         if (fTreeTR->GetBranch("TOF"))    fTreeTR->SetBranchAddress("TOF",    &trefs[4]);
292         if (fTreeTR->GetBranch("FRAME"))  fTreeTR->SetBranchAddress("FRAME",  &trefs[5]);
293         if (fTreeTR->GetBranch("MUON"))   fTreeTR->SetBranchAddress("MUON",   &trefs[6]);
294     }
295
296     Int_t np = fStack->GetNprimary();
297     Int_t nt = fTreeTR->GetEntries();
298     
299     //
300     // Loop over tracks and find the secondaries with the help of the kine tree
301     Int_t ifills = 0;
302     Int_t it     = 0;
303     Int_t itlast = 0;
304     TParticle* part;
305
306     for (Int_t ip = np - 1; ip > -1; ip--) {
307         part = fStack->Particle(ip);
308 //      printf("Particle %5d %5d %5d %5d %5d %5d \n", 
309 //             ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), 
310 //             part->GetLastDaughter(), part->TestBit(kTransportBit));
311
312         // Determine range of secondaries produced by this primary during transport     
313         Int_t dau1  = part->GetFirstDaughter();
314         if (dau1 < np) continue;  // This particle has no secondaries produced during transport
315         Int_t dau2  = -1;
316         if (dau1 > -1) {
317             Int_t inext = ip - 1;
318             while (dau2 < 0) {
319                 if (inext >= 0) {
320                     part = fStack->Particle(inext);
321                     dau2 =  part->GetFirstDaughter();
322                     if (dau2 == -1 || dau2 < np) {
323                         dau2 = -1;
324                     } else {
325                         dau2--;
326                     }
327                 } else {
328                     dau2 = fStack->GetNtrack() - 1;
329                 }
330                 inext--;
331             } // find upper bound
332         }  // dau2 < 0
333         
334
335 //      printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
336 //
337 // Loop over reference hits and find secondary label
338 // First the tricky part: find the entry in treeTR than contains the hits or
339 // make sure that no hits exist.
340 //
341         Bool_t hasHits   = kFALSE;
342         Bool_t isOutside = kFALSE;
343
344         it = itlast;
345         while (!hasHits && !isOutside && it < nt) {
346             fTreeTR->GetEntry(it++);
347             for (Int_t ib = 0; ib < 7; ib++) {
348                 if (!trefs[ib]) continue;
349                 Int_t nh = trefs[ib]->GetEntries();
350                 for (Int_t ih = 0; ih < nh; ih++) {
351                     AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
352                     Int_t label = tr->Label();
353                     if (label >= dau1 && label <= dau2) {
354                         hasHits = kTRUE;
355                         itlast = it - 1;
356                         break;
357                     }
358                     if (label > dau2 || label < ip) {
359                         isOutside = kTRUE;
360                         itlast = it - 1;
361                         break;
362                     }
363                 } // hits
364                 if (hasHits || isOutside) break;
365             } // branches
366         } // entries
367
368         if (!hasHits) {
369             // Write empty entries
370             for (Int_t id = dau1; (id <= dau2); id++) {
371                 fTmpTreeTR->Fill();
372                 ifills++;
373             } 
374         } else {
375             // Collect all hits
376             fTreeTR->GetEntry(itlast);
377             for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
378                 for (Int_t ib = 0; ib < 7; ib++) {
379                     if (!trefs[ib]) continue;
380                     Int_t nh = trefs[ib]->GetEntries();
381                     for (Int_t ih = 0; ih < nh; ih++) {
382                         AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
383                         Int_t label = tr->Label();
384                         // Skip primaries
385                         if (label == ip) continue;
386                         if (label > dau2 || label < dau1) 
387                             printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n", 
388                                    itlast, label, dau1, dau2);
389                         if (label == id) {
390                             // secondary found
391                             tr->SetDetectorId(ib-1);
392                             Int_t nref =  fTRBuffer->GetEntriesFast();
393                             TClonesArray &lref = *fTRBuffer;
394                             new(lref[nref]) AliTrackReference(*tr);
395                         }
396                     } // hits
397                 } // branches
398                 fTmpTreeTR->Fill();
399                 fTRBuffer->Delete();
400                 ifills++;
401             } // daughters
402         } // has hits
403     } // tracks
404
405     //
406     // Now loop again and write the primaries
407     //
408     it = nt - 1;
409     for (Int_t ip = 0; ip < np; ip++) {
410         Int_t labmax = -1;
411         while (labmax < ip && it > -1) {
412             fTreeTR->GetEntry(it--);
413             for (Int_t ib = 0; ib < 7; ib++) {
414                 if (!trefs[ib]) continue;
415                 Int_t nh = trefs[ib]->GetEntries();
416                 // 
417                 // Loop over reference hits and find primary labels
418                 for (Int_t ih = 0; ih < nh; ih++) {
419                     AliTrackReference* tr = (AliTrackReference*)  trefs[ib]->At(ih);
420                     Int_t label = tr->Label();
421                     if (label < np && label > labmax) {
422                         labmax = label;
423                     }
424                     
425                     if (label == ip) {
426                         tr->SetDetectorId(ib-1);
427                         Int_t nref = fTRBuffer->GetEntriesFast();
428                         TClonesArray &lref = *fTRBuffer;
429                         new(lref[nref]) AliTrackReference(*tr);
430                     }
431                 } // hits
432             } // branches
433         } // entries
434         it++;
435         fTmpTreeTR->Fill();
436         fTRBuffer->Delete();
437         ifills++;
438     } // tracks
439     // Check
440
441
442     // Clean-up
443     delete fTreeTR; fTreeTR = 0;
444     
445     for (Int_t ib = 0; ib < 7; ib++) {
446         if (trefs[ib]) {
447             trefs[ib]->Clear();
448             delete trefs[ib];
449             trefs[ib] = 0;
450         }
451     }
452
453     if (ifills != fStack->GetNtrack()) 
454         printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", 
455                ifills, fStack->GetNtrack());
456
457     fTmpTreeTR->Write();
458     fTreeTR = fTmpTreeTR;
459 }
460
461 AliVParticle* AliMCEvent::GetTrack(Int_t i) const
462 {
463     // Get MC Particle i
464     //
465
466     if (fExternal) {
467         return ((AliVParticle*) (fMCParticles->At(i)));
468     }
469     
470     //
471     // Check first if this explicitely accesses the subsidiary event
472     
473     if (i >= BgLabelOffset()) {
474         if (fSubsidiaryEvents) {
475             AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
476             return (bgEvent->GetTrack(i - BgLabelOffset()));
477         } else {
478             return 0;
479         }
480     }
481     
482     //
483     AliMCParticle *mcParticle = 0;
484     TParticle     *particle   = 0;
485     TClonesArray  *trefs      = 0;
486     Int_t          ntref      = 0;
487     TRefArray     *rarray     = 0;
488
489
490
491     // Out of range check
492     if (i < 0 || i >= fNparticles) {
493         AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
494         mcParticle = 0;
495         return (mcParticle);
496     }
497
498     
499     if (fSubsidiaryEvents) {
500         AliMCEvent*   mc;
501         Int_t idx = FindIndexAndEvent(i, mc);
502         return (mc->GetTrack(idx));
503     } 
504
505     //
506     // First check If the MC Particle has been already cached
507     if(!fMCParticleMap->At(i)) {
508         // Get particle from the stack
509         particle   = fStack->Particle(i);
510         // Get track references from Tree TR
511         if (fTreeTR) {
512             fTreeTR->GetEntry(fStack->TreeKEntry(i));
513             trefs     = fTRBuffer;
514             ntref     = trefs->GetEntriesFast();
515             rarray    = new TRefArray(ntref);
516             Int_t nen = fTrackReferences->GetEntriesFast();
517             for (Int_t j = 0; j < ntref; j++) {
518                 // Save the track references in a TClonesArray
519                 AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]);
520                 // Save the pointer in a TRefArray
521                 new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
522                 rarray->AddAt((*fTrackReferences)[nen], j);
523                 nen++;
524             } // loop over track references for entry i
525         } // if TreeTR available
526         Int_t nentries = fMCParticles->GetEntriesFast();
527         new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
528         mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
529         fMCParticleMap->AddAt(mcParticle, i);
530
531         TParticle* part = mcParticle->Particle();
532         Int_t imo  = part->GetFirstMother();
533         Int_t id1  = part->GetFirstDaughter();
534         Int_t id2  = part->GetLastDaughter();
535         if (fPrimaryOffset > 0 || fSecondaryOffset > 0) {
536             // Remapping of the mother and daughter indices
537             if (imo < fNprimaries) {
538                 mcParticle->SetMother(imo + fPrimaryOffset);
539             } else {
540                 mcParticle->SetMother(imo + fSecondaryOffset - fNprimaries);
541             }
542
543             if (id1 < fNprimaries) {
544                 mcParticle->SetFirstDaughter(id1 + fPrimaryOffset);
545                 mcParticle->SetLastDaughter (id2 + fPrimaryOffset);
546             } else {
547                 mcParticle->SetFirstDaughter(id1 + fSecondaryOffset - fNprimaries);
548                 mcParticle->SetLastDaughter (id2 + fSecondaryOffset - fNprimaries);
549             }
550
551             
552             if (i > fNprimaries) {
553                 mcParticle->SetLabel(i + fPrimaryOffset);
554             } else {
555                 mcParticle->SetLabel(i + fSecondaryOffset - fNprimaries);
556             }
557             
558             
559         } else {
560             mcParticle->SetFirstDaughter(id1);
561             mcParticle->SetLastDaughter (id2);
562             mcParticle->SetMother       (imo);
563         }
564
565     } else {
566         mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
567     }
568     
569     
570     return mcParticle;
571 }
572
573     AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}
574
575
576 void AliMCEvent::AddSubsidiaryEvent(AliMCEvent* event) 
577 {
578     // Add a subsidiary event to the list; for example merged background event.
579     if (!fSubsidiaryEvents) {
580         TList* events = new TList();
581         events->Add(new AliMCEvent(*this));
582         fSubsidiaryEvents = events;
583     }
584     
585     fSubsidiaryEvents->Add(event);
586 }
587
588 Int_t AliMCEvent::FindIndexAndEvent(Int_t oldidx, AliMCEvent*& event) const
589 {
590     // Find the index and event in case of composed events like signal + background
591     TIter next(fSubsidiaryEvents);
592     next.Reset();
593      if (oldidx < fNprimaries) {
594         while((event = (AliMCEvent*)next())) {
595             if (oldidx < (event->GetPrimaryOffset() + event->GetNumberOfPrimaries())) break;
596         }
597         return (oldidx - event->GetPrimaryOffset());
598     } else {
599         while((event = (AliMCEvent*)next())) {
600             if (oldidx < (event->GetSecondaryOffset() + (event->GetNumberOfTracks() - event->GetNumberOfPrimaries()))) break;
601         }
602         return (oldidx - event->GetSecondaryOffset() + event->GetNumberOfPrimaries());
603     }
604 }
605
606 Int_t AliMCEvent::BgLabelToIndex(Int_t label)
607 {
608     // Convert a background label to an absolute index
609     if (fSubsidiaryEvents) {
610         AliMCEvent* bgEvent = (AliMCEvent*) (fSubsidiaryEvents->At(1));
611         label -= BgLabelOffset();
612         if (label < bgEvent->GetNumberOfPrimaries()) {
613             label += bgEvent->GetPrimaryOffset();
614         } else {
615             label += (bgEvent->GetSecondaryOffset() - fNprimaries);
616         }
617     }
618     return (label);
619 }
620
621
622 Bool_t AliMCEvent::IsPhysicalPrimary(Int_t i) 
623 {
624 //
625 // Delegate to subevent if necesarry 
626
627     
628     if (!fSubsidiaryEvents) {
629         return fStack->IsPhysicalPrimary(i);
630     } else {
631         AliMCEvent* evt = 0;
632         Int_t idx = FindIndexAndEvent(i, evt);
633         return (evt->IsPhysicalPrimary(idx));
634     }
635 }
636
637 void AliMCEvent::InitEvent()
638 {
639 //
640 // Initialize the subsidiary event structure
641     if (fSubsidiaryEvents) {
642         TIter next(fSubsidiaryEvents);
643         AliMCEvent* evt;
644         fNprimaries = 0;
645         fNparticles = 0;
646         
647         while((evt = (AliMCEvent*)next())) {
648             fNprimaries += evt->GetNumberOfPrimaries(); 
649             fNparticles += evt->GetNumberOfTracks();    
650         }
651         
652         Int_t ioffp = 0;
653         Int_t ioffs = fNprimaries;
654         next.Reset();
655         
656         while((evt = (AliMCEvent*)next())) {
657             evt->SetPrimaryOffset(ioffp);
658             evt->SetSecondaryOffset(ioffs);
659             ioffp += evt->GetNumberOfPrimaries();
660             ioffs += (evt->GetNumberOfTracks() - evt->GetNumberOfPrimaries());      
661         }
662     }
663 }
664 void AliMCEvent::PreReadAll()                              
665 {
666     // Preread the MC information
667     Int_t i;
668     // secondaries
669     for (i = fStack->GetNprimary(); i < fStack->GetNtrack(); i++) 
670     {
671         GetTrack(i);
672     }
673     // primaries
674     for (i = 0; i < fStack->GetNprimary(); i++) 
675     {
676         GetTrack(i);
677     }
678     
679     
680 }
681
682
683 ClassImp(AliMCEvent)