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