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