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