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