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