7e9c6b53e0bb4408a57d5e1cf81aa81a9a7fa8e5
[u/mrichter/AliRoot.git] / STEER / AliRun.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 /*
17 $Log$
18 Revision 1.62  2001/04/06 11:12:33  morsch
19 Clear fParticles after each event. (Ivana Hrivnacova)
20
21 Revision 1.61  2001/03/30 07:04:10  morsch
22 Call fGenerator->FinishRun() for final print-outs, cross-section and weight calculations.
23
24 Revision 1.60  2001/03/21 18:22:30  hristov
25 fParticleFileMap fix (I.Hrivnacova)
26
27 Revision 1.59  2001/03/12 17:47:03  hristov
28 Changes needed on Sun with CC 5.0
29
30 Revision 1.58  2001/03/09 14:27:26  morsch
31 Fix for multiple events per file: inhibit decrease of size of  fParticleFileMap.
32
33 Revision 1.57  2001/02/23 17:40:23  buncic
34 All trees needed for simulation created in RunMC(). TreeR and its branches
35 are now created in new RunReco() method.
36
37 Revision 1.56  2001/02/14 15:45:20  hristov
38 Algorithmic way of getting entry index in fParticleMap. Protection of fParticleFileMap (I.Hrivnacova)
39
40 Revision 1.55  2001/02/12 15:52:54  buncic
41 Removed OpenBaseFile().
42
43 Revision 1.54  2001/02/07 10:39:05  hristov
44 Remove default value for argument
45
46 Revision 1.53  2001/02/06 11:02:26  hristov
47 New SetTrack interface added, added check for unfilled particles in FinishEvent (I.Hrivnacova)
48
49 Revision 1.52  2001/02/05 16:22:25  buncic
50 Added TreeS to GetEvent().
51
52 Revision 1.51  2001/02/02 15:16:20  morsch
53 SetHighWaterMark method added to mark last particle in event.
54
55 Revision 1.50  2001/01/27 10:32:00  hristov
56 Leave the loop when primaries are filled (I.Hrivnacova)
57
58 Revision 1.49  2001/01/26 19:58:48  hristov
59 Major upgrade of AliRoot code
60
61 Revision 1.48  2001/01/17 10:50:50  hristov
62 Corrections to destructors
63
64 Revision 1.47  2000/12/18 10:44:01  morsch
65 Possibility to set field map by passing pointer to objet of type AliMagF via
66 SetField().
67 Example:
68 gAlice->SetField(new AliMagFCM("Map2", "$(ALICE_ROOT)/data/field01.dat",2,1.,10.));
69
70 Revision 1.46  2000/12/14 19:29:27  fca
71 galice.cuts was not read any more
72
73 Revision 1.45  2000/11/30 07:12:49  alibrary
74 Introducing new Rndm and QA classes
75
76 Revision 1.44  2000/10/26 13:58:59  morsch
77 Add possibility to choose the lego generator (of type AliGeneratorLego or derived) when running
78 RunLego(). Default is the base class AliGeneratorLego.
79
80 Revision 1.43  2000/10/09 09:43:17  fca
81 Special remapping of hits for TPC and TRD. End-of-primary action introduced
82
83 Revision 1.42  2000/10/02 21:28:14  fca
84 Removal of useless dependecies via forward declarations
85
86 Revision 1.41  2000/07/13 16:19:09  fca
87 Mainly coding conventions + some small bug fixes
88
89 Revision 1.40  2000/07/12 08:56:25  fca
90 Coding convention correction and warning removal
91
92 Revision 1.39  2000/07/11 18:24:59  fca
93 Coding convention corrections + few minor bug fixes
94
95 Revision 1.38  2000/06/20 13:05:45  fca
96 Writing down the TREE headers before job starts
97
98 Revision 1.37  2000/06/09 20:05:11  morsch
99 Introduce possibility to chose magnetic field version 3: AliMagFDM + field02.dat
100
101 Revision 1.36  2000/06/08 14:03:58  hristov
102 Only one initializer for a default argument
103
104 Revision 1.35  2000/06/07 10:13:14  hristov
105 Delete only existent objects.
106
107 Revision 1.34  2000/05/18 10:45:38  fca
108 Delete Particle Factory properly
109
110 Revision 1.33  2000/05/16 13:10:40  fca
111 New method IsNewTrack and fix for a problem in Father-Daughter relations
112
113 Revision 1.32  2000/04/27 10:38:21  fca
114 Correct termination of Lego Run and introduce Lego getter in AliRun
115
116 Revision 1.31  2000/04/26 10:17:32  fca
117 Changes in Lego for G4 compatibility
118
119 Revision 1.30  2000/04/18 19:11:40  fca
120 Introduce variable Config.C function signature
121
122 Revision 1.29  2000/04/07 11:12:34  fca
123 G4 compatibility changes
124
125 Revision 1.28  2000/04/05 06:51:06  fca
126 Workaround for an HP compiler problem
127
128 Revision 1.27  2000/03/22 18:08:07  fca
129 Rationalisation of the virtual MC interfaces
130
131 Revision 1.26  2000/03/22 13:42:26  fca
132 SetGenerator does not replace an existing generator, ResetGenerator does
133
134 Revision 1.25  2000/02/23 16:25:22  fca
135 AliVMC and AliGeant3 classes introduced
136 ReadEuclid moved from AliRun to AliModule
137
138 Revision 1.24  2000/01/19 17:17:20  fca
139 Introducing a list of lists of hits -- more hits allowed for detector now
140
141 Revision 1.23  1999/12/03 11:14:31  fca
142 Fixing previous wrong checking
143
144 Revision 1.21  1999/11/25 10:40:08  fca
145 Fixing daughters information also in primary tracks
146
147 Revision 1.20  1999/10/04 18:08:49  fca
148 Adding protection against inconsistent Euclid files
149
150 Revision 1.19  1999/09/29 07:50:40  fca
151 Introduction of the Copyright and cvs Log
152
153 */
154
155 ///////////////////////////////////////////////////////////////////////////////
156 //                                                                           //
157 //  Control class for Alice C++                                              //
158 //  Only one single instance of this class exists.                           //
159 //  The object is created in main program aliroot                            //
160 //  and is pointed by the global gAlice.                                     //
161 //                                                                           //
162 //   -Supports the list of all Alice Detectors (fModules).                 //
163 //   -Supports the list of particles (fParticles).                           //
164 //   -Supports the Trees.                                                    //
165 //   -Supports the geometry.                                                 //
166 //   -Supports the event display.                                            //
167 //Begin_Html
168 /*
169 <img src="picts/AliRunClass.gif">
170 */
171 //End_Html
172 //Begin_Html
173 /*
174 <img src="picts/alirun.gif">
175 */
176 //End_Html
177 //                                                                           //
178 ///////////////////////////////////////////////////////////////////////////////
179
180 #include <stdlib.h>
181 #include <stdio.h>
182 #include <string.h>
183 #include <iostream.h>
184  
185 #include <TFile.h>
186 #include <TRandom.h>
187 #include <TBRIK.h> 
188 #include <TNode.h> 
189 #include <TCint.h> 
190 #include <TSystem.h>
191 #include <TObjectTable.h>
192 #include <TTree.h>
193 #include <TGeometry.h>
194 #include <TROOT.h>
195 #include <TBrowser.h>
196 #include <TFolder.h>
197
198 #include "TParticle.h"
199 #include "AliRun.h"
200 #include "AliDisplay.h"
201 #include "AliMC.h"
202 #include "AliLego.h"
203 #include "AliMagFC.h"
204 #include "AliMagFCM.h"
205 #include "AliMagFDM.h"
206 #include "AliHit.h"
207 #include "TRandom3.h"
208 #include "AliMCQA.h"
209 #include "AliGenerator.h"
210 #include "AliLegoGenerator.h"
211 #include "AliConfig.h"
212 #include "AliStack.h"
213 #include "AliGenEventHeader.h"
214 #include "AliHeader.h"
215
216 #include "AliDetector.h"
217
218 AliRun *gAlice;
219
220 ClassImp(AliRun)
221
222 //_____________________________________________________________________________
223 AliRun::AliRun()
224 {
225   //
226   // Default constructor for AliRun
227   //
228   fHeader    = 0;
229   fRun       = 0;
230   fEvent     = 0;
231   fStack     = 0;
232   fModules   = 0;
233   fGenerator = 0;
234   fTreeD     = 0;
235   fTreeH     = 0;
236   fTreeE     = 0;
237   fTreeR     = 0;
238   fTreeS     = 0;
239   fGeometry  = 0;
240   fDisplay   = 0;
241   fField     = 0;
242   fMC        = 0;
243   fNdets     = 0;
244   fImedia    = 0;
245   fTrRmax    = 1.e10;
246   fTrZmax    = 1.e10;
247   fInitDone  = kFALSE;
248   fLego      = 0;
249   fPDGDB     = 0;        //Particle factory object!
250   fHitLists  = 0;
251   fConfigFunction    = "\0";
252   fRandom = 0;
253   fMCQA = 0;
254   fTransParName = "\0";
255   fBaseFileName = ".\0";
256   fDebug        = 0;
257 }
258
259 //_____________________________________________________________________________
260 AliRun::AliRun(const char *name, const char *title)
261   : TNamed(name,title)
262 {
263   //
264   //  Constructor for the main processor.
265   //  Creates the geometry
266   //  Creates the list of Detectors.
267   //  Creates the list of particles.
268   //
269   Int_t i;
270   
271   gAlice     = this;
272   fTreeD     = 0;
273   fTreeH     = 0;
274   fTreeE     = 0;
275   fTreeR     = 0;
276   fTreeS     = 0;
277   fTrRmax    = 1.e10;
278   fTrZmax    = 1.e10;
279   fGenerator = 0;
280   fInitDone  = kFALSE;
281   fLego      = 0;
282   fField     = 0;
283   fConfigFunction    = "Config();";
284
285   // Set random number generator
286   gRandom = fRandom = new TRandom3();
287
288   if (gSystem->Getenv("CONFIG_SEED")) {
289      gRandom->SetSeed((UInt_t)atoi(gSystem->Getenv("CONFIG_SEED")));
290   }
291   
292   gROOT->GetListOfBrowsables()->Add(this,name);
293   //
294   // Particle stack
295   fStack = new AliStack(10000);
296   // create the support list for the various Detectors
297   fModules = new TObjArray(77);
298   //
299   // Create the TNode geometry for the event display
300   
301   BuildSimpleGeometry();
302   
303   fHeader    = new AliHeader();
304   fRun       = 0;
305   fEvent     = 0;
306   //
307   fDisplay = 0;
308   //
309   // Create default mag field
310   SetField();
311   //
312   fMC      = gMC;
313   //
314   // Prepare the tracking medium lists
315   fImedia = new TArrayI(1000);
316   for(i=0;i<1000;i++) (*fImedia)[i]=-99;
317   //
318   // Make particles
319   fPDGDB     = TDatabasePDG::Instance();        //Particle factory object!
320
321   AliConfig::Instance()->Add(fPDGDB); 
322   //
323   // Create HitLists list
324   fHitLists  = new TList();
325   //
326   SetTransPar();
327   fBaseFileName = ".\0";
328   //
329   fDebug        = 0;
330 }
331
332
333 //_____________________________________________________________________________
334 AliRun::~AliRun()
335 {
336   //
337   // Default AliRun destructor
338   //
339   delete fImedia;
340   delete fField;
341   delete fMC;
342   delete fGeometry;
343   delete fDisplay;
344   delete fGenerator;
345   delete fLego;
346   delete fTreeD;
347   delete fTreeH;
348   delete fTreeE;
349   delete fTreeR;
350   delete fTreeS;
351   if (fModules) {
352     fModules->Delete();
353     delete fModules;
354   }
355   delete fStack;
356   delete fHitLists;
357   delete fPDGDB;
358   delete fMCQA;
359   delete fHeader;
360 }
361
362 //_____________________________________________________________________________
363 void AliRun::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
364 {
365   //
366   //  Add a hit to detector id
367   //
368   TObjArray &dets = *fModules;
369   if(dets[id]) ((AliModule*) dets[id])->AddHit(track,vol,hits);
370 }
371
372 //_____________________________________________________________________________
373 void AliRun::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
374 {
375   //
376   // Add digit to detector id
377   //
378   TObjArray &dets = *fModules;
379   if(dets[id]) ((AliModule*) dets[id])->AddDigit(tracks,digits);
380 }
381
382 //_____________________________________________________________________________
383 void AliRun::Browse(TBrowser *b)
384 {
385   //
386   // Called when the item "Run" is clicked on the left pane
387   // of the Root browser.
388   // It displays the Root Trees and all detectors.
389   //
390   if(!fStack) fStack=fHeader->Stack();
391   TTree*  pTreeK = fStack->TreeK();
392     
393   if (pTreeK) b->Add(pTreeK,pTreeK->GetName());
394   if (fTreeH) b->Add(fTreeH,fTreeH->GetName());
395   if (fTreeD) b->Add(fTreeD,fTreeD->GetName());
396   if (fTreeE) b->Add(fTreeE,fTreeE->GetName());
397   if (fTreeR) b->Add(fTreeR,fTreeR->GetName());
398   if (fTreeS) b->Add(fTreeS,fTreeS->GetName());
399   
400   TIter next(fModules);
401   AliModule *detector;
402   while((detector = (AliModule*)next())) {
403     b->Add(detector,detector->GetName());
404   }
405   b->Add(fMCQA,"AliMCQA");
406 }
407
408 //_____________________________________________________________________________
409 void AliRun::Build()
410 {
411   //
412   // Initialize Alice geometry
413   // Dummy routine
414   //
415 }
416  
417 //_____________________________________________________________________________
418 void AliRun::BuildSimpleGeometry()
419 {
420   //
421   // Create a simple TNode geometry used by Root display engine
422   //
423   // Initialise geometry
424   //
425   fGeometry = new TGeometry("AliceGeom","Galice Geometry for Hits");
426   new TMaterial("void","Vacuum",0,0,0);  //Everything is void
427   TBRIK *brik = new TBRIK("S_alice","alice volume","void",2000,2000,3000);
428   brik->SetVisibility(0);
429   new TNode("alice","alice","S_alice");
430 }
431
432 //_____________________________________________________________________________
433 void AliRun::CleanDetectors()
434 {
435   //
436   // Clean Detectors at the end of event
437   //
438   TIter next(fModules);
439   AliModule *detector;
440   while((detector = (AliModule*)next())) {
441     detector->FinishEvent();
442   }
443 }
444
445 //_____________________________________________________________________________
446 Int_t AliRun::DistancetoPrimitive(Int_t, Int_t)
447 {
448   //
449   // Return the distance from the mouse to the AliRun object
450   // Dummy routine
451   //
452   return 9999;
453 }
454
455 //_____________________________________________________________________________
456 void AliRun::DumpPart (Int_t i) const
457 {
458   //
459   // Dumps particle i in the stack
460   //
461     fStack->DumpPart(i);
462 }
463
464 //_____________________________________________________________________________
465 void AliRun::DumpPStack () const
466 {
467   //
468   // Dumps the particle stack
469   //
470     fStack->DumpPStack();
471 }
472
473 //_____________________________________________________________________________
474 void  AliRun::SetField(AliMagF* magField)
475 {
476     // Set Magnetic Field Map
477     fField = magField;
478     fField->ReadField();
479 }
480
481 //_____________________________________________________________________________
482 void AliRun::SetField(Int_t type, Int_t version, Float_t scale,
483                       Float_t maxField, char* filename)
484 {
485   //
486   //  Set magnetic field parameters
487   //  type      Magnetic field transport flag 0=no field, 2=helix, 3=Runge Kutta
488   //  version   Magnetic field map version (only 1 active now)
489   //  scale     Scale factor for the magnetic field
490   //  maxField  Maximum value for the magnetic field
491
492   //
493   // --- Sanity check on mag field flags
494   if(fField) delete fField;
495   if(version==1) {
496     fField = new AliMagFC("Map1"," ",type,scale,maxField);
497   } else if(version<=2) {
498     fField = new AliMagFCM("Map2-3",filename,type,scale,maxField);
499     fField->ReadField();
500   } else if(version==3) {
501     fField = new AliMagFDM("Map4",filename,type,scale,maxField);
502     fField->ReadField();
503   } else {
504     Warning("SetField","Invalid map %d\n",version);
505   }
506 }
507  
508 //_____________________________________________________________________________
509 void AliRun::PreTrack()
510 {
511      TObjArray &dets = *fModules;
512      AliModule *module;
513
514      for(Int_t i=0; i<=fNdets; i++)
515        if((module = (AliModule*)dets[i]))
516          module->PreTrack();
517
518      fMCQA->PreTrack();
519 }
520
521 //_____________________________________________________________________________
522 void AliRun::PostTrack()
523 {
524      TObjArray &dets = *fModules;
525      AliModule *module;
526
527      for(Int_t i=0; i<=fNdets; i++)
528        if((module = (AliModule*)dets[i]))
529          module->PostTrack();
530 }
531
532 //_____________________________________________________________________________
533 void AliRun::FinishPrimary()
534 {
535   //
536   // Called  at the end of each primary track
537   //
538   
539   //  static Int_t count=0;
540   //  const Int_t times=10;
541   // This primary is finished, purify stack
542   fStack->PurifyKine();
543
544   TIter next(fModules);
545   AliModule *detector;
546   while((detector = (AliModule*)next())) {
547     detector->FinishPrimary();
548   }
549
550   // Write out hits if any
551   if (gAlice->TreeH()) {
552     gAlice->TreeH()->Fill();
553   }
554   
555   //
556   //  if(++count%times==1) gObjectTable->Print();
557 }
558
559 //_____________________________________________________________________________
560 void AliRun::BeginPrimary()
561 {
562   //
563   // Called  at the beginning of each primary track
564   //
565   
566   // Reset Hits info
567   gAlice->ResetHits();
568
569 }
570
571 //_____________________________________________________________________________
572 void AliRun::FinishEvent()
573 {
574   //
575   // Called at the end of the event.
576   //
577   
578   //
579   if(fLego) fLego->FinishEvent();
580
581   //Update the energy deposit tables
582   Int_t i;
583   for(i=0;i<fEventEnergy.GetSize();i++) {
584     fSummEnergy[i]+=fEventEnergy[i];
585     fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
586   }
587
588
589   
590   // Update Header information 
591
592   fHeader->SetNprimary(fStack->GetNprimary());
593   fHeader->SetNtrack(fStack->GetNtrack());  
594
595   
596   // Write out the kinematics
597   fStack->FinishEvent();
598    
599   // Write out the event Header information
600   if (fTreeE) {
601       fHeader->SetStack(fStack);
602       fTreeE->Fill();
603   }
604   
605   
606   // Write Tree headers
607   TTree* pTreeK = fStack->TreeK();
608   if (pTreeK) pTreeK->Write(0,TObject::kOverwrite);
609   if (fTreeH) fTreeH->Write(0,TObject::kOverwrite);
610   
611   ++fEvent;
612 }
613
614 //_____________________________________________________________________________
615 void AliRun::FinishRun()
616 {
617   //
618   // Called at the end of the run.
619   //
620
621   //
622   if(fLego) fLego->FinishRun();
623   
624   // Clean detector information
625   TIter next(fModules);
626   AliModule *detector;
627   while((detector = (AliModule*)next())) {
628     detector->FinishRun();
629   }
630   
631   //Output energy summary tables
632   EnergySummary();
633
634   TFile *file = fTreeE->GetCurrentFile();
635
636   file->cd();
637   
638   fTreeE->Write(0,TObject::kOverwrite);
639   
640   // Write AliRun info and all detectors parameters
641   Write(0,TObject::kOverwrite);
642   
643   // Clean tree information
644
645   fStack->FinishRun();
646   
647   if (fTreeH) {
648     delete fTreeH; fTreeH = 0;
649   }
650   if (fTreeD) {
651     delete fTreeD; fTreeD = 0;
652   }
653   if (fTreeR) {
654     delete fTreeR; fTreeR = 0;
655   }
656   if (fTreeE) {
657     delete fTreeE; fTreeE = 0;
658   }
659   if (fTreeS) {
660     delete fTreeS; fTreeS = 0;
661   }
662     
663   // Close output file
664   file->Write();
665 }
666
667 //_____________________________________________________________________________
668 void AliRun::FlagTrack(Int_t track)
669 {
670   // Delegate to stack
671   //
672     fStack->FlagTrack(track);
673 }
674  
675 //_____________________________________________________________________________
676 void AliRun::EnergySummary()
677 {
678   //
679   // Print summary of deposited energy
680   //
681
682   Int_t ndep=0;
683   Float_t edtot=0;
684   Float_t ed, ed2;
685   Int_t kn, i, left, j, id;
686   const Float_t kzero=0;
687   Int_t ievent=fHeader->GetEvent()+1;
688   //
689   // Energy loss information
690   if(ievent) {
691     printf("***************** Energy Loss Information per event (GEV) *****************\n");
692     for(kn=1;kn<fEventEnergy.GetSize();kn++) {
693       ed=fSummEnergy[kn];
694       if(ed>0) {
695         fEventEnergy[ndep]=kn;
696         if(ievent>1) {
697           ed=ed/ievent;
698           ed2=fSum2Energy[kn];
699           ed2=ed2/ievent;
700           ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
701         } else 
702           ed2=99;
703         fSummEnergy[ndep]=ed;
704         fSum2Energy[ndep]=TMath::Min((Float_t) 99.,TMath::Max(ed2,kzero));
705         edtot+=ed;
706         ndep++;
707       }
708     }
709     for(kn=0;kn<(ndep-1)/3+1;kn++) {
710       left=ndep-kn*3;
711       for(i=0;i<(3<left?3:left);i++) {
712         j=kn*3+i;
713         id=Int_t (fEventEnergy[j]+0.1);
714         printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
715       }
716       printf("\n");
717     }
718     //
719     // Relative energy loss in different detectors
720     printf("******************** Relative Energy Loss per event ********************\n");
721     printf("Total energy loss per event %10.3f GeV\n",edtot);
722     for(kn=0;kn<(ndep-1)/5+1;kn++) {
723       left=ndep-kn*5;
724       for(i=0;i<(5<left?5:left);i++) {
725         j=kn*5+i;
726         id=Int_t (fEventEnergy[j]+0.1);
727         printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
728       }
729       printf("\n");
730     }
731     for(kn=0;kn<75;kn++) printf("*"); 
732     printf("\n");
733   }
734   //
735   // Reset the TArray's
736   //  fEventEnergy.Set(0);
737   //  fSummEnergy.Set(0);
738   //  fSum2Energy.Set(0);
739 }
740
741 //_____________________________________________________________________________
742 AliModule *AliRun::GetModule(const char *name) const
743 {
744   //
745   // Return pointer to detector from name
746   //
747   return (AliModule*)fModules->FindObject(name);
748 }
749  
750 //_____________________________________________________________________________
751 AliDetector *AliRun::GetDetector(const char *name) const
752 {
753   //
754   // Return pointer to detector from name
755   //
756   return (AliDetector*)fModules->FindObject(name);
757 }
758  
759 //_____________________________________________________________________________
760 Int_t AliRun::GetModuleID(const char *name) const
761 {
762   //
763   // Return galice internal detector identifier from name
764   //
765   Int_t i=-1;
766   TObject *mod=fModules->FindObject(name);
767   if(mod) i=fModules->IndexOf(mod);
768   return i;
769 }
770  
771 //_____________________________________________________________________________
772 Int_t AliRun::GetEvent(Int_t event)
773 {
774   //
775   // Connect the Trees Kinematics and Hits for event # event
776   // Set branch addresses
777   //
778
779   // Reset existing structures
780   ResetHits();
781   ResetDigits();
782   ResetSDigits();
783   
784   // Delete Trees already connected
785   if (fTreeH) delete fTreeH;
786   if (fTreeD) delete fTreeD;
787   if (fTreeR) delete fTreeR;
788   if (fTreeS) delete fTreeS;
789
790  // Create the particle stack
791   if (fHeader) delete fHeader; 
792   fHeader = 0;
793   
794   // Get header from file
795   if(fTreeE) {
796       fTreeE->SetBranchAddress("Header", &fHeader);
797       fTreeE->GetEntry(event);
798   }  
799   else 
800       Error("GetEvent","Cannot find Header Tree (TE)\n");
801
802   // Get the stack from the header
803   if (fStack) delete fStack;
804   fStack = fHeader->Stack();
805   fStack->GetEvent(event);
806   //
807   TFile *file = fTreeE->GetCurrentFile();
808   char treeName[20];
809   file->cd();
810
811   // Get Hits Tree header from file
812   sprintf(treeName,"TreeH%d",event);
813   fTreeH = (TTree*)gDirectory->Get(treeName);
814   if (!fTreeH) {
815     Error("GetEvent","cannot find Hits Tree for event:%d\n",event);
816   }
817   
818   // Get Digits Tree header from file
819   sprintf(treeName,"TreeD%d",event);
820   fTreeD = (TTree*)gDirectory->Get(treeName);
821   if (!fTreeD) {
822     // Warning("GetEvent","cannot find Digits Tree for event:%d\n",event);
823   }
824
825   file->cd();
826
827   // Get SDigits Tree header from file
828   sprintf(treeName,"TreeS%d",event);
829   fTreeS = (TTree*)gDirectory->Get(treeName);
830   if (!fTreeS) {
831     // Warning("GetEvent","cannot find SDigits Tree for event:%d\n",event);
832   }
833
834   file->cd();
835   
836   // Get Reconstruct Tree header from file
837   sprintf(treeName,"TreeR%d",event);
838   fTreeR = (TTree*)gDirectory->Get(treeName);
839   if (!fTreeR) {
840     //    printf("WARNING: cannot find Reconstructed Tree for event:%d\n",event);
841   }
842
843   file->cd();
844  
845   // Set Trees branch addresses
846   TIter next(fModules);
847   AliModule *detector;
848   while((detector = (AliModule*)next())) {
849     detector->SetTreeAddress();
850   }
851
852   return fStack->GetNtrack();
853 }
854
855 //_____________________________________________________________________________
856 TGeometry *AliRun::GetGeometry()
857 {
858   //
859   // Import Alice geometry from current file
860   // Return pointer to geometry object
861   //
862   if (!fGeometry) fGeometry = (TGeometry*)gDirectory->Get("AliceGeom");
863   //
864   // Unlink and relink nodes in detectors
865   // This is bad and there must be a better way...
866   //
867   
868   TIter next(fModules);
869   AliModule *detector;
870   while((detector = (AliModule*)next())) {
871     TList *dnodes=detector->Nodes();
872     Int_t j;
873     TNode *node, *node1;
874     for ( j=0; j<dnodes->GetSize(); j++) {
875       node = (TNode*) dnodes->At(j);
876       node1 = fGeometry->GetNode(node->GetName());
877       dnodes->Remove(node);
878       dnodes->AddAt(node1,j);
879     }
880   }
881   return fGeometry;
882 }
883
884 //_____________________________________________________________________________
885 void AliRun::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom,
886                           Float_t &e, Float_t *vpos, Float_t *polar,
887                           Float_t &tof)
888 {
889   // Delegate to stack
890   //
891     fStack->GetNextTrack(mtrack, ipart, pmom, e, vpos, polar, tof);  
892 }
893
894 //_____________________________________________________________________________
895 Int_t AliRun::GetPrimary(Int_t track) const
896 {
897   //
898   // return number of primary that has generated track
899   //
900     return fStack->GetPrimary(track);
901 }
902  
903 //_____________________________________________________________________________
904 void AliRun::InitMC(const char *setup)
905 {
906   //
907   // Initialize the Alice setup
908   //
909
910   if(fInitDone) {
911     Warning("Init","Cannot initialise AliRun twice!\n");
912     return;
913   }
914     
915   gROOT->LoadMacro(setup);
916   gInterpreter->ProcessLine(fConfigFunction.Data());
917
918
919   gMC->DefineParticles();  //Create standard MC particles
920
921   TObject *objfirst, *objlast;
922
923   fNdets = fModules->GetLast()+1;
924
925   //
926   //=================Create Materials and geometry
927   gMC->Init();
928
929   // Added also after in case of interactive initialisation of modules
930   fNdets = fModules->GetLast()+1;
931
932    TIter next(fModules);
933    AliModule *detector;
934    while((detector = (AliModule*)next())) {
935       detector->SetTreeAddress();
936       objlast = gDirectory->GetList()->Last();
937       
938       // Add Detector histograms in Detector list of histograms
939       if (objlast) objfirst = gDirectory->GetList()->After(objlast);
940       else         objfirst = gDirectory->GetList()->First();
941       while (objfirst) {
942         detector->Histograms()->Add(objfirst);
943         objfirst = gDirectory->GetList()->After(objfirst);
944       }
945    }
946    ReadTransPar(); //Read the cuts for all materials
947    
948    MediaTable(); //Build the special IMEDIA table
949    
950    //Initialise geometry deposition table
951    fEventEnergy.Set(gMC->NofVolumes()+1);
952    fSummEnergy.Set(gMC->NofVolumes()+1);
953    fSum2Energy.Set(gMC->NofVolumes()+1);
954    
955    //Compute cross-sections
956    gMC->BuildPhysics();
957    
958    //Write Geometry object to current file.
959    fGeometry->Write();
960    
961    fInitDone = kTRUE;
962
963    fMCQA = new AliMCQA(fNdets);
964
965    //
966    // Save stuff at the beginning of the file to avoid file corruption
967    Write();
968 }
969
970 //_____________________________________________________________________________
971 void AliRun::MediaTable()
972 {
973   //
974   // Built media table to get from the media number to
975   // the detector id
976   //
977   Int_t kz, nz, idt, lz, i, k, ind;
978   //  Int_t ibeg;
979   TObjArray &dets = *gAlice->Detectors();
980   AliModule *det;
981   //
982   // For all detectors
983   for (kz=0;kz<fNdets;kz++) {
984     // If detector is defined
985     if((det=(AliModule*) dets[kz])) {
986         TArrayI &idtmed = *(det->GetIdtmed()); 
987         for(nz=0;nz<100;nz++) {
988         // Find max and min material number
989         if((idt=idtmed[nz])) {
990           det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
991           det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
992         }
993       }
994       if(det->LoMedium() > det->HiMedium()) {
995         det->LoMedium() = 0;
996         det->HiMedium() = 0;
997       } else {
998         if(det->HiMedium() > fImedia->GetSize()) {
999           Error("MediaTable","Increase fImedia from %d to %d",
1000                 fImedia->GetSize(),det->HiMedium());
1001           return;
1002         }
1003         // Tag all materials in rage as belonging to detector kz
1004         for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
1005           (*fImedia)[lz]=kz;
1006         }
1007       }
1008     }
1009   }
1010   //
1011   // Print summary table
1012   printf(" Traking media ranges:\n");
1013   for(i=0;i<(fNdets-1)/6+1;i++) {
1014     for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
1015       ind=i*6+k;
1016       det=(AliModule*)dets[ind];
1017       if(det)
1018         printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
1019                det->HiMedium());
1020       else
1021         printf(" %6s: %3d -> %3d;","NULL",0,0);
1022     }
1023     printf("\n");
1024   }
1025 }
1026
1027 //____________________________________________________________________________
1028 void AliRun::SetGenerator(AliGenerator *generator)
1029 {
1030   //
1031   // Load the event generator
1032   //
1033   if(!fGenerator) fGenerator = generator;
1034 }
1035
1036 //____________________________________________________________________________
1037 void AliRun::ResetGenerator(AliGenerator *generator)
1038 {
1039   //
1040   // Load the event generator
1041   //
1042   if(fGenerator)
1043     if(generator)
1044       Warning("ResetGenerator","Replacing generator %s with %s\n",
1045               fGenerator->GetName(),generator->GetName());
1046     else
1047       Warning("ResetGenerator","Replacing generator %s with NULL\n",
1048               fGenerator->GetName());
1049   fGenerator = generator;
1050 }
1051
1052 //____________________________________________________________________________
1053 void AliRun::SetTransPar(char *filename)
1054 {
1055   fTransParName = filename;
1056 }
1057
1058 //____________________________________________________________________________
1059 void AliRun::SetBaseFile(char *filename)
1060 {
1061   fBaseFileName = filename;
1062 }
1063
1064 //____________________________________________________________________________
1065 void AliRun::ReadTransPar()
1066 {
1067   //
1068   // Read filename to set the transport parameters
1069   //
1070
1071
1072   const Int_t kncuts=10;
1073   const Int_t knflags=11;
1074   const Int_t knpars=kncuts+knflags;
1075   const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
1076                                "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
1077                                "BREM","COMP","DCAY","DRAY","HADR","LOSS",
1078                                "MULS","PAIR","PHOT","RAYL"};
1079   char line[256];
1080   char detName[7];
1081   char* filtmp;
1082   Float_t cut[kncuts];
1083   Int_t flag[knflags];
1084   Int_t i, itmed, iret, ktmed, kz;
1085   FILE *lun;
1086   //
1087   // See whether the file is there
1088   filtmp=gSystem->ExpandPathName(fTransParName.Data());
1089   lun=fopen(filtmp,"r");
1090   delete [] filtmp;
1091   if(!lun) {
1092     Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
1093     return;
1094   }
1095   //
1096   if(fDebug) {
1097     printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1098     printf(" *%59s\n","*");
1099     printf(" *       Please check carefully what you are doing!%10s\n","*");
1100     printf(" *%59s\n","*");
1101   }
1102   //
1103   while(1) {
1104     // Initialise cuts and flags
1105     for(i=0;i<kncuts;i++) cut[i]=-99;
1106     for(i=0;i<knflags;i++) flag[i]=-99;
1107     itmed=0;
1108     for(i=0;i<256;i++) line[i]='\0';
1109     // Read up to the end of line excluded
1110     iret=fscanf(lun,"%[^\n]",line);
1111     if(iret<0) {
1112       //End of file
1113       fclose(lun);
1114       if(fDebug){
1115         printf(" *%59s\n","*");
1116         printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1117       }
1118       return;
1119     }
1120     // Read the end of line
1121     fscanf(lun,"%*c");
1122     if(!iret) continue;
1123     if(line[0]=='*') continue;
1124     // Read the numbers
1125     iret=sscanf(line,"%s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d",
1126                 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
1127                 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
1128                 &flag[8],&flag[9],&flag[10]);
1129     if(!iret) continue;
1130     if(iret<0) {
1131       //reading error
1132       Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
1133       continue;
1134     }
1135     // Check that the module exist
1136     AliModule *mod = GetModule(detName);
1137     if(mod) {
1138       // Get the array of media numbers
1139       TArrayI &idtmed = *mod->GetIdtmed();
1140       // Check that the tracking medium code is valid
1141       if(0<=itmed && itmed < 100) {
1142         ktmed=idtmed[itmed];
1143         if(!ktmed) {
1144           Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
1145           continue;
1146         }
1147         // Set energy thresholds
1148         for(kz=0;kz<kncuts;kz++) {
1149           if(cut[kz]>=0) {
1150             if(fDebug) printf(" *  %-6s set to %10.3E for tracking medium code %4d for %s\n",
1151                    kpars[kz],cut[kz],itmed,mod->GetName());
1152             gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
1153           }
1154         }
1155         // Set transport mechanisms
1156         for(kz=0;kz<knflags;kz++) {
1157           if(flag[kz]>=0) {
1158             if(fDebug) printf(" *  %-6s set to %10d for tracking medium code %4d for %s\n",
1159                    kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
1160             gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
1161           }
1162         }
1163       } else {
1164         Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
1165         continue;
1166       }
1167     } else {
1168       if(fDebug) printf("%s::ReadTransParModule: %s not present\n",ClassName(),detName);
1169       continue;
1170     }
1171   }
1172 }
1173
1174
1175 //_____________________________________________________________________________
1176 void AliRun::MakeTree(Option_t *option, const char *file)
1177 {
1178   //
1179   //  Create the ROOT trees
1180   //  Loop on all detectors to create the Root branch (if any)
1181   //
1182
1183   char hname[30];
1184   //
1185   // Analyse options
1186   const char *oK = strstr(option,"K");
1187   const char *oH = strstr(option,"H");
1188   const char *oE = strstr(option,"E");
1189   const char *oD = strstr(option,"D");
1190   const char *oR = strstr(option,"R");
1191   const char *oS = strstr(option,"S");
1192   //
1193
1194   TDirectory *cwd = gDirectory;
1195
1196   TBranch *branch = 0;
1197   
1198   if (oK) fStack->MakeTree(fEvent, file);
1199
1200   if (oE && !fTreeE) {
1201     fTreeE = new TTree("TE","Header");
1202     branch = fTreeE->Branch("Header", "AliHeader", &fHeader, 4000, 0);         
1203     branch->SetAutoDelete(kFALSE);           
1204     TFolder *folder = (TFolder *)gROOT->FindObjectAny("/Folders/RunMC/Event/Header");
1205     if (folder) folder->Add(fHeader);
1206 //    branch = fTreeE->Branch("Stack","AliStack", &fStack, 4000, 0);
1207 //    branch->SetAutoDelete(kFALSE);         
1208 //    if (folder) folder->Add(fStack);
1209     fTreeE->Write(0,TObject::kOverwrite);
1210   }
1211   
1212   if (file && branch) {
1213     char * outFile = new char[strlen(gAlice->GetBaseFile())+strlen(file)+2];
1214     sprintf(outFile,"%s/%s",GetBaseFile(),file);
1215     branch->SetFile(outFile);
1216     TIter next( branch->GetListOfBranches());
1217     while ((branch=(TBranch*)next())) {
1218        branch->SetFile(outFile);
1219     } 
1220     if (GetDebug()>1)
1221         printf("* MakeBranch * Diverting Branch %s to file %s\n", branch->GetName(),file);
1222     cwd->cd();
1223     delete outFile;
1224   }
1225   
1226   if (oH && !fTreeH) {
1227     sprintf(hname,"TreeH%d",fEvent);
1228     fTreeH = new TTree(hname,"Hits");
1229     fTreeH->SetAutoSave(1000000000); //no autosave
1230     fTreeH->Write(0,TObject::kOverwrite);
1231   }
1232   if (oD && !fTreeD) {
1233     sprintf(hname,"TreeD%d",fEvent);
1234     fTreeD = new TTree(hname,"Digits");
1235     fTreeD->Write(0,TObject::kOverwrite);
1236   }
1237   if (oS && !fTreeS) {
1238     sprintf(hname,"TreeS%d",fEvent);
1239     fTreeS = new TTree(hname,"SDigits");
1240     fTreeS->Write(0,TObject::kOverwrite);
1241   }
1242   if (oR && !fTreeR) {
1243     sprintf(hname,"TreeR%d",fEvent);
1244     fTreeR = new TTree(hname,"Reconstruction");
1245     fTreeR->Write(0,TObject::kOverwrite);
1246   }
1247
1248   //
1249   // Create a branch for hits/digits for each detector
1250   // Each branch is a TClonesArray. Each data member of the Hits classes
1251   // will be in turn a subbranch of the detector master branch
1252   TIter next(fModules);
1253   AliModule *detector;
1254   while((detector = (AliModule*)next())) {
1255      if (oH) detector->MakeBranch(option,file);
1256   }
1257 }
1258
1259 //_____________________________________________________________________________
1260 inline TParticle* AliRun::Particle(Int_t i)
1261 {
1262     return fStack->Particle(i);
1263 }
1264
1265 //_____________________________________________________________________________
1266 void AliRun::BeginEvent()
1267 {
1268     // Clean-up previous event
1269     // Energy scores  
1270     fEventEnergy.Reset();  
1271     // Clean detector information
1272     CleanDetectors();
1273     // Reset stack info
1274     fStack->Reset();
1275
1276  
1277   //
1278   //  Reset all Detectors & kinematics & trees
1279   //
1280   char hname[30];
1281   //
1282   // Initialise event header
1283   fHeader->Reset(fRun,fEvent);
1284   //
1285   fStack->BeginEvent(fEvent);
1286
1287   //
1288   if(fLego) {
1289     fLego->BeginEvent();
1290     return;
1291   }
1292
1293   //
1294
1295   ResetHits();
1296   ResetDigits();
1297   ResetSDigits();
1298
1299
1300   if(fTreeH) {
1301     fTreeH->Reset();
1302     sprintf(hname,"TreeH%d",fEvent);
1303     fTreeH->SetName(hname);
1304   }
1305   if(fTreeD) {
1306     fTreeD->Reset();
1307     sprintf(hname,"TreeD%d",fEvent);
1308     fTreeD->SetName(hname);
1309     fTreeD->Write(0,TObject::kOverwrite);
1310   }
1311   if(fTreeS) {
1312     fTreeS->Reset();
1313     sprintf(hname,"TreeS%d",fEvent);
1314     fTreeS->SetName(hname);
1315     fTreeS->Write(0,TObject::kOverwrite);
1316   }
1317   if(fTreeR) {
1318     fTreeR->Reset();
1319     sprintf(hname,"TreeR%d",fEvent);
1320     fTreeR->SetName(hname);
1321     fTreeR->Write(0,TObject::kOverwrite);
1322   }
1323 }
1324 //_____________________________________________________________________________
1325 void AliRun::ResetDigits()
1326 {
1327   //
1328   //  Reset all Detectors digits
1329   //
1330   TIter next(fModules);
1331   AliModule *detector;
1332   while((detector = (AliModule*)next())) {
1333      detector->ResetDigits();
1334   }
1335 }
1336
1337 //_____________________________________________________________________________
1338 void AliRun::ResetSDigits()
1339 {
1340   //
1341   //  Reset all Detectors digits
1342   //
1343   TIter next(fModules);
1344   AliModule *detector;
1345   while((detector = (AliModule*)next())) {
1346      detector->ResetSDigits();
1347   }
1348 }
1349
1350 //_____________________________________________________________________________
1351 void AliRun::ResetHits()
1352 {
1353   //
1354   //  Reset all Detectors hits
1355   //
1356   TIter next(fModules);
1357   AliModule *detector;
1358   while((detector = (AliModule*)next())) {
1359      detector->ResetHits();
1360   }
1361 }
1362
1363 //_____________________________________________________________________________
1364 void AliRun::ResetPoints()
1365 {
1366   //
1367   // Reset all Detectors points
1368   //
1369   TIter next(fModules);
1370   AliModule *detector;
1371   while((detector = (AliModule*)next())) {
1372      detector->ResetPoints();
1373   }
1374 }
1375
1376 //_____________________________________________________________________________
1377 void AliRun::RunMC(Int_t nevent, const char *setup)
1378 {
1379   //
1380   // Main function to be called to process a galice run
1381   // example
1382   //    Root > gAlice.Run(); 
1383   // a positive number of events will cause the finish routine
1384   // to be called
1385   //
1386
1387   // check if initialisation has been done
1388   if (!fInitDone) InitMC(setup);
1389   
1390   // Create the Root Tree with one branch per detector
1391
1392    MakeTree("ESDR");
1393
1394   if (gSystem->Getenv("CONFIG_SPLIT_FILE")) {
1395      MakeTree("K","Kine.root");
1396      MakeTree("H","Hits.root");
1397   } else {
1398      MakeTree("KH");
1399   }
1400
1401   gMC->ProcessRun(nevent);
1402
1403   // End of this run, close files
1404   if(nevent>0) FinishRun();
1405 }
1406
1407 //_____________________________________________________________________________
1408 void AliRun::RunReco(const char *selected)
1409 {
1410   //
1411   // Main function to be called to reconstruct Alice event
1412   //  
1413    for (Int_t nevent=0; nevent<gAlice->TreeE()->GetEntries(); nevent++) {
1414      GetEvent(nevent);
1415      // MakeTree("R");
1416      Digits2Reco(selected);
1417    }
1418 }
1419
1420 //_____________________________________________________________________________
1421
1422 void AliRun::Hits2Digits(const char *selected)
1423 {
1424
1425    // Convert Hits to sumable digits
1426    // 
1427    for (Int_t nevent=0; nevent<gAlice->TreeE()->GetEntries(); nevent++) {
1428      GetEvent(nevent);
1429     // MakeTree("D");
1430      Hits2SDigits(selected);
1431      SDigits2Digits(selected);
1432    }  
1433 }
1434
1435
1436 //_____________________________________________________________________________
1437
1438 void AliRun::Tree2Tree(Option_t *option, const char *selected)
1439 {
1440   //
1441   // Function to transform the content of
1442   //  
1443   // - TreeH to TreeS (option "S")
1444   // - TreeS to TreeD (option "D")
1445   // - TreeD to TreeR (option "R")
1446   // 
1447   // If multiple options are specified ("SDR"), transformation will be done in sequence for
1448   // selected detector and for all detectors if none is selected (detector string 
1449   // can contain blank separated list of detector names). 
1450
1451
1452    const char *oS = strstr(option,"S");
1453    const char *oD = strstr(option,"D");
1454    const char *oR = strstr(option,"R");
1455    
1456    TObjArray *detectors = Detectors();
1457
1458    TIter next(detectors);
1459
1460    AliDetector *detector = 0;
1461
1462    TDirectory *cwd = gDirectory;
1463
1464    char outFile[32];
1465    
1466    while((detector = (AliDetector*)next())) {
1467      if (selected) 
1468        if (strcmp(detector->GetName(),selected)) continue;
1469      if (detector->IsActive()){ 
1470        if (gSystem->Getenv("CONFIG_SPLIT_FILE")) {          
1471           if (oS) {
1472             sprintf(outFile,"SDigits.%s.root",detector->GetName());
1473             detector->MakeBranch("S",outFile);
1474           }    
1475           if (oD) {
1476             sprintf(outFile,"Digits.%s.root",detector->GetName());
1477             detector->MakeBranch("D",outFile);
1478           }    
1479           if (oR) {
1480             sprintf(outFile,"Reco.%s.root",detector->GetName());
1481             detector->MakeBranch("R",outFile);
1482           }    
1483        } else {
1484           detector->MakeBranch(option);
1485        }
1486        
1487        cwd->cd(); 
1488
1489        if (oS) {
1490           cout << "Hits2SDigits: Processing " << detector->GetName() << "..." << endl;
1491           detector->Hits2SDigits(); 
1492        }  
1493        if (oD) {
1494           cout << "SDigits2Digits: Processing " << detector->GetName() << "..." << endl;
1495           detector->SDigits2Digits();
1496        } 
1497        if (oR) {
1498           cout << "Digits2Reco: Processing " << detector->GetName() << "..." << endl;
1499           detector->Digits2Reco(); 
1500        }
1501        
1502        cwd->cd();        
1503      }   
1504    }
1505 }
1506
1507
1508 //_____________________________________________________________________________
1509 void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min,
1510                      Float_t c1max,Int_t nc2,Float_t c2min,Float_t c2max,
1511                      Float_t rmin,Float_t rmax,Float_t zmax, AliLegoGenerator* gener)
1512 {
1513   //
1514   // Generates lego plots of:
1515   //    - radiation length map phi vs theta
1516   //    - radiation length map phi vs eta
1517   //    - interaction length map
1518   //    - g/cm2 length map
1519   //
1520   //  ntheta    bins in theta, eta
1521   //  themin    minimum angle in theta (degrees)
1522   //  themax    maximum angle in theta (degrees)
1523   //  nphi      bins in phi
1524   //  phimin    minimum angle in phi (degrees)
1525   //  phimax    maximum angle in phi (degrees)
1526   //  rmin      minimum radius
1527   //  rmax      maximum radius
1528   //  
1529   //
1530   //  The number of events generated = ntheta*nphi
1531   //  run input parameters in macro setup (default="Config.C")
1532   //
1533   //  Use macro "lego.C" to visualize the 3 lego plots in spherical coordinates
1534   //Begin_Html
1535   /*
1536     <img src="picts/AliRunLego1.gif">
1537   */
1538   //End_Html
1539   //Begin_Html
1540   /*
1541     <img src="picts/AliRunLego2.gif">
1542   */
1543   //End_Html
1544   //Begin_Html
1545   /*
1546     <img src="picts/AliRunLego3.gif">
1547   */
1548   //End_Html
1549   //
1550
1551   // check if initialisation has been done
1552   if (!fInitDone) InitMC(setup);
1553   //Save current generator
1554   AliGenerator *gen=Generator();
1555
1556   // Set new generator
1557   if (!gener) gener  = new AliLegoGenerator();
1558   ResetGenerator(gener);
1559   //
1560   // Configure Generator
1561   gener->SetRadiusRange(rmin, rmax);
1562   gener->SetZMax(zmax);
1563   gener->SetCoor1Range(nc1, c1min, c1max);
1564   gener->SetCoor2Range(nc2, c2min, c2max);
1565   
1566   
1567   //Create Lego object  
1568   fLego = new AliLego("lego",gener);
1569
1570   //Prepare MC for Lego Run
1571   gMC->InitLego();
1572   
1573   //Run Lego Object
1574
1575   gMC->ProcessRun(nc1*nc2+1);
1576   
1577   // Create only the Root event Tree
1578   MakeTree("E");
1579   
1580   // End of this run, close files
1581   FinishRun();
1582   // Restore current generator
1583   ResetGenerator(gen);
1584   // Delete Lego Object
1585   delete fLego; fLego=0;
1586 }
1587
1588 //_____________________________________________________________________________
1589 void AliRun::SetConfigFunction(const char * config) 
1590 {
1591   //
1592   // Set the signature of the function contained in Config.C to configure
1593   // the run
1594   //
1595   fConfigFunction=config;
1596 }
1597
1598 //_____________________________________________________________________________
1599 void AliRun::SetCurrentTrack(Int_t track)
1600
1601   //
1602   // Set current track number
1603   //
1604     fStack->SetCurrentTrack(track); 
1605 }
1606  
1607 //_____________________________________________________________________________
1608 void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
1609                       Float_t *vpos, Float_t *polar, Float_t tof,
1610                       AliMCProcess mech, Int_t &ntr, Float_t weight)
1611
1612 // Delegate to stack
1613 //
1614     fStack->SetTrack(done, parent, pdg, pmom, vpos, polar, tof,
1615                      mech, ntr, weight);
1616 }
1617
1618 //_____________________________________________________________________________
1619 void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg,
1620                       Double_t px, Double_t py, Double_t pz, Double_t e,
1621                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
1622                       Double_t polx, Double_t poly, Double_t polz,
1623                       AliMCProcess mech, Int_t &ntr, Float_t weight)
1624
1625   // Delegate to stack
1626   //
1627     fStack->SetTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
1628                    polx, poly, polz, mech, ntr, weight);
1629     
1630 }
1631
1632 //_____________________________________________________________________________
1633 void AliRun::SetHighWaterMark(const Int_t nt)
1634 {
1635     //
1636     // Set high water mark for last track in event
1637     fStack->SetHighWaterMark(nt);
1638 }
1639
1640 //_____________________________________________________________________________
1641 void AliRun::KeepTrack(const Int_t track)
1642
1643   //
1644   // Delegate to stack
1645   //
1646     fStack->KeepTrack(track);
1647 }
1648  
1649 //_____________________________________________________________________________
1650 void AliRun::StepManager(Int_t id) 
1651 {
1652   //
1653   // Called at every step during transport
1654   //
1655
1656   //
1657   // --- If lego option, do it and leave 
1658   if (fLego)
1659     fLego->StepManager();
1660   else {
1661     Int_t copy;
1662     //Update energy deposition tables
1663     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
1664   
1665     //Call the appropriate stepping routine;
1666     AliModule *det = (AliModule*)fModules->At(id);
1667     if(det) {
1668       fMCQA->StepManager(id);
1669       det->StepManager();
1670     }
1671   }
1672 }
1673
1674 //_____________________________________________________________________________
1675 void AliRun::Streamer(TBuffer &R__b)
1676 {
1677    // Stream an object of class AliRun.
1678
1679    if (R__b.IsReading()) {
1680       if (!gAlice) gAlice = this;
1681
1682       AliRun::Class()->ReadBuffer(R__b, this);
1683       //
1684       gROOT->GetListOfBrowsables()->Add(this,"Run");
1685
1686       fTreeE = (TTree*)gDirectory->Get("TE");
1687       if (fTreeE) {
1688           fTreeE->SetBranchAddress("Header", &fHeader);
1689       }
1690       
1691       else    Error("Streamer","cannot find Header Tree\n");
1692       fTreeE->GetEntry(0);
1693
1694       gRandom = fRandom;
1695    } else {
1696       AliRun::Class()->WriteBuffer(R__b, this);
1697    }
1698 }
1699
1700
1701 //___________________________________________________________________________
1702 Int_t AliRun::CurrentTrack() const {
1703   //
1704   // Returns current track
1705   //
1706   return fStack->CurrentTrack();
1707 }
1708
1709 //___________________________________________________________________________
1710 Int_t AliRun::GetNtrack() const {
1711   //
1712   // Returns number of tracks in stack
1713   //
1714   return fStack->GetNtrack();
1715 }
1716
1717 //___________________________________________________________________________
1718 TObjArray* AliRun::Particles() {
1719   //
1720   // Returns pointer to Particles array
1721   //
1722   return fStack->Particles();
1723 }
1724
1725 //___________________________________________________________________________
1726 TTree* AliRun::TreeK() {
1727   //
1728   // Returns pointer to the TreeK array
1729   //
1730   return fStack->TreeK();
1731 }