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