Typo in conversion factor corrected. (J. Klein)
[u/mrichter/AliRoot.git] / STEER / AliMC.cxx
CommitLineData
5d12ce38 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
d0d4a6b3 18// This class is extracted from the AliRun class
19// and contains all the MC-related functionality
20// The number of dependencies has to be reduced...
21// Author: F.Carminati
22// Federico.Carminati@cern.ch
23
c75cfb51 24#include <RVersion.h>
5d12ce38 25#include <TBrowser.h>
7ca4655f 26#include <TClonesArray.h>
27#include <TGeoManager.h>
5d12ce38 28#include <TStopwatch.h>
29#include <TSystem.h>
30#include <TVirtualMC.h>
9ac3aec9 31#include <TParticle.h>
ced249e6 32#include <TROOT.h>
3cac0921 33#include <TFile.h>
b60e0f5e 34
21bf7095 35#include "AliLog.h"
5d12ce38 36#include "AliDetector.h"
37#include "AliGenerator.h"
38#include "AliHeader.h"
39#include "AliLego.h"
40#include "AliMC.h"
41#include "AliMCQA.h"
42#include "AliRun.h"
942a9039 43#include "AliHit.h"
5d12ce38 44#include "AliStack.h"
0561efeb 45#include "AliMagF.h"
5d12ce38 46#include "AliTrackReference.h"
39b3f8ba 47#include "AliSimulation.h"
392808e8 48#include "AliGeomManager.h"
ced249e6 49#include "AliCDBManager.h"
50#include "AliCDBStorage.h"
51#include "AliCDBEntry.h"
5d12ce38 52
53
54ClassImp(AliMC)
55
56//_______________________________________________________________________
57AliMC::AliMC() :
58 fGenerator(0),
59 fEventEnergy(0),
60 fSummEnergy(0),
61 fSum2Energy(0),
62 fTrRmax(1.e10),
63 fTrZmax(1.e10),
90e48c0c 64 fRDecayMax(1.e10),
38ca2ad6 65 fRDecayMin(-1.),
90e48c0c 66 fDecayPdg(0),
5d12ce38 67 fImedia(0),
68 fTransParName("\0"),
69 fMCQA(0),
70 fHitLists(0),
3cac0921 71 fTmpTreeTR(0),
72 fTmpFileTR(0),
73 fTrackReferences(0),
74 fTmpTrackReferences(0)
5d12ce38 75
76{
d0d4a6b3 77 //default constructor
0561efeb 78 DecayLimits();
5d12ce38 79}
80
81//_______________________________________________________________________
82AliMC::AliMC(const char *name, const char *title) :
83 TVirtualMCApplication(name, title),
84 fGenerator(0),
85 fEventEnergy(0),
86 fSummEnergy(0),
87 fSum2Energy(0),
88 fTrRmax(1.e10),
89 fTrZmax(1.e10),
90e48c0c 90 fRDecayMax(1.e10),
38ca2ad6 91 fRDecayMin(-1.),
90e48c0c 92 fDecayPdg(0),
5d12ce38 93 fImedia(new TArrayI(1000)),
94 fTransParName("\0"),
95 fMCQA(0),
96 fHitLists(new TList()),
3cac0921 97 fTmpTreeTR(0),
98 fTmpFileTR(0),
99 fTrackReferences(new TClonesArray("AliTrackReference", 100)),
100 fTmpTrackReferences(new TClonesArray("AliTrackReference", 100))
5d12ce38 101{
d0d4a6b3 102 //constructor
5d12ce38 103 // Set transport parameters
104 SetTransPar();
0561efeb 105 DecayLimits();
5d12ce38 106 // Prepare the tracking medium lists
107 for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
5d12ce38 108}
109
110//_______________________________________________________________________
111AliMC::AliMC(const AliMC &mc) :
112 TVirtualMCApplication(mc),
113 fGenerator(0),
114 fEventEnergy(0),
115 fSummEnergy(0),
116 fSum2Energy(0),
117 fTrRmax(1.e10),
118 fTrZmax(1.e10),
90e48c0c 119 fRDecayMax(1.e10),
38ca2ad6 120 fRDecayMin(-1.),
90e48c0c 121 fDecayPdg(0),
5d12ce38 122 fImedia(0),
123 fTransParName("\0"),
124 fMCQA(0),
125 fHitLists(0),
3cac0921 126 fTmpTreeTR(0),
127 fTmpFileTR(0),
128 fTrackReferences(0),
129 fTmpTrackReferences(0)
5d12ce38 130{
131 //
90e48c0c 132 // Copy constructor for AliMC
5d12ce38 133 //
134 mc.Copy(*this);
135}
136
137//_______________________________________________________________________
138AliMC::~AliMC()
139{
d0d4a6b3 140 //destructor
5d12ce38 141 delete fGenerator;
142 delete fImedia;
143 delete fMCQA;
144 delete fHitLists;
145 // Delete track references
146 if (fTrackReferences) {
147 fTrackReferences->Delete();
148 delete fTrackReferences;
149 fTrackReferences = 0;
150 }
151
3cac0921 152if (fTmpTrackReferences) {
153 fTmpTrackReferences->Delete();
154 delete fTmpTrackReferences;
155 fTmpTrackReferences = 0;
156 }
157
5d12ce38 158}
159
160//_______________________________________________________________________
6c4904c2 161void AliMC::Copy(TObject &) const
5d12ce38 162{
d0d4a6b3 163 //dummy Copy function
21bf7095 164 AliFatal("Not implemented!");
5d12ce38 165}
166
167//_______________________________________________________________________
168void AliMC::ConstructGeometry()
169{
170 //
4a9de4af 171 // Either load geometry from file or create it through usual
172 // loop on detectors. In the first case the method
173 // AliModule::CreateMaterials() only builds fIdtmed and is postponed
174 // at InitGeometry().
5d12ce38 175 //
176
ced249e6 177 if(gAlice->IsRootGeometry()){ //load geometry either from CDB or from file
178 if(gAlice->IsGeomFromCDB()){
179 AliInfo("Loading geometry from CDB default storage");
180 AliCDBPath path("GRP","Geometry","Data");
181 AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
182 if(!entry) AliFatal("Unable to load geometry from CDB!");
183 entry->SetOwner(0);
184 gGeoManager = (TGeoManager*) entry->GetObject();
185 if (!gGeoManager) AliFatal("TGeoManager object not found in the specified CDB entry!");
4a9de4af 186 }else{
ced249e6 187 // Load geometry
188 const char *geomfilename = gAlice->GetGeometryFileName();
189 if(gSystem->ExpandPathName(geomfilename)){
190 AliInfo(Form("Loading geometry from file:\n %40s",geomfilename));
191 TGeoManager::Import(geomfilename);
192 }else{
193 AliInfo(Form("Geometry file %40s not found!\n",geomfilename));
194 return;
195 }
4a9de4af 196 }
197 }else{
198 // Create modules, materials, geometry
5d12ce38 199 TStopwatch stw;
200 TIter next(gAlice->Modules());
201 AliModule *detector;
21bf7095 202 AliDebug(1, "Geometry creation:");
5d12ce38 203 while((detector = dynamic_cast<AliModule*>(next()))) {
204 stw.Start();
205 // Initialise detector materials and geometry
206 detector->CreateMaterials();
207 detector->CreateGeometry();
21bf7095 208 AliInfo(Form("%10s R:%.2fs C:%.2fs",
209 detector->GetName(),stw.RealTime(),stw.CpuTime()));
5d12ce38 210 }
4a9de4af 211 }
212
5d12ce38 213}
214
215//_______________________________________________________________________
39b3f8ba 216Bool_t AliMC::MisalignGeometry()
217{
218// Call misalignment code if AliSimulation object was defined.
0717eed2 219
ced249e6 220 if(!gAlice->IsRootGeometry()){
221 //Set alignable volumes for the whole geometry
222 SetAllAlignableVolumes();
223 }
0717eed2 224 // Misalign geometry via AliSimulation instance
39b3f8ba 225 if (!AliSimulation::GetInstance()) return kFALSE;
392808e8 226 AliGeomManager::SetGeometry(gGeoManager);
39b3f8ba 227 return AliSimulation::GetInstance()->MisalignGeometry(gAlice->GetRunLoader());
228}
229
230//_______________________________________________________________________
661663fa 231void AliMC::ConstructOpGeometry()
232{
233 //
234 // Loop all detector modules and call DefineOpticalProperties() method
235 //
236
237 TIter next(gAlice->Modules());
238 AliModule *detector;
239 AliInfo("Optical properties definition");
240 while((detector = dynamic_cast<AliModule*>(next()))) {
241 // Initialise detector optical properties
242 detector->DefineOpticalProperties();
243 }
244}
245
246//_______________________________________________________________________
5d12ce38 247void AliMC::InitGeometry()
248{
249 //
27cf3cdc 250 // Initialize detectors
5d12ce38 251 //
252
4a9de4af 253 AliInfo("Initialisation:");
254 TStopwatch stw;
255 TIter next(gAlice->Modules());
256 AliModule *detector;
257 while((detector = dynamic_cast<AliModule*>(next()))) {
258 stw.Start();
27cf3cdc 259 // Initialise detector geometry
4a9de4af 260 if(gAlice->IsRootGeometry()) detector->CreateMaterials();
261 detector->Init();
4a9de4af 262 AliInfo(Form("%10s R:%.2fs C:%.2fs",
263 detector->GetName(),stw.RealTime(),stw.CpuTime()));
264 }
4787b401 265}
266
267//_______________________________________________________________________
268void AliMC::SetAllAlignableVolumes()
269{
270 //
271 // Add alignable volumes (TGeoPNEntries) looping on all
272 // active modules
273 //
274
275 AliInfo(Form("Setting entries for all alignable volumes of active detectors"));
276 AliModule *detector;
277 TIter next(gAlice->Modules());
278 while((detector = dynamic_cast<AliModule*>(next()))) {
279 detector->AddAlignableVolumes();
280 }
5d12ce38 281}
282
283//_______________________________________________________________________
2057aecb 284void AliMC::GeneratePrimaries()
5d12ce38 285{
286 //
287 // Generate primary particles and fill them in the stack.
288 //
289
290 Generator()->Generate();
291}
292
293//_______________________________________________________________________
294void AliMC::SetGenerator(AliGenerator *generator)
295{
296 //
297 // Load the event generator
298 //
299 if(!fGenerator) fGenerator = generator;
300}
301
302//_______________________________________________________________________
303void AliMC::ResetGenerator(AliGenerator *generator)
304{
305 //
306 // Load the event generator
307 //
9a213abc 308 if(fGenerator) {
309 if(generator) {
21bf7095 310 AliWarning(Form("Replacing generator %s with %s",
9a213abc 311 fGenerator->GetName(),generator->GetName()));
312 }
313 else {
21bf7095 314 AliWarning(Form("Replacing generator %s with NULL",
315 fGenerator->GetName()));
9a213abc 316 }
317 }
5d12ce38 318 fGenerator = generator;
319}
320
321//_______________________________________________________________________
322void AliMC::FinishRun()
323{
324 // Clean generator information
21bf7095 325 AliDebug(1, "fGenerator->FinishRun()");
5d12ce38 326 fGenerator->FinishRun();
327
328 //Output energy summary tables
21bf7095 329 AliDebug(1, "EnergySummary()");
330 ToAliDebug(1, EnergySummary());
5d12ce38 331}
332
333//_______________________________________________________________________
334void AliMC::BeginPrimary()
335{
336 //
337 // Called at the beginning of each primary track
338 //
339
340 // Reset Hits info
341 ResetHits();
342 ResetTrackReferences();
5d12ce38 343}
344
345//_______________________________________________________________________
346void AliMC::PreTrack()
347{
d0d4a6b3 348 // Actions before the track's transport
5d12ce38 349 TObjArray &dets = *gAlice->Modules();
350 AliModule *module;
351
352 for(Int_t i=0; i<=gAlice->GetNdets(); i++)
353 if((module = dynamic_cast<AliModule*>(dets[i])))
354 module->PreTrack();
355
356 fMCQA->PreTrack();
357}
358
359//_______________________________________________________________________
360void AliMC::Stepping()
361{
362 //
363 // Called at every step during transport
364 //
0561efeb 365
38ca2ad6 366 Int_t id = DetFromMate(gMC->CurrentMedium());
5d12ce38 367 if (id < 0) return;
368
0561efeb 369
370 if ( gMC->IsNewTrack() &&
371 gMC->TrackTime() == 0. &&
38ca2ad6 372 fRDecayMin >= 0. &&
0561efeb 373 fRDecayMax > fRDecayMin &&
374 gMC->TrackPid() == fDecayPdg )
375 {
376 FixParticleDecaytime();
377 }
378
379
380
5d12ce38 381 //
382 // --- If lego option, do it and leave
383 if (gAlice->Lego())
384 gAlice->Lego()->StepManager();
385 else {
386 Int_t copy;
387 //Update energy deposition tables
388 AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
ab03c084 389 //
390 // write tracke reference for track which is dissapearing - MI
b60e0f5e 391 if (gMC->IsTrackDisappeared()) {
9ac3aec9 392 if (gMC->Etot()>0.05) AddTrackReference(GetCurrentTrackNumber(),
393 AliTrackReference::kDisappeared);
b60e0f5e 394 }
5d12ce38 395
396 //Call the appropriate stepping routine;
397 AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
398 if(det && det->StepManagerIsEnabled()) {
e569033b 399 if(AliLog::GetGlobalDebugLevel()>0) fMCQA->StepManager(id);
5d12ce38 400 det->StepManager();
401 }
402 }
403}
404
405//_______________________________________________________________________
406void AliMC::EnergySummary()
407{
9ac3aec9 408 //e
5d12ce38 409 // Print summary of deposited energy
410 //
411
412 Int_t ndep=0;
413 Float_t edtot=0;
414 Float_t ed, ed2;
415 Int_t kn, i, left, j, id;
416 const Float_t kzero=0;
417 Int_t ievent=gAlice->GetRunLoader()->GetHeader()->GetEvent()+1;
418 //
419 // Energy loss information
420 if(ievent) {
421 printf("***************** Energy Loss Information per event (GEV) *****************\n");
422 for(kn=1;kn<fEventEnergy.GetSize();kn++) {
423 ed=fSummEnergy[kn];
424 if(ed>0) {
425 fEventEnergy[ndep]=kn;
426 if(ievent>1) {
427 ed=ed/ievent;
428 ed2=fSum2Energy[kn];
429 ed2=ed2/ievent;
430 ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
431 } else
432 ed2=99;
433 fSummEnergy[ndep]=ed;
434 fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
435 edtot+=ed;
436 ndep++;
437 }
438 }
439 for(kn=0;kn<(ndep-1)/3+1;kn++) {
440 left=ndep-kn*3;
441 for(i=0;i<(3<left?3:left);i++) {
442 j=kn*3+i;
443 id=Int_t (fEventEnergy[j]+0.1);
444 printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
445 }
446 printf("\n");
447 }
448 //
449 // Relative energy loss in different detectors
450 printf("******************** Relative Energy Loss per event ********************\n");
451 printf("Total energy loss per event %10.3f GeV\n",edtot);
452 for(kn=0;kn<(ndep-1)/5+1;kn++) {
453 left=ndep-kn*5;
454 for(i=0;i<(5<left?5:left);i++) {
455 j=kn*5+i;
456 id=Int_t (fEventEnergy[j]+0.1);
457 printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
458 }
459 printf("\n");
460 }
461 for(kn=0;kn<75;kn++) printf("*");
462 printf("\n");
463 }
464 //
465 // Reset the TArray's
466 // fEventEnergy.Set(0);
467 // fSummEnergy.Set(0);
468 // fSum2Energy.Set(0);
469}
470
471//_____________________________________________________________________________
472void AliMC::BeginEvent()
473{
474 //
475 // Clean-up previous event
476 // Energy scores
21bf7095 477 AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
478 AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
479 AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
480 AliDebug(1, " BEGINNING EVENT ");
481 AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
482 AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
483 AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
5d12ce38 484
3cac0921 485 AliRunLoader *runloader=gAlice->GetRunLoader();
5d12ce38 486
487 /*******************************/
488 /* Clean after eventual */
489 /* previous event */
490 /*******************************/
491
492
493 //Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors),
494 gAlice->SetEventNrInRun(gAlice->GetEventNrInRun()+1);
495 runloader->SetEventNumber(gAlice->GetEventNrInRun());// sets new files, cleans the previous event stuff, if necessary, etc.,
21bf7095 496 AliDebug(1, Form("EventNr is %d",gAlice->GetEventNrInRun()));
5d12ce38 497
498 fEventEnergy.Reset();
499 // Clean detector information
500
501 if (runloader->Stack())
3cac0921 502 runloader->Stack()->Reset();//clean stack -> tree is unloaded
5d12ce38 503 else
3cac0921 504 runloader->MakeStack();//or make a new one
504b172d 505
9ac3aec9 506
3cac0921 507 if(gAlice->Lego() == 0x0)
508 {
509 AliDebug(1, "fRunLoader->MakeTree(K)");
510 runloader->MakeTree("K");
511 }
512
21bf7095 513 AliDebug(1, "gMC->SetStack(fRunLoader->Stack())");
5d12ce38 514 gMC->SetStack(gAlice->GetRunLoader()->Stack());//Was in InitMC - but was moved here
515 //because we don't have guarantee that
516 //stack pointer is not going to change from event to event
517 //since it bellobgs to header and is obtained via RunLoader
518 //
519 // Reset all Detectors & kinematics & make/reset trees
520 //
521
522 runloader->GetHeader()->Reset(gAlice->GetRunNumber(),gAlice->GetEvNumber(),
3cac0921 523 gAlice->GetEventNrInRun());
5d12ce38 524// fRunLoader->WriteKinematics("OVERWRITE"); is there any reason to rewrite here since MakeTree does so
525
504b172d 526 if(gAlice->Lego())
3cac0921 527 {
528 gAlice->Lego()->BeginEvent();
529 return;
530 }
504b172d 531
5d12ce38 532
21bf7095 533 AliDebug(1, "ResetHits()");
5d12ce38 534 ResetHits();
504b172d 535
21bf7095 536 AliDebug(1, "fRunLoader->MakeTree(H)");
5d12ce38 537 runloader->MakeTree("H");
504b172d 538
5d12ce38 539
5d12ce38 540
3cac0921 541 MakeTmpTrackRefsTree();
5d12ce38 542 //create new branches and SetAdresses
543 TIter next(gAlice->Modules());
544 AliModule *detector;
545 while((detector = (AliModule*)next()))
546 {
3cac0921 547 AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName()));
548 detector->MakeBranch("H");
5d12ce38 549 }
550}
551
552//_______________________________________________________________________
553void AliMC::ResetHits()
554{
555 //
556 // Reset all Detectors hits
557 //
558 TIter next(gAlice->Modules());
559 AliModule *detector;
560 while((detector = dynamic_cast<AliModule*>(next()))) {
561 detector->ResetHits();
562 }
563}
564
565//_______________________________________________________________________
566void AliMC::PostTrack()
567{
d0d4a6b3 568 // Posts tracks for each module
4a9de4af 569 TObjArray &dets = *gAlice->Modules();
570 AliModule *module;
571
572 for(Int_t i=0; i<=gAlice->GetNdets(); i++)
573 if((module = dynamic_cast<AliModule*>(dets[i])))
574 module->PostTrack();
5d12ce38 575}
576
577//_______________________________________________________________________
578void AliMC::FinishPrimary()
579{
580 //
581 // Called at the end of each primary track
582 //
583 AliRunLoader *runloader=gAlice->GetRunLoader();
584 // static Int_t count=0;
585 // const Int_t times=10;
586 // This primary is finished, purify stack
c75cfb51 587#if ROOT_VERSION_CODE > 262152
942a9039 588 if (!(gMC->SecondariesAreOrdered())) {
589 runloader->Stack()->ReorderKine();
590 RemapHits();
591 }
c75cfb51 592#endif
5d12ce38 593 runloader->Stack()->PurifyKine();
942a9039 594 RemapHits();
504b172d 595
5d12ce38 596 TIter next(gAlice->Modules());
597 AliModule *detector;
598 while((detector = dynamic_cast<AliModule*>(next()))) {
599 detector->FinishPrimary();
504b172d 600 AliLoader* loader = detector->GetLoader();
601 if(loader)
602 {
603 TTree* treeH = loader->TreeH();
604 if (treeH) treeH->Fill(); //can be Lego run and treeH can not exist
605 }
5d12ce38 606 }
607
608 // Write out track references if any
3cac0921 609 if (fTmpTreeTR) fTmpTreeTR->Fill();
5d12ce38 610}
611
942a9039 612void AliMC::RemapHits()
613{
614//
615// Remaps the track labels of the hits
616 AliRunLoader *runloader=gAlice->GetRunLoader();
617 AliStack* stack = runloader->Stack();
618 TList* hitLists = GetHitLists();
619 TIter next(hitLists);
620 TCollection *hitList;
621
622 while((hitList = dynamic_cast<TCollection*>(next()))) {
623 TIter nexthit(hitList);
624 AliHit *hit;
625 while((hit = dynamic_cast<AliHit*>(nexthit()))) {
626 hit->SetTrack(stack->TrackLabel(hit->GetTrack()));
627 }
628 }
629
630 //
631 // This for detectors which have a special mapping mechanism
632 // for hits, such as TPC and TRD
633 //
634
635 TObjArray* modules = gAlice->Modules();
636 TIter nextmod(modules);
9ac3aec9 637 AliDetector *detector;
638 while((detector = dynamic_cast<AliDetector*>(nextmod()))) {
942a9039 639 detector->RemapTrackHitIDs(stack->TrackLabelMap());
942a9039 640 }
641 //
642 RemapTrackReferencesIDs(stack->TrackLabelMap());
643}
644
5d12ce38 645//_______________________________________________________________________
646void AliMC::FinishEvent()
647{
648 //
649 // Called at the end of the event.
650 //
651
652 //
38ca2ad6 653
5d12ce38 654 if(gAlice->Lego()) gAlice->Lego()->FinishEvent();
655
656 TIter next(gAlice->Modules());
657 AliModule *detector;
658 while((detector = dynamic_cast<AliModule*>(next()))) {
659 detector->FinishEvent();
660 }
661
662 //Update the energy deposit tables
663 Int_t i;
664 for(i=0;i<fEventEnergy.GetSize();i++)
665 {
666 fSummEnergy[i]+=fEventEnergy[i];
667 fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
668 }
669
670 AliRunLoader *runloader=gAlice->GetRunLoader();
671
672 AliHeader* header = runloader->GetHeader();
673 AliStack* stack = runloader->Stack();
674 if ( (header == 0x0) || (stack == 0x0) )
675 {//check if we got header and stack. If not cry and exit aliroot
21bf7095 676 AliFatal("Can not get the stack or header from LOADER");
5d12ce38 677 return;//never reached
678 }
679 // Update Header information
680 header->SetNprimary(stack->GetNprimary());
681 header->SetNtrack(stack->GetNtrack());
682
5d12ce38 683 // Write out the kinematics
38ca2ad6 684 if (!gAlice->Lego()) stack->FinishEvent();
3cac0921 685
686 // Synchronize the TreeTR with TreeK
687 ReorderAndExpandTreeTR();
5d12ce38 688
689 // Write out the event Header information
690 TTree* treeE = runloader->TreeE();
691 if (treeE)
692 {
693 header->SetStack(stack);
694 treeE->Fill();
695 }
696 else
697 {
21bf7095 698 AliError("Can not get TreeE from RL");
5d12ce38 699 }
700
504b172d 701 if(gAlice->Lego() == 0x0)
702 {
703 runloader->WriteKinematics("OVERWRITE");
704 runloader->WriteTrackRefs("OVERWRITE");
705 runloader->WriteHits("OVERWRITE");
706 }
707
21bf7095 708 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
709 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
710 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
711 AliDebug(1, " FINISHING EVENT ");
712 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
713 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
714 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
5d12ce38 715}
716
717//_______________________________________________________________________
718void AliMC::Field(const Double_t* x, Double_t* b) const
719{
d0d4a6b3 720 // Calculates field "b" at point "x"
5d12ce38 721 gAlice->Field(x,b);
722}
723
724//_______________________________________________________________________
725void AliMC::Init()
726{
d0d4a6b3 727 // MC initialization
5d12ce38 728
729 //=================Create Materials and geometry
730 gMC->Init();
4e889ff4 731 // Set alignable volumes for the whole geometry (with old root)
732#if ROOT_VERSION_CODE < 331527
733 SetAllAlignableVolumes();
734#endif
5d12ce38 735 //Read the cuts for all materials
736 ReadTransPar();
737 //Build the special IMEDIA table
738 MediaTable();
739
740 //Compute cross-sections
741 gMC->BuildPhysics();
742
5d12ce38 743 //Initialise geometry deposition table
744 fEventEnergy.Set(gMC->NofVolumes()+1);
745 fSummEnergy.Set(gMC->NofVolumes()+1);
746 fSum2Energy.Set(gMC->NofVolumes()+1);
747
748 //
749 fMCQA = new AliMCQA(gAlice->GetNdets());
750
751 // Register MC in configuration
752 AliConfig::Instance()->Add(gMC);
753
754}
755
756//_______________________________________________________________________
757void AliMC::MediaTable()
758{
759 //
760 // Built media table to get from the media number to
761 // the detector id
762 //
763
764 Int_t kz, nz, idt, lz, i, k, ind;
765 // Int_t ibeg;
766 TObjArray &dets = *gAlice->Detectors();
767 AliModule *det;
768 Int_t ndets=gAlice->GetNdets();
769 //
770 // For all detectors
771 for (kz=0;kz<ndets;kz++) {
772 // If detector is defined
773 if((det=dynamic_cast<AliModule*>(dets[kz]))) {
774 TArrayI &idtmed = *(det->GetIdtmed());
775 for(nz=0;nz<100;nz++) {
0561efeb 776
5d12ce38 777 // Find max and min material number
778 if((idt=idtmed[nz])) {
779 det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
780 det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
781 }
782 }
783 if(det->LoMedium() > det->HiMedium()) {
784 det->LoMedium() = 0;
785 det->HiMedium() = 0;
786 } else {
787 if(det->HiMedium() > fImedia->GetSize()) {
21bf7095 788 AliError(Form("Increase fImedia from %d to %d",
789 fImedia->GetSize(),det->HiMedium()));
5d12ce38 790 return;
791 }
792 // Tag all materials in rage as belonging to detector kz
793 for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
794 (*fImedia)[lz]=kz;
795 }
796 }
797 }
798 }
799 //
800 // Print summary table
21bf7095 801 AliInfo("Tracking media ranges:");
802 ToAliInfo(
5d12ce38 803 for(i=0;i<(ndets-1)/6+1;i++) {
804 for(k=0;k< (6<ndets-i*6?6:ndets-i*6);k++) {
805 ind=i*6+k;
806 det=dynamic_cast<AliModule*>(dets[ind]);
807 if(det)
808 printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
809 det->HiMedium());
810 else
811 printf(" %6s: %3d -> %3d;","NULL",0,0);
812 }
813 printf("\n");
814 }
21bf7095 815 )
5d12ce38 816}
817
818//_______________________________________________________________________
819void AliMC::ReadTransPar()
820{
821 //
822 // Read filename to set the transport parameters
823 //
824
825
826 const Int_t kncuts=10;
827 const Int_t knflags=11;
828 const Int_t knpars=kncuts+knflags;
829 const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
830 "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
831 "BREM","COMP","DCAY","DRAY","HADR","LOSS",
832 "MULS","PAIR","PHOT","RAYL"};
833 char line[256];
834 char detName[7];
835 char* filtmp;
836 Float_t cut[kncuts];
837 Int_t flag[knflags];
838 Int_t i, itmed, iret, ktmed, kz;
839 FILE *lun;
840 //
841 // See whether the file is there
842 filtmp=gSystem->ExpandPathName(fTransParName.Data());
843 lun=fopen(filtmp,"r");
844 delete [] filtmp;
845 if(!lun) {
21bf7095 846 AliWarning(Form("File %s does not exist!",fTransParName.Data()));
5d12ce38 847 return;
848 }
849 //
5d12ce38 850 while(1) {
851 // Initialise cuts and flags
852 for(i=0;i<kncuts;i++) cut[i]=-99;
853 for(i=0;i<knflags;i++) flag[i]=-99;
854 itmed=0;
855 for(i=0;i<256;i++) line[i]='\0';
856 // Read up to the end of line excluded
857 iret=fscanf(lun,"%[^\n]",line);
858 if(iret<0) {
859 //End of file
860 fclose(lun);
5d12ce38 861 return;
862 }
863 // Read the end of line
864 fscanf(lun,"%*c");
865 if(!iret) continue;
866 if(line[0]=='*') continue;
867 // Read the numbers
868 iret=sscanf(line,"%s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d",
869 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
870 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
871 &flag[8],&flag[9],&flag[10]);
872 if(!iret) continue;
873 if(iret<0) {
874 //reading error
21bf7095 875 AliWarning(Form("Error reading file %s",fTransParName.Data()));
5d12ce38 876 continue;
877 }
878 // Check that the module exist
879 AliModule *mod = gAlice->GetModule(detName);
880 if(mod) {
881 // Get the array of media numbers
882 TArrayI &idtmed = *mod->GetIdtmed();
883 // Check that the tracking medium code is valid
884 if(0<=itmed && itmed < 100) {
885 ktmed=idtmed[itmed];
886 if(!ktmed) {
21bf7095 887 AliWarning(Form("Invalid tracking medium code %d for %s",itmed,mod->GetName()));
5d12ce38 888 continue;
889 }
890 // Set energy thresholds
891 for(kz=0;kz<kncuts;kz++) {
892 if(cut[kz]>=0) {
21bf7095 893 AliDebug(2, Form("%-6s set to %10.3E for tracking medium code %4d for %s",
894 kpars[kz],cut[kz],itmed,mod->GetName()));
5d12ce38 895 gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
896 }
897 }
898 // Set transport mechanisms
899 for(kz=0;kz<knflags;kz++) {
900 if(flag[kz]>=0) {
21bf7095 901 AliDebug(2, Form("%-6s set to %10d for tracking medium code %4d for %s",
902 kpars[kncuts+kz],flag[kz],itmed,mod->GetName()));
5d12ce38 903 gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
904 }
905 }
906 } else {
21bf7095 907 AliWarning(Form("Invalid medium code %d",itmed));
5d12ce38 908 continue;
909 }
910 } else {
21bf7095 911 AliDebug(1, Form("%s not present",detName));
5d12ce38 912 continue;
913 }
914 }
915}
916
917//_______________________________________________________________________
918void AliMC::SetTransPar(const char *filename)
919{
920 //
921 // Sets the file name for transport parameters
922 //
923 fTransParName = filename;
924}
925
926//_______________________________________________________________________
0054628d 927void AliMC::Browse(TBrowser *b)
5d12ce38 928{
929 //
930 // Called when the item "Run" is clicked on the left pane
931 // of the Root browser.
932 // It displays the Root Trees and all detectors.
933 //
934 //detectors are in folders anyway
935 b->Add(fMCQA,"AliMCQA");
936}
937
5d12ce38 938//_______________________________________________________________________
939void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
940{
941 //
942 // Add a hit to detector id
943 //
944 TObjArray &dets = *gAlice->Modules();
945 if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
946}
947
948//_______________________________________________________________________
949void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
950{
951 //
952 // Add digit to detector id
953 //
954 TObjArray &dets = *gAlice->Modules();
955 if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
956}
957
958//_______________________________________________________________________
959Int_t AliMC::GetCurrentTrackNumber() const {
960 //
961 // Returns current track
962 //
963 return gAlice->GetRunLoader()->Stack()->GetCurrentTrackNumber();
964}
965
966//_______________________________________________________________________
967void AliMC::DumpPart (Int_t i) const
968{
969 //
970 // Dumps particle i in the stack
971 //
972 AliRunLoader * runloader = gAlice->GetRunLoader();
973 if (runloader->Stack())
974 runloader->Stack()->DumpPart(i);
975}
976
977//_______________________________________________________________________
978void AliMC::DumpPStack () const
979{
980 //
981 // Dumps the particle stack
982 //
983 AliRunLoader * runloader = gAlice->GetRunLoader();
984 if (runloader->Stack())
985 runloader->Stack()->DumpPStack();
986}
987
988//_______________________________________________________________________
989Int_t AliMC::GetNtrack() const {
990 //
991 // Returns number of tracks in stack
992 //
993 Int_t ntracks = -1;
994 AliRunLoader * runloader = gAlice->GetRunLoader();
995 if (runloader->Stack())
996 ntracks = runloader->Stack()->GetNtrack();
997 return ntracks;
998}
999
1000//_______________________________________________________________________
1001Int_t AliMC::GetPrimary(Int_t track) const
1002{
1003 //
1004 // return number of primary that has generated track
1005 //
1006 Int_t nprimary = -999;
1007 AliRunLoader * runloader = gAlice->GetRunLoader();
1008 if (runloader->Stack())
1009 nprimary = runloader->Stack()->GetPrimary(track);
1010 return nprimary;
1011}
1012
1013//_______________________________________________________________________
1014TParticle* AliMC::Particle(Int_t i) const
1015{
d0d4a6b3 1016 // Returns the i-th particle from the stack taking into account
1017 // the remaping done by PurifyKine
5d12ce38 1018 AliRunLoader * runloader = gAlice->GetRunLoader();
1019 if (runloader)
1020 if (runloader->Stack())
1021 return runloader->Stack()->Particle(i);
1022 return 0x0;
1023}
1024
1025//_______________________________________________________________________
1026TObjArray* AliMC::Particles() const {
1027 //
1028 // Returns pointer to Particles array
1029 //
1030 AliRunLoader * runloader = gAlice->GetRunLoader();
1031 if (runloader)
1032 if (runloader->Stack())
1033 return runloader->Stack()->Particles();
1034 return 0x0;
1035}
1036
1037//_______________________________________________________________________
1038void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
1039 Float_t *vpos, Float_t *polar, Float_t tof,
2057aecb 1040 TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
5d12ce38 1041{
1042// Delegate to stack
1043//
1044 AliRunLoader * runloader = gAlice->GetRunLoader();
1045 if (runloader)
1046 if (runloader->Stack())
1047 runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
1048 mech, ntr, weight, is);
1049}
1050
1051//_______________________________________________________________________
1052void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg,
1053 Double_t px, Double_t py, Double_t pz, Double_t e,
1054 Double_t vx, Double_t vy, Double_t vz, Double_t tof,
1055 Double_t polx, Double_t poly, Double_t polz,
2057aecb 1056 TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
5d12ce38 1057{
1058 // Delegate to stack
1059 //
1060 AliRunLoader * runloader = gAlice->GetRunLoader();
1061 if (runloader)
1062 if (runloader->Stack())
1063 runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
1064 polx, poly, polz, mech, ntr, weight, is);
1065}
1066
1067//_______________________________________________________________________
2057aecb 1068void AliMC::SetHighWaterMark(Int_t nt) const
5d12ce38 1069{
1070 //
1071 // Set high water mark for last track in event
1072 AliRunLoader * runloader = gAlice->GetRunLoader();
1073 if (runloader)
1074 if (runloader->Stack())
1075 runloader->Stack()->SetHighWaterMark(nt);
1076}
1077
1078//_______________________________________________________________________
2057aecb 1079void AliMC::KeepTrack(Int_t track) const
5d12ce38 1080{
1081 //
1082 // Delegate to stack
1083 //
1084 AliRunLoader * runloader = gAlice->GetRunLoader();
1085 if (runloader)
1086 if (runloader->Stack())
1087 runloader->Stack()->KeepTrack(track);
1088}
1089
1090//_______________________________________________________________________
2057aecb 1091void AliMC::FlagTrack(Int_t track) const
5d12ce38 1092{
1093 // Delegate to stack
1094 //
1095 AliRunLoader * runloader = gAlice->GetRunLoader();
1096 if (runloader)
1097 if (runloader->Stack())
1098 runloader->Stack()->FlagTrack(track);
1099}
1100
1101//_______________________________________________________________________
2057aecb 1102void AliMC::SetCurrentTrack(Int_t track) const
5d12ce38 1103{
1104 //
1105 // Set current track number
1106 //
1107 AliRunLoader * runloader = gAlice->GetRunLoader();
1108 if (runloader)
1109 if (runloader->Stack())
1110 runloader->Stack()->SetCurrentTrack(track);
1111}
1112
1113//_______________________________________________________________________
9ac3aec9 1114AliTrackReference* AliMC::AddTrackReference(Int_t label, Int_t id)
1115{
5d12ce38 1116 //
1117 // add a trackrefernce to the list
1118 if (!fTrackReferences) {
21bf7095 1119 AliError("Container trackrefernce not active");
9ac3aec9 1120 return NULL;
5d12ce38 1121 }
9ac3aec9 1122 Int_t primary = GetPrimary(label);
1123 Particle(primary)->SetBit(kKeepBit);
1124
3cac0921 1125 Int_t nref = fTmpTrackReferences->GetEntriesFast();
1126 TClonesArray &lref = *fTmpTrackReferences;
9ac3aec9 1127 return new(lref[nref]) AliTrackReference(label, id);
5d12ce38 1128}
1129
1130
1131
1132//_______________________________________________________________________
1133void AliMC::ResetTrackReferences()
1134{
1135 //
1136 // Reset all references
1137 //
3cac0921 1138 if (fTmpTrackReferences) fTmpTrackReferences->Clear();
5d12ce38 1139}
1140
1141void AliMC::RemapTrackReferencesIDs(Int_t *map)
1142{
1143 //
1144 // Remapping track reference
1145 // Called at finish primary
1146 //
3cac0921 1147 if (!fTmpTrackReferences) return;
1148 Int_t nEntries = fTmpTrackReferences->GetEntries();
4e490a96 1149 for (Int_t i=0; i < nEntries; i++){
3cac0921 1150 AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTmpTrackReferences->UncheckedAt(i));
4e490a96 1151 if (ref) {
1152 Int_t newID = map[ref->GetTrack()];
1153 if (newID>=0) ref->SetTrack(newID);
1154 else {
1155 ref->SetBit(kNotDeleted,kFALSE);
3cac0921 1156 fTmpTrackReferences->RemoveAt(i);
4e490a96 1157 }
1158 } // if ref
5d12ce38 1159 }
3cac0921 1160 fTmpTrackReferences->Compress();
5d12ce38 1161}
0561efeb 1162
1163void AliMC::FixParticleDecaytime()
1164{
1165 //
1166 // Fix the particle decay time according to rmin and rmax for decays
1167 //
1168
1169 TLorentzVector p;
1170 gMC->TrackMomentum(p);
1171 Double_t tmin, tmax;
1172 Double_t b;
1173
1174 // Transverse velocity
1175 Double_t vt = p.Pt() / p.E();
1176
1177 if ((b = gAlice->Field()->SolenoidField()) > 0.) { // [kG]
1178
1179 // Radius of helix
1180
1181 Double_t rho = p.Pt() / 0.0003 / b; // [cm]
1182
1183 // Revolution frequency
1184
1185 Double_t omega = vt / rho;
1186
1187 // Maximum and minimum decay time
1188 //
1189 // Check for curlers first
1190 if (fRDecayMax * fRDecayMax / rho / rho / 2. > 1.) return;
1191
1192 //
1193
1194 tmax = TMath::ACos(1. - fRDecayMax * fRDecayMax / rho / rho / 2.) / omega; // [ct]
1195 tmin = TMath::ACos(1. - fRDecayMin * fRDecayMin / rho / rho / 2.) / omega; // [ct]
1196 } else {
1197 tmax = fRDecayMax / vt; // [ct]
1198 tmin = fRDecayMin / vt; // [ct]
1199 }
1200
1201 //
1202 // Dial t using the two limits
1203 Double_t t = tmin + (tmax - tmin) * gRandom->Rndm(); // [ct]
1204 //
1205 //
1206 // Force decay time in transport code
1207 //
fa7ca4da 1208 gMC->ForceDecayTime(t / 2.99792458e10);
0561efeb 1209}
3cac0921 1210
1211void AliMC::MakeTmpTrackRefsTree()
1212{
1213 // Make the temporary track reference tree
1214 fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
1215 fTmpTreeTR = new TTree("TreeTR", "Track References");
1216 if (!fTmpTrackReferences) fTmpTrackReferences = new TClonesArray("AliTrackReference", 100);
1217 fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTmpTrackReferences, 4000);
1218}
1219
1220void AliMC::ReorderAndExpandTreeTR()
1221{
1222//
1223// Reorder and expand the temporary track reference tree in order to match the kinematics tree
1224//
1225
1226 AliRunLoader *rl = gAlice->GetRunLoader();
1227//
1228// TreeTR
1229 AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
1230 rl->MakeTrackRefsContainer();
1231 TTree * treeTR = rl->TreeTR();
1232 if (treeTR){
1233 // make branch for central track references
1234 if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference",0);
1235 TBranch *branch;
1236 branch = treeTR->Branch("TrackReferences",&fTrackReferences);
1237 branch->SetAddress(&fTrackReferences);
1238 }
1239
1240 AliStack* stack = rl->Stack();
1241 Int_t np = stack->GetNprimary();
1242 Int_t nt = fTmpTreeTR->GetEntries();
3cac0921 1243 //
1244 // Loop over tracks and find the secondaries with the help of the kine tree
5ce21fc4 1245 Int_t ifills = 0;
1246 Int_t it = 0;
1247 for (Int_t ip = np - 1; ip > -1; ip--) {
1248 TParticle *part = stack->Particle(ip);
2c4e7c77 1249// printf("Particle %5d %5d %5d %5d %5d \n", ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), part->GetLastDaughter());
1250
5ce21fc4 1251 // Skip primaries that have not been transported
3cac0921 1252 Int_t dau1 = part->GetFirstDaughter();
1253 Int_t dau2 = -1;
2c4e7c77 1254 if ((dau1 > -1 && dau1 < np) || part->GetStatusCode() > 1) continue;
5ce21fc4 1255 //
1256 fTmpTreeTR->GetEntry(it++);
1257 Int_t nh = fTmpTrackReferences->GetEntries();
1258 // Determine range of secondaries produced by this primary
3cac0921 1259 if (dau1 > -1) {
5ce21fc4 1260 Int_t inext = ip - 1;
3cac0921 1261 while (dau2 < 0) {
1262 if (inext >= 0) {
1263 part = stack->Particle(inext);
5ce21fc4 1264 dau2 = part->GetFirstDaughter();
1265 if (dau2 == -1 || dau2 < np) {
1266 dau2 = -1;
1267 } else {
1268 dau2--;
1269 }
3cac0921 1270 } else {
1271 dau2 = stack->GetNtrack() - 1;
1272 }
1273 inext--;
5ce21fc4 1274 } // find upper bound
1275 } // dau2 < 0
2c4e7c77 1276// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
1277
3cac0921 1278 //
1279 // Loop over reference hits and find secondary label
1280 for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
1281 for (Int_t ih = 0; ih < nh; ih++) {
1282 AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences->At(ih);
1283 Int_t label = tr->Label();
5ce21fc4 1284 // Skip primaries
1285 if (label == ip) continue;
2c4e7c77 1286 if (label > dau2 || label < dau1)
1287 AliWarning(Form("Track Reference Label out of range !: %5d %5d %5d \n", label, dau1, dau2));
3cac0921 1288 if (label == id) {
1289 // secondary found
1290 Int_t nref = fTrackReferences->GetEntriesFast();
1291 TClonesArray &lref = *fTrackReferences;
1292 new(lref[nref]) AliTrackReference(*tr);
1293 }
1294 } // hits
1295 treeTR->Fill();
1296 fTrackReferences->Clear();
5ce21fc4 1297 ifills++;
3cac0921 1298 } // daughters
1299 } // tracks
3cac0921 1300 //
1301 // Now loop again and write the primaries
5ce21fc4 1302 it = nt - 1;
1303 for (Int_t ip = 0; ip < np; ip++) {
1304 TParticle* part = stack->Particle(ip);
2c4e7c77 1305 if ((part->GetFirstDaughter() == -1 && part->GetStatusCode() <= 1) || part->GetFirstDaughter() >= np) {
5ce21fc4 1306 // Skip particles that have not been transported
1307 fTmpTreeTR->GetEntry(it--);
1308 Int_t nh = fTmpTrackReferences->GetEntries();
1309 //
1310 // Loop over reference hits and find primary labels
1311 for (Int_t ih = 0; ih < nh; ih++) {
1312 AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences->At(ih);
1313 Int_t label = tr->Label();
1314 if (label == ip) {
1315 Int_t nref = fTrackReferences->GetEntriesFast();
1316 TClonesArray &lref = *fTrackReferences;
1317 new(lref[nref]) AliTrackReference(*tr);
1318 }
1319 }
3cac0921 1320 }
1321 treeTR->Fill();
1322 fTrackReferences->Clear();
5ce21fc4 1323 ifills++;
3cac0921 1324 } // tracks
5ce21fc4 1325 // Check
1326 if (ifills != stack->GetNtrack())
1327 AliWarning(Form("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, stack->GetNtrack()));
3cac0921 1328//
1329// Clean-up
1330 delete fTmpTreeTR;
1331 fTmpFileTR->Close();
1332 delete fTmpFileTR;
1333 delete fTmpTrackReferences;
1334 fTmpTrackReferences = 0;
1335 gSystem->Exec("rm -rf TrackRefsTmp.root");
1336}
1337