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