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