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