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