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