Use SetInputEvent()
[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();
ac2b85c0 594
942a9039 595 RemapHits();
504b172d 596
5d12ce38 597 TIter next(gAlice->Modules());
598 AliModule *detector;
599 while((detector = dynamic_cast<AliModule*>(next()))) {
600 detector->FinishPrimary();
504b172d 601 AliLoader* loader = detector->GetLoader();
602 if(loader)
603 {
604 TTree* treeH = loader->TreeH();
605 if (treeH) treeH->Fill(); //can be Lego run and treeH can not exist
606 }
5d12ce38 607 }
608
609 // Write out track references if any
3cac0921 610 if (fTmpTreeTR) fTmpTreeTR->Fill();
5d12ce38 611}
612
942a9039 613void AliMC::RemapHits()
614{
615//
616// Remaps the track labels of the hits
617 AliRunLoader *runloader=gAlice->GetRunLoader();
618 AliStack* stack = runloader->Stack();
619 TList* hitLists = GetHitLists();
620 TIter next(hitLists);
621 TCollection *hitList;
622
623 while((hitList = dynamic_cast<TCollection*>(next()))) {
624 TIter nexthit(hitList);
625 AliHit *hit;
626 while((hit = dynamic_cast<AliHit*>(nexthit()))) {
627 hit->SetTrack(stack->TrackLabel(hit->GetTrack()));
628 }
629 }
630
631 //
632 // This for detectors which have a special mapping mechanism
633 // for hits, such as TPC and TRD
634 //
ac2b85c0 635
942a9039 636
637 TObjArray* modules = gAlice->Modules();
638 TIter nextmod(modules);
ac2b85c0 639 AliModule *module;
640 while((module = (AliModule*) nextmod())) {
641 AliDetector* det = dynamic_cast<AliDetector*> (module);
642 if (det) det->RemapTrackHitIDs(stack->TrackLabelMap());
942a9039 643 }
644 //
645 RemapTrackReferencesIDs(stack->TrackLabelMap());
646}
647
5d12ce38 648//_______________________________________________________________________
649void AliMC::FinishEvent()
650{
651 //
652 // Called at the end of the event.
653 //
654
655 //
38ca2ad6 656
5d12ce38 657 if(gAlice->Lego()) gAlice->Lego()->FinishEvent();
658
659 TIter next(gAlice->Modules());
660 AliModule *detector;
661 while((detector = dynamic_cast<AliModule*>(next()))) {
662 detector->FinishEvent();
663 }
664
665 //Update the energy deposit tables
666 Int_t i;
667 for(i=0;i<fEventEnergy.GetSize();i++)
668 {
669 fSummEnergy[i]+=fEventEnergy[i];
670 fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
671 }
672
673 AliRunLoader *runloader=gAlice->GetRunLoader();
674
675 AliHeader* header = runloader->GetHeader();
676 AliStack* stack = runloader->Stack();
677 if ( (header == 0x0) || (stack == 0x0) )
678 {//check if we got header and stack. If not cry and exit aliroot
21bf7095 679 AliFatal("Can not get the stack or header from LOADER");
5d12ce38 680 return;//never reached
681 }
682 // Update Header information
683 header->SetNprimary(stack->GetNprimary());
684 header->SetNtrack(stack->GetNtrack());
685
5d12ce38 686 // Write out the kinematics
38ca2ad6 687 if (!gAlice->Lego()) stack->FinishEvent();
3cac0921 688
689 // Synchronize the TreeTR with TreeK
172bc1fb 690 if (fTmpTreeTR) ReorderAndExpandTreeTR();
5d12ce38 691
692 // Write out the event Header information
693 TTree* treeE = runloader->TreeE();
694 if (treeE)
695 {
696 header->SetStack(stack);
697 treeE->Fill();
698 }
699 else
700 {
21bf7095 701 AliError("Can not get TreeE from RL");
5d12ce38 702 }
703
504b172d 704 if(gAlice->Lego() == 0x0)
705 {
706 runloader->WriteKinematics("OVERWRITE");
707 runloader->WriteTrackRefs("OVERWRITE");
708 runloader->WriteHits("OVERWRITE");
709 }
710
21bf7095 711 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
712 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
713 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
714 AliDebug(1, " FINISHING EVENT ");
715 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
716 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
717 AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
5d12ce38 718}
719
720//_______________________________________________________________________
721void AliMC::Field(const Double_t* x, Double_t* b) const
722{
d0d4a6b3 723 // Calculates field "b" at point "x"
5d12ce38 724 gAlice->Field(x,b);
725}
726
727//_______________________________________________________________________
728void AliMC::Init()
729{
d0d4a6b3 730 // MC initialization
5d12ce38 731
732 //=================Create Materials and geometry
733 gMC->Init();
4e889ff4 734 // Set alignable volumes for the whole geometry (with old root)
735#if ROOT_VERSION_CODE < 331527
736 SetAllAlignableVolumes();
737#endif
5d12ce38 738 //Read the cuts for all materials
739 ReadTransPar();
740 //Build the special IMEDIA table
741 MediaTable();
742
743 //Compute cross-sections
744 gMC->BuildPhysics();
745
5d12ce38 746 //Initialise geometry deposition table
747 fEventEnergy.Set(gMC->NofVolumes()+1);
748 fSummEnergy.Set(gMC->NofVolumes()+1);
749 fSum2Energy.Set(gMC->NofVolumes()+1);
750
751 //
752 fMCQA = new AliMCQA(gAlice->GetNdets());
753
754 // Register MC in configuration
755 AliConfig::Instance()->Add(gMC);
756
757}
758
759//_______________________________________________________________________
760void AliMC::MediaTable()
761{
762 //
763 // Built media table to get from the media number to
764 // the detector id
765 //
766
767 Int_t kz, nz, idt, lz, i, k, ind;
768 // Int_t ibeg;
769 TObjArray &dets = *gAlice->Detectors();
770 AliModule *det;
771 Int_t ndets=gAlice->GetNdets();
772 //
773 // For all detectors
774 for (kz=0;kz<ndets;kz++) {
775 // If detector is defined
776 if((det=dynamic_cast<AliModule*>(dets[kz]))) {
777 TArrayI &idtmed = *(det->GetIdtmed());
778 for(nz=0;nz<100;nz++) {
0561efeb 779
5d12ce38 780 // Find max and min material number
781 if((idt=idtmed[nz])) {
782 det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
783 det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
784 }
785 }
786 if(det->LoMedium() > det->HiMedium()) {
787 det->LoMedium() = 0;
788 det->HiMedium() = 0;
789 } else {
790 if(det->HiMedium() > fImedia->GetSize()) {
21bf7095 791 AliError(Form("Increase fImedia from %d to %d",
792 fImedia->GetSize(),det->HiMedium()));
5d12ce38 793 return;
794 }
795 // Tag all materials in rage as belonging to detector kz
796 for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
797 (*fImedia)[lz]=kz;
798 }
799 }
800 }
801 }
802 //
803 // Print summary table
21bf7095 804 AliInfo("Tracking media ranges:");
805 ToAliInfo(
5d12ce38 806 for(i=0;i<(ndets-1)/6+1;i++) {
807 for(k=0;k< (6<ndets-i*6?6:ndets-i*6);k++) {
808 ind=i*6+k;
809 det=dynamic_cast<AliModule*>(dets[ind]);
810 if(det)
811 printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
812 det->HiMedium());
813 else
814 printf(" %6s: %3d -> %3d;","NULL",0,0);
815 }
816 printf("\n");
817 }
21bf7095 818 )
5d12ce38 819}
820
821//_______________________________________________________________________
822void AliMC::ReadTransPar()
823{
824 //
825 // Read filename to set the transport parameters
826 //
827
828
829 const Int_t kncuts=10;
830 const Int_t knflags=11;
831 const Int_t knpars=kncuts+knflags;
832 const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
833 "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
834 "BREM","COMP","DCAY","DRAY","HADR","LOSS",
835 "MULS","PAIR","PHOT","RAYL"};
836 char line[256];
837 char detName[7];
838 char* filtmp;
839 Float_t cut[kncuts];
840 Int_t flag[knflags];
841 Int_t i, itmed, iret, ktmed, kz;
842 FILE *lun;
843 //
844 // See whether the file is there
845 filtmp=gSystem->ExpandPathName(fTransParName.Data());
846 lun=fopen(filtmp,"r");
847 delete [] filtmp;
848 if(!lun) {
21bf7095 849 AliWarning(Form("File %s does not exist!",fTransParName.Data()));
5d12ce38 850 return;
851 }
852 //
5d12ce38 853 while(1) {
854 // Initialise cuts and flags
855 for(i=0;i<kncuts;i++) cut[i]=-99;
856 for(i=0;i<knflags;i++) flag[i]=-99;
857 itmed=0;
858 for(i=0;i<256;i++) line[i]='\0';
859 // Read up to the end of line excluded
860 iret=fscanf(lun,"%[^\n]",line);
861 if(iret<0) {
862 //End of file
863 fclose(lun);
5d12ce38 864 return;
865 }
866 // Read the end of line
867 fscanf(lun,"%*c");
868 if(!iret) continue;
869 if(line[0]=='*') continue;
870 // Read the numbers
871 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",
872 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
873 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
874 &flag[8],&flag[9],&flag[10]);
875 if(!iret) continue;
876 if(iret<0) {
877 //reading error
21bf7095 878 AliWarning(Form("Error reading file %s",fTransParName.Data()));
5d12ce38 879 continue;
880 }
881 // Check that the module exist
882 AliModule *mod = gAlice->GetModule(detName);
883 if(mod) {
884 // Get the array of media numbers
885 TArrayI &idtmed = *mod->GetIdtmed();
886 // Check that the tracking medium code is valid
887 if(0<=itmed && itmed < 100) {
888 ktmed=idtmed[itmed];
889 if(!ktmed) {
21bf7095 890 AliWarning(Form("Invalid tracking medium code %d for %s",itmed,mod->GetName()));
5d12ce38 891 continue;
892 }
893 // Set energy thresholds
894 for(kz=0;kz<kncuts;kz++) {
895 if(cut[kz]>=0) {
21bf7095 896 AliDebug(2, Form("%-6s set to %10.3E for tracking medium code %4d for %s",
897 kpars[kz],cut[kz],itmed,mod->GetName()));
5d12ce38 898 gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
899 }
900 }
901 // Set transport mechanisms
902 for(kz=0;kz<knflags;kz++) {
903 if(flag[kz]>=0) {
21bf7095 904 AliDebug(2, Form("%-6s set to %10d for tracking medium code %4d for %s",
905 kpars[kncuts+kz],flag[kz],itmed,mod->GetName()));
5d12ce38 906 gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
907 }
908 }
909 } else {
21bf7095 910 AliWarning(Form("Invalid medium code %d",itmed));
5d12ce38 911 continue;
912 }
913 } else {
21bf7095 914 AliDebug(1, Form("%s not present",detName));
5d12ce38 915 continue;
916 }
917 }
918}
919
920//_______________________________________________________________________
921void AliMC::SetTransPar(const char *filename)
922{
923 //
924 // Sets the file name for transport parameters
925 //
926 fTransParName = filename;
927}
928
929//_______________________________________________________________________
0054628d 930void AliMC::Browse(TBrowser *b)
5d12ce38 931{
932 //
933 // Called when the item "Run" is clicked on the left pane
934 // of the Root browser.
935 // It displays the Root Trees and all detectors.
936 //
937 //detectors are in folders anyway
938 b->Add(fMCQA,"AliMCQA");
939}
940
5d12ce38 941//_______________________________________________________________________
942void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
943{
944 //
945 // Add a hit to detector id
946 //
947 TObjArray &dets = *gAlice->Modules();
948 if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
949}
950
951//_______________________________________________________________________
952void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
953{
954 //
955 // Add digit to detector id
956 //
957 TObjArray &dets = *gAlice->Modules();
958 if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
959}
960
961//_______________________________________________________________________
962Int_t AliMC::GetCurrentTrackNumber() const {
963 //
964 // Returns current track
965 //
966 return gAlice->GetRunLoader()->Stack()->GetCurrentTrackNumber();
967}
968
969//_______________________________________________________________________
970void AliMC::DumpPart (Int_t i) const
971{
972 //
973 // Dumps particle i in the stack
974 //
975 AliRunLoader * runloader = gAlice->GetRunLoader();
976 if (runloader->Stack())
977 runloader->Stack()->DumpPart(i);
978}
979
980//_______________________________________________________________________
981void AliMC::DumpPStack () const
982{
983 //
984 // Dumps the particle stack
985 //
986 AliRunLoader * runloader = gAlice->GetRunLoader();
987 if (runloader->Stack())
988 runloader->Stack()->DumpPStack();
989}
990
991//_______________________________________________________________________
992Int_t AliMC::GetNtrack() const {
993 //
994 // Returns number of tracks in stack
995 //
996 Int_t ntracks = -1;
997 AliRunLoader * runloader = gAlice->GetRunLoader();
998 if (runloader->Stack())
999 ntracks = runloader->Stack()->GetNtrack();
1000 return ntracks;
1001}
1002
1003//_______________________________________________________________________
1004Int_t AliMC::GetPrimary(Int_t track) const
1005{
1006 //
1007 // return number of primary that has generated track
1008 //
1009 Int_t nprimary = -999;
1010 AliRunLoader * runloader = gAlice->GetRunLoader();
1011 if (runloader->Stack())
1012 nprimary = runloader->Stack()->GetPrimary(track);
1013 return nprimary;
1014}
1015
1016//_______________________________________________________________________
1017TParticle* AliMC::Particle(Int_t i) const
1018{
d0d4a6b3 1019 // Returns the i-th particle from the stack taking into account
1020 // the remaping done by PurifyKine
5d12ce38 1021 AliRunLoader * runloader = gAlice->GetRunLoader();
1022 if (runloader)
1023 if (runloader->Stack())
1024 return runloader->Stack()->Particle(i);
1025 return 0x0;
1026}
1027
1028//_______________________________________________________________________
1029TObjArray* AliMC::Particles() const {
1030 //
1031 // Returns pointer to Particles array
1032 //
1033 AliRunLoader * runloader = gAlice->GetRunLoader();
1034 if (runloader)
1035 if (runloader->Stack())
1036 return runloader->Stack()->Particles();
1037 return 0x0;
1038}
1039
1040//_______________________________________________________________________
1041void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
1042 Float_t *vpos, Float_t *polar, Float_t tof,
2057aecb 1043 TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
5d12ce38 1044{
1045// Delegate to stack
1046//
1047 AliRunLoader * runloader = gAlice->GetRunLoader();
1048 if (runloader)
1049 if (runloader->Stack())
1050 runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
1051 mech, ntr, weight, is);
1052}
1053
1054//_______________________________________________________________________
1055void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg,
1056 Double_t px, Double_t py, Double_t pz, Double_t e,
1057 Double_t vx, Double_t vy, Double_t vz, Double_t tof,
1058 Double_t polx, Double_t poly, Double_t polz,
2057aecb 1059 TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
5d12ce38 1060{
1061 // Delegate to stack
1062 //
1063 AliRunLoader * runloader = gAlice->GetRunLoader();
1064 if (runloader)
1065 if (runloader->Stack())
1066 runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
1067 polx, poly, polz, mech, ntr, weight, is);
1068}
1069
1070//_______________________________________________________________________
2057aecb 1071void AliMC::SetHighWaterMark(Int_t nt) const
5d12ce38 1072{
1073 //
1074 // Set high water mark for last track in event
1075 AliRunLoader * runloader = gAlice->GetRunLoader();
1076 if (runloader)
1077 if (runloader->Stack())
1078 runloader->Stack()->SetHighWaterMark(nt);
1079}
1080
1081//_______________________________________________________________________
2057aecb 1082void AliMC::KeepTrack(Int_t track) const
5d12ce38 1083{
1084 //
1085 // Delegate to stack
1086 //
1087 AliRunLoader * runloader = gAlice->GetRunLoader();
1088 if (runloader)
1089 if (runloader->Stack())
1090 runloader->Stack()->KeepTrack(track);
1091}
1092
1093//_______________________________________________________________________
2057aecb 1094void AliMC::FlagTrack(Int_t track) const
5d12ce38 1095{
1096 // Delegate to stack
1097 //
1098 AliRunLoader * runloader = gAlice->GetRunLoader();
1099 if (runloader)
1100 if (runloader->Stack())
1101 runloader->Stack()->FlagTrack(track);
1102}
1103
1104//_______________________________________________________________________
2057aecb 1105void AliMC::SetCurrentTrack(Int_t track) const
5d12ce38 1106{
1107 //
1108 // Set current track number
1109 //
1110 AliRunLoader * runloader = gAlice->GetRunLoader();
1111 if (runloader)
1112 if (runloader->Stack())
1113 runloader->Stack()->SetCurrentTrack(track);
1114}
1115
1116//_______________________________________________________________________
9ac3aec9 1117AliTrackReference* AliMC::AddTrackReference(Int_t label, Int_t id)
1118{
5d12ce38 1119 //
1120 // add a trackrefernce to the list
1121 if (!fTrackReferences) {
21bf7095 1122 AliError("Container trackrefernce not active");
9ac3aec9 1123 return NULL;
5d12ce38 1124 }
9ac3aec9 1125 Int_t primary = GetPrimary(label);
1126 Particle(primary)->SetBit(kKeepBit);
1127
3cac0921 1128 Int_t nref = fTmpTrackReferences->GetEntriesFast();
1129 TClonesArray &lref = *fTmpTrackReferences;
9ac3aec9 1130 return new(lref[nref]) AliTrackReference(label, id);
5d12ce38 1131}
1132
1133
1134
1135//_______________________________________________________________________
1136void AliMC::ResetTrackReferences()
1137{
1138 //
1139 // Reset all references
1140 //
3cac0921 1141 if (fTmpTrackReferences) fTmpTrackReferences->Clear();
5d12ce38 1142}
1143
1144void AliMC::RemapTrackReferencesIDs(Int_t *map)
1145{
1146 //
1147 // Remapping track reference
1148 // Called at finish primary
1149 //
3cac0921 1150 if (!fTmpTrackReferences) return;
1151 Int_t nEntries = fTmpTrackReferences->GetEntries();
4e490a96 1152 for (Int_t i=0; i < nEntries; i++){
3cac0921 1153 AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTmpTrackReferences->UncheckedAt(i));
4e490a96 1154 if (ref) {
1155 Int_t newID = map[ref->GetTrack()];
1156 if (newID>=0) ref->SetTrack(newID);
1157 else {
1158 ref->SetBit(kNotDeleted,kFALSE);
3cac0921 1159 fTmpTrackReferences->RemoveAt(i);
4e490a96 1160 }
1161 } // if ref
5d12ce38 1162 }
3cac0921 1163 fTmpTrackReferences->Compress();
5d12ce38 1164}
0561efeb 1165
1166void AliMC::FixParticleDecaytime()
1167{
1168 //
1169 // Fix the particle decay time according to rmin and rmax for decays
1170 //
1171
1172 TLorentzVector p;
1173 gMC->TrackMomentum(p);
1174 Double_t tmin, tmax;
1175 Double_t b;
1176
1177 // Transverse velocity
1178 Double_t vt = p.Pt() / p.E();
1179
1180 if ((b = gAlice->Field()->SolenoidField()) > 0.) { // [kG]
1181
1182 // Radius of helix
1183
1184 Double_t rho = p.Pt() / 0.0003 / b; // [cm]
1185
1186 // Revolution frequency
1187
1188 Double_t omega = vt / rho;
1189
1190 // Maximum and minimum decay time
1191 //
1192 // Check for curlers first
1193 if (fRDecayMax * fRDecayMax / rho / rho / 2. > 1.) return;
1194
1195 //
1196
1197 tmax = TMath::ACos(1. - fRDecayMax * fRDecayMax / rho / rho / 2.) / omega; // [ct]
1198 tmin = TMath::ACos(1. - fRDecayMin * fRDecayMin / rho / rho / 2.) / omega; // [ct]
1199 } else {
1200 tmax = fRDecayMax / vt; // [ct]
1201 tmin = fRDecayMin / vt; // [ct]
1202 }
1203
1204 //
1205 // Dial t using the two limits
1206 Double_t t = tmin + (tmax - tmin) * gRandom->Rndm(); // [ct]
1207 //
1208 //
1209 // Force decay time in transport code
1210 //
fa7ca4da 1211 gMC->ForceDecayTime(t / 2.99792458e10);
0561efeb 1212}
3cac0921 1213
1214void AliMC::MakeTmpTrackRefsTree()
1215{
1216 // Make the temporary track reference tree
1217 fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
1218 fTmpTreeTR = new TTree("TreeTR", "Track References");
1219 if (!fTmpTrackReferences) fTmpTrackReferences = new TClonesArray("AliTrackReference", 100);
1220 fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTmpTrackReferences, 4000);
1221}
1222
1223void AliMC::ReorderAndExpandTreeTR()
1224{
1225//
1226// Reorder and expand the temporary track reference tree in order to match the kinematics tree
1227//
1228
1229 AliRunLoader *rl = gAlice->GetRunLoader();
1230//
1231// TreeTR
1232 AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
1233 rl->MakeTrackRefsContainer();
1234 TTree * treeTR = rl->TreeTR();
1235 if (treeTR){
1236 // make branch for central track references
1237 if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference",0);
1238 TBranch *branch;
1239 branch = treeTR->Branch("TrackReferences",&fTrackReferences);
1240 branch->SetAddress(&fTrackReferences);
1241 }
1242
1243 AliStack* stack = rl->Stack();
1244 Int_t np = stack->GetNprimary();
1245 Int_t nt = fTmpTreeTR->GetEntries();
3cac0921 1246 //
1247 // Loop over tracks and find the secondaries with the help of the kine tree
5ce21fc4 1248 Int_t ifills = 0;
1249 Int_t it = 0;
1250 for (Int_t ip = np - 1; ip > -1; ip--) {
1251 TParticle *part = stack->Particle(ip);
2c4e7c77 1252// printf("Particle %5d %5d %5d %5d %5d \n", ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), part->GetLastDaughter());
1253
5ce21fc4 1254 // Skip primaries that have not been transported
3cac0921 1255 Int_t dau1 = part->GetFirstDaughter();
1256 Int_t dau2 = -1;
f4414561 1257 // if ((dau1 > -1 && dau1 < np) || part->GetStatusCode() > 1) continue;
1258 if (!part->TestBit(kTransportBit)) continue;
5ce21fc4 1259 //
1260 fTmpTreeTR->GetEntry(it++);
1261 Int_t nh = fTmpTrackReferences->GetEntries();
1262 // Determine range of secondaries produced by this primary
3cac0921 1263 if (dau1 > -1) {
5ce21fc4 1264 Int_t inext = ip - 1;
3cac0921 1265 while (dau2 < 0) {
1266 if (inext >= 0) {
1267 part = stack->Particle(inext);
5ce21fc4 1268 dau2 = part->GetFirstDaughter();
1269 if (dau2 == -1 || dau2 < np) {
1270 dau2 = -1;
1271 } else {
1272 dau2--;
1273 }
3cac0921 1274 } else {
1275 dau2 = stack->GetNtrack() - 1;
1276 }
1277 inext--;
5ce21fc4 1278 } // find upper bound
1279 } // dau2 < 0
2c4e7c77 1280// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
1281
3cac0921 1282 //
1283 // Loop over reference hits and find secondary label
1284 for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
1285 for (Int_t ih = 0; ih < nh; ih++) {
1286 AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences->At(ih);
1287 Int_t label = tr->Label();
5ce21fc4 1288 // Skip primaries
1289 if (label == ip) continue;
2c4e7c77 1290 if (label > dau2 || label < dau1)
1291 AliWarning(Form("Track Reference Label out of range !: %5d %5d %5d \n", label, dau1, dau2));
3cac0921 1292 if (label == id) {
1293 // secondary found
1294 Int_t nref = fTrackReferences->GetEntriesFast();
1295 TClonesArray &lref = *fTrackReferences;
1296 new(lref[nref]) AliTrackReference(*tr);
1297 }
1298 } // hits
1299 treeTR->Fill();
1300 fTrackReferences->Clear();
5ce21fc4 1301 ifills++;
3cac0921 1302 } // daughters
1303 } // tracks
3cac0921 1304 //
1305 // Now loop again and write the primaries
5ce21fc4 1306 it = nt - 1;
1307 for (Int_t ip = 0; ip < np; ip++) {
1308 TParticle* part = stack->Particle(ip);
f4414561 1309// if ((part->GetFirstDaughter() == -1 && part->GetStatusCode() <= 1) || part->GetFirstDaughter() >= np)
c226988d 1310 if (part->TestBit(kTransportBit))
f4414561 1311 {
5ce21fc4 1312 // Skip particles that have not been transported
1313 fTmpTreeTR->GetEntry(it--);
1314 Int_t nh = fTmpTrackReferences->GetEntries();
1315 //
1316 // Loop over reference hits and find primary labels
1317 for (Int_t ih = 0; ih < nh; ih++) {
1318 AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences->At(ih);
1319 Int_t label = tr->Label();
1320 if (label == ip) {
1321 Int_t nref = fTrackReferences->GetEntriesFast();
1322 TClonesArray &lref = *fTrackReferences;
1323 new(lref[nref]) AliTrackReference(*tr);
1324 }
1325 }
3cac0921 1326 }
1327 treeTR->Fill();
1328 fTrackReferences->Clear();
5ce21fc4 1329 ifills++;
3cac0921 1330 } // tracks
5ce21fc4 1331 // Check
1332 if (ifills != stack->GetNtrack())
1333 AliWarning(Form("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, stack->GetNtrack()));
3cac0921 1334//
1335// Clean-up
1336 delete fTmpTreeTR;
1337 fTmpFileTR->Close();
1338 delete fTmpFileTR;
1339 delete fTmpTrackReferences;
1340 fTmpTrackReferences = 0;
1341 gSystem->Exec("rm -rf TrackRefsTmp.root");
1342}
1343