Implementation of copy constructor and asignment operator
[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
18#include <TBrowser.h>
19#include <TStopwatch.h>
20#include <TSystem.h>
21#include <TVirtualMC.h>
22
23#include "AliDetector.h"
24#include "AliGenerator.h"
25#include "AliHeader.h"
26#include "AliLego.h"
27#include "AliMC.h"
28#include "AliMCQA.h"
29#include "AliRun.h"
30#include "AliStack.h"
31#include "AliTrackReference.h"
32
33
34ClassImp(AliMC)
35
36//_______________________________________________________________________
37AliMC::AliMC() :
38 fGenerator(0),
39 fEventEnergy(0),
40 fSummEnergy(0),
41 fSum2Energy(0),
42 fTrRmax(1.e10),
43 fTrZmax(1.e10),
44 fImedia(0),
45 fTransParName("\0"),
46 fMCQA(0),
47 fHitLists(0),
48 fTrackReferences(0)
49
50{
51}
52
53//_______________________________________________________________________
54AliMC::AliMC(const char *name, const char *title) :
55 TVirtualMCApplication(name, title),
56 fGenerator(0),
57 fEventEnergy(0),
58 fSummEnergy(0),
59 fSum2Energy(0),
60 fTrRmax(1.e10),
61 fTrZmax(1.e10),
62 fImedia(new TArrayI(1000)),
63 fTransParName("\0"),
64 fMCQA(0),
65 fHitLists(new TList()),
66 fTrackReferences(new TClonesArray("AliTrackReference", 100))
67{
68
69 // Set transport parameters
70 SetTransPar();
71
72 // Prepare the tracking medium lists
73 for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
74
75}
76
77//_______________________________________________________________________
78AliMC::AliMC(const AliMC &mc) :
79 TVirtualMCApplication(mc),
80 fGenerator(0),
81 fEventEnergy(0),
82 fSummEnergy(0),
83 fSum2Energy(0),
84 fTrRmax(1.e10),
85 fTrZmax(1.e10),
86 fImedia(0),
87 fTransParName("\0"),
88 fMCQA(0),
89 fHitLists(0),
90 fTrackReferences(0)
91{
92 //
93 // Copy constructor for AliRun
94 //
95 mc.Copy(*this);
96}
97
98//_______________________________________________________________________
99AliMC::~AliMC()
100{
101 delete fGenerator;
102 delete fImedia;
103 delete fMCQA;
104 delete fHitLists;
105 // Delete track references
106 if (fTrackReferences) {
107 fTrackReferences->Delete();
108 delete fTrackReferences;
109 fTrackReferences = 0;
110 }
111
112}
113
114//_______________________________________________________________________
115void AliMC::Copy(AliMC &) const
116{
117 Fatal("Copy","Not implemented!\n");
118}
119
120//_______________________________________________________________________
121void AliMC::ConstructGeometry()
122{
123 //
124 // Create modules, materials, geometry
125 //
126
127 TStopwatch stw;
128 TIter next(gAlice->Modules());
129 AliModule *detector;
130 if (GetDebug()) Info("ConstructGeometry","Geometry creation:");
131 while((detector = dynamic_cast<AliModule*>(next()))) {
132 stw.Start();
133 // Initialise detector materials and geometry
134 detector->CreateMaterials();
135 detector->CreateGeometry();
136 printf("%10s R:%.2fs C:%.2fs\n",
137 detector->GetName(),stw.RealTime(),stw.CpuTime());
138 }
139}
140
141//_______________________________________________________________________
142void AliMC::InitGeometry()
143{
144 //
145 // Initialize detectors and display geometry
146 //
147
148 printf("Initialisation:\n");
149 TStopwatch stw;
150 TIter next(gAlice->Modules());
151 AliModule *detector;
152 while((detector = dynamic_cast<AliModule*>(next()))) {
153 stw.Start();
154 // Initialise detector and display geometry
155 detector->Init();
156 detector->BuildGeometry();
157 printf("%10s R:%.2fs C:%.2fs\n",
158 detector->GetName(),stw.RealTime(),stw.CpuTime());
159 }
160
161}
162
163//_______________________________________________________________________
164void AliMC::GeneratePrimaries()
165{
166 //
167 // Generate primary particles and fill them in the stack.
168 //
169
170 Generator()->Generate();
171}
172
173//_______________________________________________________________________
174void AliMC::SetGenerator(AliGenerator *generator)
175{
176 //
177 // Load the event generator
178 //
179 if(!fGenerator) fGenerator = generator;
180}
181
182//_______________________________________________________________________
183void AliMC::ResetGenerator(AliGenerator *generator)
184{
185 //
186 // Load the event generator
187 //
188 if(fGenerator)
189 if(generator)
190 Warning("ResetGenerator","Replacing generator %s with %s\n",
191 fGenerator->GetName(),generator->GetName());
192 else
193 Warning("ResetGenerator","Replacing generator %s with NULL\n",
194 fGenerator->GetName());
195 fGenerator = generator;
196}
197
198//_______________________________________________________________________
199void AliMC::FinishRun()
200{
201 // Clean generator information
202 if (GetDebug()) Info("FinishRun"," fGenerator->FinishRun()");
203 fGenerator->FinishRun();
204
205 //Output energy summary tables
206 if (GetDebug()) Info("FinishRun"," EnergySummary()");
207 EnergySummary();
208
209}
210
211//_______________________________________________________________________
212void AliMC::BeginPrimary()
213{
214 //
215 // Called at the beginning of each primary track
216 //
217
218 // Reset Hits info
219 ResetHits();
220 ResetTrackReferences();
221
222}
223
224//_______________________________________________________________________
225void AliMC::PreTrack()
226{
227 TObjArray &dets = *gAlice->Modules();
228 AliModule *module;
229
230 for(Int_t i=0; i<=gAlice->GetNdets(); i++)
231 if((module = dynamic_cast<AliModule*>(dets[i])))
232 module->PreTrack();
233
234 fMCQA->PreTrack();
235}
236
237//_______________________________________________________________________
238void AliMC::Stepping()
239{
240 //
241 // Called at every step during transport
242 //
243
244 Int_t id = DetFromMate(gMC->GetMedium());
245 if (id < 0) return;
246
247 //
248 // --- If lego option, do it and leave
249 if (gAlice->Lego())
250 gAlice->Lego()->StepManager();
251 else {
252 Int_t copy;
253 //Update energy deposition tables
254 AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
255
256 //Call the appropriate stepping routine;
257 AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
258 if(det && det->StepManagerIsEnabled()) {
259 fMCQA->StepManager(id);
260 det->StepManager();
261 }
262 }
263}
264
265//_______________________________________________________________________
266void AliMC::EnergySummary()
267{
268 //
269 // Print summary of deposited energy
270 //
271
272 Int_t ndep=0;
273 Float_t edtot=0;
274 Float_t ed, ed2;
275 Int_t kn, i, left, j, id;
276 const Float_t kzero=0;
277 Int_t ievent=gAlice->GetRunLoader()->GetHeader()->GetEvent()+1;
278 //
279 // Energy loss information
280 if(ievent) {
281 printf("***************** Energy Loss Information per event (GEV) *****************\n");
282 for(kn=1;kn<fEventEnergy.GetSize();kn++) {
283 ed=fSummEnergy[kn];
284 if(ed>0) {
285 fEventEnergy[ndep]=kn;
286 if(ievent>1) {
287 ed=ed/ievent;
288 ed2=fSum2Energy[kn];
289 ed2=ed2/ievent;
290 ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
291 } else
292 ed2=99;
293 fSummEnergy[ndep]=ed;
294 fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
295 edtot+=ed;
296 ndep++;
297 }
298 }
299 for(kn=0;kn<(ndep-1)/3+1;kn++) {
300 left=ndep-kn*3;
301 for(i=0;i<(3<left?3:left);i++) {
302 j=kn*3+i;
303 id=Int_t (fEventEnergy[j]+0.1);
304 printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
305 }
306 printf("\n");
307 }
308 //
309 // Relative energy loss in different detectors
310 printf("******************** Relative Energy Loss per event ********************\n");
311 printf("Total energy loss per event %10.3f GeV\n",edtot);
312 for(kn=0;kn<(ndep-1)/5+1;kn++) {
313 left=ndep-kn*5;
314 for(i=0;i<(5<left?5:left);i++) {
315 j=kn*5+i;
316 id=Int_t (fEventEnergy[j]+0.1);
317 printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
318 }
319 printf("\n");
320 }
321 for(kn=0;kn<75;kn++) printf("*");
322 printf("\n");
323 }
324 //
325 // Reset the TArray's
326 // fEventEnergy.Set(0);
327 // fSummEnergy.Set(0);
328 // fSum2Energy.Set(0);
329}
330
331//_____________________________________________________________________________
332void AliMC::BeginEvent()
333{
334 //
335 // Clean-up previous event
336 // Energy scores
337 if (GetDebug())
338 {
339 Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
340 Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
341 Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
342 Info("BeginEvent"," BEGINNING EVENT ");
343 Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
344 Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
345 Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
346 }
347
348 AliRunLoader *runloader=gAlice->GetRunLoader();
349
350 /*******************************/
351 /* Clean after eventual */
352 /* previous event */
353 /*******************************/
354
355
356 //Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors),
357 gAlice->SetEventNrInRun(gAlice->GetEventNrInRun()+1);
358 runloader->SetEventNumber(gAlice->GetEventNrInRun());// sets new files, cleans the previous event stuff, if necessary, etc.,
359 if (GetDebug()) Info("BeginEvent","EventNr is %d",gAlice->GetEventNrInRun());
360
361 fEventEnergy.Reset();
362 // Clean detector information
363
364 if (runloader->Stack())
365 runloader->Stack()->Reset();//clean stack -> tree is unloaded
366 else
367 runloader->MakeStack();//or make a new one
504b172d 368
369 if(gAlice->Lego() == 0x0)
370 {
371 if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTree(K)");
372 runloader->MakeTree("K");
373 }
374
5d12ce38 375 if (GetDebug()) Info("BeginEvent"," gMC->SetStack(fRunLoader->Stack())");
376 gMC->SetStack(gAlice->GetRunLoader()->Stack());//Was in InitMC - but was moved here
377 //because we don't have guarantee that
378 //stack pointer is not going to change from event to event
379 //since it bellobgs to header and is obtained via RunLoader
380 //
381 // Reset all Detectors & kinematics & make/reset trees
382 //
383
384 runloader->GetHeader()->Reset(gAlice->GetRunNumber(),gAlice->GetEvNumber(),
385 gAlice->GetEventNrInRun());
386// fRunLoader->WriteKinematics("OVERWRITE"); is there any reason to rewrite here since MakeTree does so
387
504b172d 388 if(gAlice->Lego())
389 {
390 gAlice->Lego()->BeginEvent();
391 return;
392 }
393
5d12ce38 394
395 if (GetDebug()) Info("BeginEvent"," ResetHits()");
396 ResetHits();
504b172d 397
5d12ce38 398 if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTree(H)");
399 runloader->MakeTree("H");
504b172d 400
401 if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTrackRefsContainer()");
402 runloader->MakeTrackRefsContainer();//for insurance
5d12ce38 403
5d12ce38 404
405 //create new branches and SetAdresses
406 TIter next(gAlice->Modules());
407 AliModule *detector;
408 while((detector = (AliModule*)next()))
409 {
410 if (GetDebug()) Info("BeginEvent"," %s->MakeBranch(H)",detector->GetName());
411 detector->MakeBranch("H");
412 if (GetDebug()) Info("BeginEvent"," %s->MakeBranchTR()",detector->GetName());
413 detector->MakeBranchTR();
414 if (GetDebug()) Info("BeginEvent"," %s->SetTreeAddress()",detector->GetName());
415 detector->SetTreeAddress();
416 }
417}
418
419//_______________________________________________________________________
420void AliMC::ResetHits()
421{
422 //
423 // Reset all Detectors hits
424 //
425 TIter next(gAlice->Modules());
426 AliModule *detector;
427 while((detector = dynamic_cast<AliModule*>(next()))) {
428 detector->ResetHits();
429 }
430}
431
432//_______________________________________________________________________
433void AliMC::PostTrack()
434{
435 TObjArray &dets = *gAlice->Modules();
436 AliModule *module;
437
438 for(Int_t i=0; i<=gAlice->GetNdets(); i++)
439 if((module = dynamic_cast<AliModule*>(dets[i])))
440 module->PostTrack();
441}
442
443//_______________________________________________________________________
444void AliMC::FinishPrimary()
445{
446 //
447 // Called at the end of each primary track
448 //
449 AliRunLoader *runloader=gAlice->GetRunLoader();
450 // static Int_t count=0;
451 // const Int_t times=10;
452 // This primary is finished, purify stack
453 runloader->Stack()->PurifyKine();
504b172d 454
5d12ce38 455 TIter next(gAlice->Modules());
456 AliModule *detector;
457 while((detector = dynamic_cast<AliModule*>(next()))) {
458 detector->FinishPrimary();
504b172d 459 AliLoader* loader = detector->GetLoader();
460 if(loader)
461 {
462 TTree* treeH = loader->TreeH();
463 if (treeH) treeH->Fill(); //can be Lego run and treeH can not exist
464 }
5d12ce38 465 }
466
467 // Write out track references if any
468 if (runloader->TreeTR()) runloader->TreeTR()->Fill();
469}
470
471//_______________________________________________________________________
472void AliMC::FinishEvent()
473{
474 //
475 // Called at the end of the event.
476 //
477
478 //
479 if(gAlice->Lego()) gAlice->Lego()->FinishEvent();
480
481 TIter next(gAlice->Modules());
482 AliModule *detector;
483 while((detector = dynamic_cast<AliModule*>(next()))) {
484 detector->FinishEvent();
485 }
486
487 //Update the energy deposit tables
488 Int_t i;
489 for(i=0;i<fEventEnergy.GetSize();i++)
490 {
491 fSummEnergy[i]+=fEventEnergy[i];
492 fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
493 }
494
495 AliRunLoader *runloader=gAlice->GetRunLoader();
496
497 AliHeader* header = runloader->GetHeader();
498 AliStack* stack = runloader->Stack();
499 if ( (header == 0x0) || (stack == 0x0) )
500 {//check if we got header and stack. If not cry and exit aliroot
501 Fatal("AliRun","Can not get the stack or header from LOADER");
502 return;//never reached
503 }
504 // Update Header information
505 header->SetNprimary(stack->GetNprimary());
506 header->SetNtrack(stack->GetNtrack());
507
508
509 // Write out the kinematics
510 stack->FinishEvent();
511
512 // Write out the event Header information
513 TTree* treeE = runloader->TreeE();
514 if (treeE)
515 {
516 header->SetStack(stack);
517 treeE->Fill();
518 }
519 else
520 {
521 Error("FinishEvent","Can not get TreeE from RL");
522 }
523
504b172d 524 if(gAlice->Lego() == 0x0)
525 {
526 runloader->WriteKinematics("OVERWRITE");
527 runloader->WriteTrackRefs("OVERWRITE");
528 runloader->WriteHits("OVERWRITE");
529 }
530
5d12ce38 531 if (GetDebug())
532 {
533 Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
534 Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
535 Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
536 Info("FinishEvent"," FINISHING EVENT ");
537 Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
538 Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
539 Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
540 }
541}
542
543//_______________________________________________________________________
544void AliMC::Field(const Double_t* x, Double_t* b) const
545{
546 gAlice->Field(x,b);
547}
548
549//_______________________________________________________________________
550void AliMC::Init()
551{
552
553 //=================Create Materials and geometry
554 gMC->Init();
555 gMC->DefineParticles(); //Create standard MC particles
556
557 //Read the cuts for all materials
558 ReadTransPar();
559 //Build the special IMEDIA table
560 MediaTable();
561
562 //Compute cross-sections
563 gMC->BuildPhysics();
564
565 //Write Geometry object to current file.
566 gAlice->GetRunLoader()->WriteGeometry();
567
568 //Initialise geometry deposition table
569 fEventEnergy.Set(gMC->NofVolumes()+1);
570 fSummEnergy.Set(gMC->NofVolumes()+1);
571 fSum2Energy.Set(gMC->NofVolumes()+1);
572
573 //
574 fMCQA = new AliMCQA(gAlice->GetNdets());
575
576 // Register MC in configuration
577 AliConfig::Instance()->Add(gMC);
578
579}
580
581//_______________________________________________________________________
582void AliMC::MediaTable()
583{
584 //
585 // Built media table to get from the media number to
586 // the detector id
587 //
588
589 Int_t kz, nz, idt, lz, i, k, ind;
590 // Int_t ibeg;
591 TObjArray &dets = *gAlice->Detectors();
592 AliModule *det;
593 Int_t ndets=gAlice->GetNdets();
594 //
595 // For all detectors
596 for (kz=0;kz<ndets;kz++) {
597 // If detector is defined
598 if((det=dynamic_cast<AliModule*>(dets[kz]))) {
599 TArrayI &idtmed = *(det->GetIdtmed());
600 for(nz=0;nz<100;nz++) {
601 // Find max and min material number
602 if((idt=idtmed[nz])) {
603 det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
604 det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
605 }
606 }
607 if(det->LoMedium() > det->HiMedium()) {
608 det->LoMedium() = 0;
609 det->HiMedium() = 0;
610 } else {
611 if(det->HiMedium() > fImedia->GetSize()) {
612 Error("MediaTable","Increase fImedia from %d to %d",
613 fImedia->GetSize(),det->HiMedium());
614 return;
615 }
616 // Tag all materials in rage as belonging to detector kz
617 for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
618 (*fImedia)[lz]=kz;
619 }
620 }
621 }
622 }
623 //
624 // Print summary table
625 printf(" Traking media ranges:\n");
626 for(i=0;i<(ndets-1)/6+1;i++) {
627 for(k=0;k< (6<ndets-i*6?6:ndets-i*6);k++) {
628 ind=i*6+k;
629 det=dynamic_cast<AliModule*>(dets[ind]);
630 if(det)
631 printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
632 det->HiMedium());
633 else
634 printf(" %6s: %3d -> %3d;","NULL",0,0);
635 }
636 printf("\n");
637 }
638}
639
640//_______________________________________________________________________
641void AliMC::ReadTransPar()
642{
643 //
644 // Read filename to set the transport parameters
645 //
646
647
648 const Int_t kncuts=10;
649 const Int_t knflags=11;
650 const Int_t knpars=kncuts+knflags;
651 const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
652 "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
653 "BREM","COMP","DCAY","DRAY","HADR","LOSS",
654 "MULS","PAIR","PHOT","RAYL"};
655 char line[256];
656 char detName[7];
657 char* filtmp;
658 Float_t cut[kncuts];
659 Int_t flag[knflags];
660 Int_t i, itmed, iret, ktmed, kz;
661 FILE *lun;
662 //
663 // See whether the file is there
664 filtmp=gSystem->ExpandPathName(fTransParName.Data());
665 lun=fopen(filtmp,"r");
666 delete [] filtmp;
667 if(!lun) {
668 Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
669 return;
670 }
671 //
672 if(fDebug) {
673 printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
674 printf(" *%59s\n","*");
675 printf(" * Please check carefully what you are doing!%10s\n","*");
676 printf(" *%59s\n","*");
677 }
678 //
679 while(1) {
680 // Initialise cuts and flags
681 for(i=0;i<kncuts;i++) cut[i]=-99;
682 for(i=0;i<knflags;i++) flag[i]=-99;
683 itmed=0;
684 for(i=0;i<256;i++) line[i]='\0';
685 // Read up to the end of line excluded
686 iret=fscanf(lun,"%[^\n]",line);
687 if(iret<0) {
688 //End of file
689 fclose(lun);
690 if(fDebug){
691 printf(" *%59s\n","*");
692 printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
693 }
694 return;
695 }
696 // Read the end of line
697 fscanf(lun,"%*c");
698 if(!iret) continue;
699 if(line[0]=='*') continue;
700 // Read the numbers
701 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",
702 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
703 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
704 &flag[8],&flag[9],&flag[10]);
705 if(!iret) continue;
706 if(iret<0) {
707 //reading error
708 Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
709 continue;
710 }
711 // Check that the module exist
712 AliModule *mod = gAlice->GetModule(detName);
713 if(mod) {
714 // Get the array of media numbers
715 TArrayI &idtmed = *mod->GetIdtmed();
716 // Check that the tracking medium code is valid
717 if(0<=itmed && itmed < 100) {
718 ktmed=idtmed[itmed];
719 if(!ktmed) {
720 Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
721 continue;
722 }
723 // Set energy thresholds
724 for(kz=0;kz<kncuts;kz++) {
725 if(cut[kz]>=0) {
726 if(fDebug) printf(" * %-6s set to %10.3E for tracking medium code %4d for %s\n",
727 kpars[kz],cut[kz],itmed,mod->GetName());
728 gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
729 }
730 }
731 // Set transport mechanisms
732 for(kz=0;kz<knflags;kz++) {
733 if(flag[kz]>=0) {
734 if(fDebug) printf(" * %-6s set to %10d for tracking medium code %4d for %s\n",
735 kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
736 gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
737 }
738 }
739 } else {
740 Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
741 continue;
742 }
743 } else {
744 if(fDebug) printf("%s::ReadTransParModule: %s not present\n",ClassName(),detName);
745 continue;
746 }
747 }
748}
749
750//_______________________________________________________________________
751void AliMC::SetTransPar(const char *filename)
752{
753 //
754 // Sets the file name for transport parameters
755 //
756 fTransParName = filename;
757}
758
759//_______________________________________________________________________
760void AliMC::Browse(TBrowser *b)
761{
762 //
763 // Called when the item "Run" is clicked on the left pane
764 // of the Root browser.
765 // It displays the Root Trees and all detectors.
766 //
767 //detectors are in folders anyway
768 b->Add(fMCQA,"AliMCQA");
769}
770
771//PH
772//_______________________________________________________________________
773void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
774{
775 //
776 // Add a hit to detector id
777 //
778 TObjArray &dets = *gAlice->Modules();
779 if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
780}
781
782//_______________________________________________________________________
783void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
784{
785 //
786 // Add digit to detector id
787 //
788 TObjArray &dets = *gAlice->Modules();
789 if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
790}
791
792//_______________________________________________________________________
793Int_t AliMC::GetCurrentTrackNumber() const {
794 //
795 // Returns current track
796 //
797 return gAlice->GetRunLoader()->Stack()->GetCurrentTrackNumber();
798}
799
800//_______________________________________________________________________
801void AliMC::DumpPart (Int_t i) const
802{
803 //
804 // Dumps particle i in the stack
805 //
806 AliRunLoader * runloader = gAlice->GetRunLoader();
807 if (runloader->Stack())
808 runloader->Stack()->DumpPart(i);
809}
810
811//_______________________________________________________________________
812void AliMC::DumpPStack () const
813{
814 //
815 // Dumps the particle stack
816 //
817 AliRunLoader * runloader = gAlice->GetRunLoader();
818 if (runloader->Stack())
819 runloader->Stack()->DumpPStack();
820}
821
822//_______________________________________________________________________
823Int_t AliMC::GetNtrack() const {
824 //
825 // Returns number of tracks in stack
826 //
827 Int_t ntracks = -1;
828 AliRunLoader * runloader = gAlice->GetRunLoader();
829 if (runloader->Stack())
830 ntracks = runloader->Stack()->GetNtrack();
831 return ntracks;
832}
833
834//_______________________________________________________________________
835Int_t AliMC::GetPrimary(Int_t track) const
836{
837 //
838 // return number of primary that has generated track
839 //
840 Int_t nprimary = -999;
841 AliRunLoader * runloader = gAlice->GetRunLoader();
842 if (runloader->Stack())
843 nprimary = runloader->Stack()->GetPrimary(track);
844 return nprimary;
845}
846
847//_______________________________________________________________________
848TParticle* AliMC::Particle(Int_t i) const
849{
850 AliRunLoader * runloader = gAlice->GetRunLoader();
851 if (runloader)
852 if (runloader->Stack())
853 return runloader->Stack()->Particle(i);
854 return 0x0;
855}
856
857//_______________________________________________________________________
858TObjArray* AliMC::Particles() const {
859 //
860 // Returns pointer to Particles array
861 //
862 AliRunLoader * runloader = gAlice->GetRunLoader();
863 if (runloader)
864 if (runloader->Stack())
865 return runloader->Stack()->Particles();
866 return 0x0;
867}
868
869//_______________________________________________________________________
870void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
871 Float_t *vpos, Float_t *polar, Float_t tof,
872 TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
873{
874// Delegate to stack
875//
876 AliRunLoader * runloader = gAlice->GetRunLoader();
877 if (runloader)
878 if (runloader->Stack())
879 runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
880 mech, ntr, weight, is);
881}
882
883//_______________________________________________________________________
884void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg,
885 Double_t px, Double_t py, Double_t pz, Double_t e,
886 Double_t vx, Double_t vy, Double_t vz, Double_t tof,
887 Double_t polx, Double_t poly, Double_t polz,
888 TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
889{
890 // Delegate to stack
891 //
892 AliRunLoader * runloader = gAlice->GetRunLoader();
893 if (runloader)
894 if (runloader->Stack())
895 runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
896 polx, poly, polz, mech, ntr, weight, is);
897}
898
899//_______________________________________________________________________
900void AliMC::SetHighWaterMark(const Int_t nt)
901{
902 //
903 // Set high water mark for last track in event
904 AliRunLoader * runloader = gAlice->GetRunLoader();
905 if (runloader)
906 if (runloader->Stack())
907 runloader->Stack()->SetHighWaterMark(nt);
908}
909
910//_______________________________________________________________________
911void AliMC::KeepTrack(const Int_t track)
912{
913 //
914 // Delegate to stack
915 //
916 AliRunLoader * runloader = gAlice->GetRunLoader();
917 if (runloader)
918 if (runloader->Stack())
919 runloader->Stack()->KeepTrack(track);
920}
921
922//_______________________________________________________________________
923void AliMC::FlagTrack(Int_t track)
924{
925 // Delegate to stack
926 //
927 AliRunLoader * runloader = gAlice->GetRunLoader();
928 if (runloader)
929 if (runloader->Stack())
930 runloader->Stack()->FlagTrack(track);
931}
932
933//_______________________________________________________________________
934void AliMC::SetCurrentTrack(Int_t track)
935{
936 //
937 // Set current track number
938 //
939 AliRunLoader * runloader = gAlice->GetRunLoader();
940 if (runloader)
941 if (runloader->Stack())
942 runloader->Stack()->SetCurrentTrack(track);
943}
944
945//_______________________________________________________________________
946void AliMC::AddTrackReference(Int_t label){
947 //
948 // add a trackrefernce to the list
949 if (!fTrackReferences) {
950 cerr<<"Container trackrefernce not active\n";
951 return;
952 }
953 Int_t nref = fTrackReferences->GetEntriesFast();
954 TClonesArray &lref = *fTrackReferences;
955 new(lref[nref]) AliTrackReference(label);
956}
957
958
959
960//_______________________________________________________________________
961void AliMC::ResetTrackReferences()
962{
963 //
964 // Reset all references
965 //
966 if (fTrackReferences) fTrackReferences->Clear();
967
968 TIter next(gAlice->Modules());
969 AliModule *detector;
970 while((detector = dynamic_cast<AliModule*>(next()))) {
971 detector->ResetTrackReferences();
972 }
973}
974
975void AliMC::RemapTrackReferencesIDs(Int_t *map)
976{
977 //
978 // Remapping track reference
979 // Called at finish primary
980 //
981 if (!fTrackReferences) return;
982 for (Int_t i=0;i<fTrackReferences->GetEntries();i++){
983 AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTrackReferences->UncheckedAt(i));
984 if (ref) {
985 Int_t newID = map[ref->GetTrack()];
986 if (newID>=0) ref->SetTrack(newID);
987 else {
988 //ref->SetTrack(-1);
989 ref->SetBit(kNotDeleted,kFALSE);
990 fTrackReferences->RemoveAt(i);
991 }
992 }
993 }
994 fTrackReferences->Compress();
995}