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