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