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