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