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