e543370b762707b36e20ef8239f79e0d1b0d0335
[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   
433   if(fLego) 
434    {
435     if (GetDebug()) Info("FinishRun"," Finish Lego");
436     fRunLoader->CdGAFile();
437     fLego->FinishRun();
438    }
439   
440   // Clean detector information
441   TIter next(fModules);
442   AliModule *detector;
443   while((detector = dynamic_cast<AliModule*>(next()))) {
444     if (GetDebug()) Info("FinishRun"," %s->FinishRun()",detector->GetName());
445     detector->FinishRun();
446   }
447   
448   //Output energy summary tables
449   if (GetDebug()) Info("FinishRun"," EnergySummary()");
450   EnergySummary();
451   
452   if (GetDebug()) Info("FinishRun"," fRunLoader->WriteHeader(OVERWRITE)");
453   fRunLoader->WriteHeader("OVERWRITE");
454
455   // Write AliRun info and all detectors parameters
456   fRunLoader->CdGAFile();
457   Write(0,TObject::kOverwrite);//write AliRun
458   fRunLoader->Write(0,TObject::kOverwrite);//write RunLoader itself
459   
460   // Clean tree information
461   if (GetDebug()) Info("FinishRun"," fRunLoader->Stack()->FinishRun()");
462   fRunLoader->Stack()->FinishRun();
463
464   // Clean detector information
465   if (GetDebug()) Info("FinishRun"," fGenerator->FinishRun()");
466   fGenerator->FinishRun();
467 }
468
469 //_______________________________________________________________________
470 void AliRun::FlagTrack(Int_t track)
471 {
472   // Delegate to stack
473   //
474     fRunLoader->Stack()->FlagTrack(track);
475 }
476  
477 //_______________________________________________________________________
478 void AliRun::EnergySummary()
479 {
480   //
481   // Print summary of deposited energy
482   //
483
484   Int_t ndep=0;
485   Float_t edtot=0;
486   Float_t ed, ed2;
487   Int_t kn, i, left, j, id;
488   const Float_t kzero=0;
489   Int_t ievent=fRunLoader->GetHeader()->GetEvent()+1;
490   //
491   // Energy loss information
492   if(ievent) {
493     printf("***************** Energy Loss Information per event (GEV) *****************\n");
494     for(kn=1;kn<fEventEnergy.GetSize();kn++) {
495       ed=fSummEnergy[kn];
496       if(ed>0) {
497         fEventEnergy[ndep]=kn;
498         if(ievent>1) {
499           ed=ed/ievent;
500           ed2=fSum2Energy[kn];
501           ed2=ed2/ievent;
502           ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
503         } else 
504           ed2=99;
505         fSummEnergy[ndep]=ed;
506         fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
507         edtot+=ed;
508         ndep++;
509       }
510     }
511     for(kn=0;kn<(ndep-1)/3+1;kn++) {
512       left=ndep-kn*3;
513       for(i=0;i<(3<left?3:left);i++) {
514         j=kn*3+i;
515         id=Int_t (fEventEnergy[j]+0.1);
516         printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
517       }
518       printf("\n");
519     }
520     //
521     // Relative energy loss in different detectors
522     printf("******************** Relative Energy Loss per event ********************\n");
523     printf("Total energy loss per event %10.3f GeV\n",edtot);
524     for(kn=0;kn<(ndep-1)/5+1;kn++) {
525       left=ndep-kn*5;
526       for(i=0;i<(5<left?5:left);i++) {
527         j=kn*5+i;
528         id=Int_t (fEventEnergy[j]+0.1);
529         printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
530       }
531       printf("\n");
532     }
533     for(kn=0;kn<75;kn++) printf("*"); 
534     printf("\n");
535   }
536   //
537   // Reset the TArray's
538   //  fEventEnergy.Set(0);
539   //  fSummEnergy.Set(0);
540   //  fSum2Energy.Set(0);
541 }
542
543 //_______________________________________________________________________
544 void AliRun::Announce() const
545 {
546   //
547   // Announce the current version of AliRoot
548   //
549   printf("%70s",
550          "****************************************************************\n");
551   printf("%6s","*");printf("%64s","*\n");
552
553   printf("%6s","*");
554   printf("    You are running AliRoot version NewIO\n");
555
556   printf("%6s","*");
557   printf("    The cvs tag for the current program is $Name$\n");
558
559   printf("%6s","*");printf("%64s","*\n");
560   printf("%70s",
561          "****************************************************************\n");
562 }
563
564 //_______________________________________________________________________
565 AliModule *AliRun::GetModule(const char *name) const
566 {
567   //
568   // Return pointer to detector from name
569   //
570   return dynamic_cast<AliModule*>(fModules->FindObject(name));
571 }
572  
573 //_______________________________________________________________________
574 AliDetector *AliRun::GetDetector(const char *name) const
575 {
576   //
577   // Return pointer to detector from name
578   //
579   return dynamic_cast<AliDetector*>(fModules->FindObject(name));
580 }
581  
582 //_______________________________________________________________________
583 Int_t AliRun::GetModuleID(const char *name) const
584 {
585   //
586   // Return galice internal detector identifier from name
587   //
588   Int_t i=-1;
589   TObject *mod=fModules->FindObject(name);
590   if(mod) i=fModules->IndexOf(mod);
591   return i;
592 }
593  
594 //_______________________________________________________________________
595 Int_t AliRun::GetEvent(Int_t event)
596 {
597 //
598 // Reloads data containers in folders # event
599 // Set branch addresses
600 //
601   if (fRunLoader == 0x0)
602    {
603      Error("GetEvent","RunLoader is not set. Can not load data.");
604      return -1;
605    }
606 /*****************************************/ 
607 /****   P R E    R E L O A D I N G    ****/
608 /*****************************************/ 
609 // Reset existing structures
610   ResetHits();
611   ResetTrackReferences();
612   ResetDigits();
613   ResetSDigits();
614
615 /*****************************************/ 
616 /****       R  E  L  O  A  D          ****/
617 /*****************************************/
618
619   fRunLoader->GetEvent(event);
620
621 /*****************************************/ 
622 /****  P O S T    R E L O A D I N G   ****/
623 /*****************************************/ 
624
625   // Set Trees branch addresses
626   TIter next(fModules);
627   AliModule *detector;
628   while((detector = dynamic_cast<AliModule*>(next()))) 
629    {
630      detector->SetTreeAddress();
631    }
632  
633   return fRunLoader->GetHeader()->GetNtrack();
634 }
635
636 //_______________________________________________________________________
637 TGeometry *AliRun::GetGeometry()
638 {
639   //
640   // Import Alice geometry from current file
641   // Return pointer to geometry object
642   //
643   if (!fGeometry) fGeometry = dynamic_cast<TGeometry*>(gDirectory->Get("AliceGeom"));
644   //
645   // Unlink and relink nodes in detectors
646   // This is bad and there must be a better way...
647   //
648   
649   TIter next(fModules);
650   AliModule *detector;
651   while((detector = dynamic_cast<AliModule*>(next()))) {
652     TList *dnodes=detector->Nodes();
653     Int_t j;
654     TNode *node, *node1;
655     for ( j=0; j<dnodes->GetSize(); j++) {
656       node = dynamic_cast<TNode*>(dnodes->At(j));
657       node1 = fGeometry->GetNode(node->GetName());
658       dnodes->Remove(node);
659       dnodes->AddAt(node1,j);
660     }
661   }
662   return fGeometry;
663 }
664
665 //_______________________________________________________________________
666 Int_t AliRun::GetPrimary(Int_t track) const
667 {
668   //
669   // return number of primary that has generated track
670   //
671     return fRunLoader->Stack()->GetPrimary(track);
672 }
673  
674 //_______________________________________________________________________
675 void AliRun::MediaTable()
676 {
677   //
678   // Built media table to get from the media number to
679   // the detector id
680   //
681
682   Int_t kz, nz, idt, lz, i, k, ind;
683   //  Int_t ibeg;
684   TObjArray &dets = *gAlice->Detectors();
685   AliModule *det;
686   //
687   // For all detectors
688   for (kz=0;kz<fNdets;kz++) {
689     // If detector is defined
690     if((det=dynamic_cast<AliModule*>(dets[kz]))) {
691         TArrayI &idtmed = *(det->GetIdtmed()); 
692         for(nz=0;nz<100;nz++) {
693         // Find max and min material number
694         if((idt=idtmed[nz])) {
695           det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
696           det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
697         }
698       }
699       if(det->LoMedium() > det->HiMedium()) {
700         det->LoMedium() = 0;
701         det->HiMedium() = 0;
702       } else {
703         if(det->HiMedium() > fImedia->GetSize()) {
704           Error("MediaTable","Increase fImedia from %d to %d",
705                 fImedia->GetSize(),det->HiMedium());
706           return;
707         }
708         // Tag all materials in rage as belonging to detector kz
709         for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
710           (*fImedia)[lz]=kz;
711         }
712       }
713     }
714   }
715   //
716   // Print summary table
717   printf(" Traking media ranges:\n");
718   for(i=0;i<(fNdets-1)/6+1;i++) {
719     for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
720       ind=i*6+k;
721       det=dynamic_cast<AliModule*>(dets[ind]);
722       if(det)
723         printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
724                det->HiMedium());
725       else
726         printf(" %6s: %3d -> %3d;","NULL",0,0);
727     }
728     printf("\n");
729   }
730 }
731
732 //_______________________________________________________________________
733 void AliRun::SetGenerator(AliGenerator *generator)
734 {
735   //
736   // Load the event generator
737   //
738   if(!fGenerator) fGenerator = generator;
739 }
740
741 //_______________________________________________________________________
742 void AliRun::ResetGenerator(AliGenerator *generator)
743 {
744   //
745   // Load the event generator
746   //
747   if(fGenerator)
748     if(generator)
749       Warning("ResetGenerator","Replacing generator %s with %s\n",
750               fGenerator->GetName(),generator->GetName());
751     else
752       Warning("ResetGenerator","Replacing generator %s with NULL\n",
753               fGenerator->GetName());
754   fGenerator = generator;
755 }
756
757 //_______________________________________________________________________
758 void AliRun::SetTransPar(const char *filename)
759 {
760   //
761   // Sets the file name for transport parameters
762   //
763   fTransParName = filename;
764 }
765
766 //_______________________________________________________________________
767 void AliRun::SetBaseFile(const char *filename)
768 {
769   fBaseFileName = filename;
770 }
771
772 //_______________________________________________________________________
773 void AliRun::ReadTransPar()
774 {
775   //
776   // Read filename to set the transport parameters
777   //
778
779
780   const Int_t kncuts=10;
781   const Int_t knflags=11;
782   const Int_t knpars=kncuts+knflags;
783   const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
784                                "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
785                                "BREM","COMP","DCAY","DRAY","HADR","LOSS",
786                                "MULS","PAIR","PHOT","RAYL"};
787   char line[256];
788   char detName[7];
789   char* filtmp;
790   Float_t cut[kncuts];
791   Int_t flag[knflags];
792   Int_t i, itmed, iret, ktmed, kz;
793   FILE *lun;
794   //
795   // See whether the file is there
796   filtmp=gSystem->ExpandPathName(fTransParName.Data());
797   lun=fopen(filtmp,"r");
798   delete [] filtmp;
799   if(!lun) {
800     Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
801     return;
802   }
803   //
804   if(fDebug) {
805     printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
806     printf(" *%59s\n","*");
807     printf(" *       Please check carefully what you are doing!%10s\n","*");
808     printf(" *%59s\n","*");
809   }
810   //
811   while(1) {
812     // Initialise cuts and flags
813     for(i=0;i<kncuts;i++) cut[i]=-99;
814     for(i=0;i<knflags;i++) flag[i]=-99;
815     itmed=0;
816     for(i=0;i<256;i++) line[i]='\0';
817     // Read up to the end of line excluded
818     iret=fscanf(lun,"%[^\n]",line);
819     if(iret<0) {
820       //End of file
821       fclose(lun);
822       if(fDebug){
823         printf(" *%59s\n","*");
824         printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
825       }
826       return;
827     }
828     // Read the end of line
829     fscanf(lun,"%*c");
830     if(!iret) continue;
831     if(line[0]=='*') continue;
832     // Read the numbers
833     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",
834                 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
835                 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
836                 &flag[8],&flag[9],&flag[10]);
837     if(!iret) continue;
838     if(iret<0) {
839       //reading error
840       Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
841       continue;
842     }
843     // Check that the module exist
844     AliModule *mod = GetModule(detName);
845     if(mod) {
846       // Get the array of media numbers
847       TArrayI &idtmed = *mod->GetIdtmed();
848       // Check that the tracking medium code is valid
849       if(0<=itmed && itmed < 100) {
850         ktmed=idtmed[itmed];
851         if(!ktmed) {
852           Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
853           continue;
854         }
855         // Set energy thresholds
856         for(kz=0;kz<kncuts;kz++) {
857           if(cut[kz]>=0) {
858             if(fDebug) printf(" *  %-6s set to %10.3E for tracking medium code %4d for %s\n",
859                    kpars[kz],cut[kz],itmed,mod->GetName());
860             gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
861           }
862         }
863         // Set transport mechanisms
864         for(kz=0;kz<knflags;kz++) {
865           if(flag[kz]>=0) {
866             if(fDebug) printf(" *  %-6s set to %10d for tracking medium code %4d for %s\n",
867                    kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
868             gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
869           }
870         }
871       } else {
872         Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
873         continue;
874       }
875     } else {
876       if(fDebug) printf("%s::ReadTransParModule: %s not present\n",ClassName(),detName);
877       continue;
878     }
879   }
880 }
881 //_____________________________________________________________________________
882
883 void AliRun::BeginEvent()
884 {
885   //
886   // Clean-up previous event
887   // Energy scores
888   if (GetDebug()) 
889    {
890      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
891      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
892      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
893      Info("BeginEvent","          BEGINNING EVENT               ");
894      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
895      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
896      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
897    }
898     
899   /*******************************/    
900   /*   Clean after eventual      */
901   /*   previous event            */
902   /*******************************/    
903
904   
905   //Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors),
906   fRunLoader->SetEventNumber(++fEventNrInRun);// sets new files, cleans the previous event stuff, if necessary, etc.,  
907   if (GetDebug()) Info("BeginEvent","EventNr is %d",fEventNrInRun);
908      
909   fEventEnergy.Reset();  
910     // Clean detector information
911   
912   if (fRunLoader->Stack())
913     fRunLoader->Stack()->Reset();//clean stack -> tree is unloaded
914   else
915     fRunLoader->MakeStack();//or make a new one
916     
917   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTree(K)");
918   fRunLoader->MakeTree("K");
919   if (GetDebug()) Info("BeginEvent","  gMC->SetStack(fRunLoader->Stack())");
920   gMC->SetStack(fRunLoader->Stack());//Was in InitMC - but was moved here 
921                                      //because we don't have guarantee that 
922                                      //stack pointer is not going to change from event to event
923                          //since it bellobgs to header and is obtained via RunLoader
924   //
925   //  Reset all Detectors & kinematics & make/reset trees
926   //
927     
928   fRunLoader->GetHeader()->Reset(fRun,fEvent,fEventNrInRun);
929 //  fRunLoader->WriteKinematics("OVERWRITE");  is there any reason to rewrite here since MakeTree does so
930
931   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTrackRefsContainer()");
932   fRunLoader->MakeTrackRefsContainer();//for insurance
933
934   if (GetDebug()) Info("BeginEvent","  ResetHits()");
935   ResetHits();
936   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTree(H)");
937   fRunLoader->MakeTree("H");
938
939   //
940   if(fLego) 
941    {
942     fLego->BeginEvent();
943     return;
944    }
945
946   //create new branches and SetAdresses
947   TIter next(fModules);
948   AliModule *detector;
949   while((detector = (AliModule*)next()))
950    {
951     if (GetDebug()) Info("BeginEvent","  %s->MakeBranch(H)",detector->GetName());
952     detector->MakeBranch("H"); 
953     if (GetDebug()) Info("BeginEvent","  %s->MakeBranchTR()",detector->GetName());
954     detector->MakeBranchTR();
955     if (GetDebug()) Info("BeginEvent","  %s->SetTreeAddress()",detector->GetName());
956     detector->SetTreeAddress();
957    }
958 }
959
960 //_______________________________________________________________________
961 TParticle* AliRun::Particle(Int_t i) const
962 {
963   if (fRunLoader)
964    if (fRunLoader->Stack())
965     return fRunLoader->Stack()->Particle(i);
966   return 0x0;   
967 }
968
969 //_______________________________________________________________________
970 void AliRun::ResetDigits()
971 {
972   //
973   //  Reset all Detectors digits
974   //
975   TIter next(fModules);
976   AliModule *detector;
977   while((detector = dynamic_cast<AliModule*>(next()))) {
978      detector->ResetDigits();
979   }
980 }
981
982 //_______________________________________________________________________
983 void AliRun::ResetSDigits()
984 {
985   //
986   //  Reset all Detectors digits
987   //
988   TIter next(fModules);
989   AliModule *detector;
990   while((detector = dynamic_cast<AliModule*>(next()))) {
991      detector->ResetSDigits();
992   }
993 }
994
995 //_______________________________________________________________________
996 void AliRun::ResetHits()
997 {
998   //
999   //  Reset all Detectors hits
1000   //
1001   TIter next(fModules);
1002   AliModule *detector;
1003   while((detector = dynamic_cast<AliModule*>(next()))) {
1004      detector->ResetHits();
1005   }
1006 }
1007 //_______________________________________________________________________
1008
1009 void AliRun::ResetTrackReferences()
1010 {
1011   //
1012   //  Reset all Detectors hits
1013   //
1014   TIter next(fModules);
1015   AliModule *detector;
1016   while((detector = dynamic_cast<AliModule*>(next()))) {
1017      detector->ResetTrackReferences();
1018   }
1019 }
1020 //_______________________________________________________________________
1021
1022 void AliRun::ResetPoints()
1023 {
1024   //
1025   // Reset all Detectors points
1026   //
1027   TIter next(fModules);
1028   AliModule *detector;
1029   while((detector = dynamic_cast<AliModule*>(next()))) {
1030      detector->ResetPoints();
1031   }
1032 }
1033 //_______________________________________________________________________
1034
1035 void AliRun::InitMC(const char *setup)
1036 {
1037   //
1038   // Initialize the Alice setup
1039   //
1040   Announce();
1041
1042   if(fInitDone) {
1043     Warning("Init","Cannot initialise AliRun twice!\n");
1044     return;
1045   }
1046     
1047   gROOT->LoadMacro(setup);
1048   gInterpreter->ProcessLine(fConfigFunction.Data());
1049
1050   // Register MC in configuration 
1051   AliConfig::Instance()->Add(gMC);
1052
1053   InitLoaders();
1054
1055   fRunLoader->MakeTree("E");
1056   fRunLoader->LoadKinematics("RECREATE");
1057   fRunLoader->LoadTrackRefs("RECREATE");
1058   fRunLoader->LoadHits("all","RECREATE");
1059   
1060   
1061   fRunLoader->CdGAFile();
1062
1063   gMC->DefineParticles();  //Create standard MC particles
1064   AliPDG::AddParticlesToPdgDataBase();  
1065
1066   TObject *objfirst, *objlast;
1067
1068   fNdets = fModules->GetLast()+1;
1069
1070   //
1071   //=================Create Materials and geometry
1072   gMC->Init();
1073
1074   // Added also after in case of interactive initialisation of modules
1075   fNdets = fModules->GetLast()+1;
1076
1077   TIter next(fModules);
1078   AliModule *detector;
1079   while((detector = dynamic_cast<AliModule*>(next()))) 
1080    {
1081      objlast = gDirectory->GetList()->Last();
1082       
1083      // Add Detector histograms in Detector list of histograms
1084      if (objlast) objfirst = gDirectory->GetList()->After(objlast);
1085      else         objfirst = gDirectory->GetList()->First();
1086      while (objfirst) 
1087       {
1088         detector->Histograms()->Add(objfirst);
1089         objfirst = gDirectory->GetList()->After(objfirst);
1090       }
1091    }
1092    ReadTransPar(); //Read the cuts for all materials
1093    
1094    MediaTable(); //Build the special IMEDIA table
1095    
1096    //Initialise geometry deposition table
1097    fEventEnergy.Set(gMC->NofVolumes()+1);
1098    fSummEnergy.Set(gMC->NofVolumes()+1);
1099    fSum2Energy.Set(gMC->NofVolumes()+1);
1100    
1101    //Compute cross-sections
1102    gMC->BuildPhysics();
1103    
1104    //Write Geometry object to current file.
1105    fRunLoader->WriteGeometry();
1106    
1107    fInitDone = kTRUE;
1108
1109    fMCQA = new AliMCQA(fNdets);
1110
1111    //
1112    // Save stuff at the beginning of the file to avoid file corruption
1113    Write();
1114    fEventNrInRun = -1; //important - we start Begin event from increasing current number in run
1115 }
1116
1117 //_______________________________________________________________________
1118
1119 void AliRun::RunMC(Int_t nevent, const char *setup)
1120 {
1121   //
1122   // Main function to be called to process a galice run
1123   // example
1124   //    Root > gAlice.Run(); 
1125   // a positive number of events will cause the finish routine
1126   // to be called
1127   //
1128   fEventsPerRun = nevent;
1129   // check if initialisation has been done
1130   if (!fInitDone) InitMC(setup);
1131   
1132   // Create the Root Tree with one branch per detector
1133   //Hits moved to begin event -> now we are crating separate tree for each event
1134
1135   gMC->ProcessRun(nevent);
1136
1137   // End of this run, close files
1138   if(nevent>0) FinishRun();
1139 }
1140
1141 //_______________________________________________________________________
1142 void AliRun::RunReco(const char *selected, Int_t first, Int_t last)
1143 {
1144   //
1145   // Main function to be called to reconstruct Alice event
1146   // 
1147    Int_t nev = fRunLoader->GetNumberOfEvents();
1148    if (GetDebug()) Info("RunReco","Found %d events",nev);
1149    Int_t nFirst = first;
1150    Int_t nLast  = (last < 0)? nev : last;
1151    
1152    for (Int_t nevent = nFirst; nevent <= nLast; nevent++) {
1153      if (GetDebug()) Info("RunReco","Processing event %d",nevent);
1154      GetEvent(nevent);
1155      Digits2Reco(selected);
1156    }
1157 }
1158
1159 //_______________________________________________________________________
1160
1161 void AliRun::Hits2Digits(const char *selected)
1162 {
1163
1164    // Convert Hits to sumable digits
1165    // 
1166    for (Int_t nevent=0; nevent<gAlice->TreeE()->GetEntries(); nevent++) {
1167      GetEvent(nevent);
1168      Hits2SDigits(selected);
1169      SDigits2Digits(selected);
1170    }  
1171 }
1172
1173
1174 //_______________________________________________________________________
1175
1176 void AliRun::Tree2Tree(Option_t *option, const char *selected)
1177 {
1178   //
1179   // Function to transform the content of
1180   //  
1181   // - TreeH to TreeS (option "S")
1182   // - TreeS to TreeD (option "D")
1183   // - TreeD to TreeR (option "R")
1184   // 
1185   // If multiple options are specified ("SDR"), transformation will be done in sequence for
1186   // selected detector and for all detectors if none is selected (detector string 
1187   // can contain blank separated list of detector names). 
1188
1189
1190    const char *oS = strstr(option,"S");
1191    const char *oD = strstr(option,"D");
1192    const char *oR = strstr(option,"R");
1193    
1194    TObjArray *detectors = Detectors();
1195
1196    TIter next(detectors);
1197
1198    AliDetector *detector = 0;
1199
1200    while((detector = dynamic_cast<AliDetector*>(next()))) {
1201      if (selected) 
1202        if (strcmp(detector->GetName(),selected)) continue;
1203      if (detector->IsActive())
1204       { 
1205        
1206        AliLoader* loader = detector->GetLoader();
1207        if (loader == 0x0) continue;
1208        
1209        if (oS) 
1210         {
1211           if (GetDebug()) Info("Tree2Tree","Processing Hits2SDigits for %s ...",detector->GetName());
1212           loader->LoadHits("read");
1213           if (loader->TreeS() == 0x0) loader->MakeTree("S");
1214           detector->MakeBranch(option);
1215           detector->SetTreeAddress();
1216           detector->Hits2SDigits();
1217           loader->UnloadHits();
1218           loader->UnloadSDigits();
1219         }  
1220        if (oD) 
1221         {
1222           if (GetDebug()) Info("Tree2Tree","Processing SDigits2Digits for %s ...",detector->GetName());
1223           loader->LoadSDigits("read");
1224           if (loader->TreeD() == 0x0) loader->MakeTree("D");
1225           detector->MakeBranch(option);
1226           detector->SetTreeAddress();
1227           detector->SDigits2Digits();
1228           loader->UnloadSDigits();
1229           loader->UnloadDigits();
1230         } 
1231        if (oR) 
1232         {
1233           if (GetDebug()) Info("Tree2Tree","Processing Digits2Reco for %s ...",detector->GetName());
1234           loader->LoadDigits("read");
1235           if (loader->TreeR() == 0x0) loader->MakeTree("R");
1236           detector->MakeBranch(option);
1237           detector->SetTreeAddress();
1238           detector->Digits2Reco(); 
1239           loader->UnloadDigits();
1240           loader->UnloadRecPoints();
1241
1242         }
1243      }   
1244    }
1245 }
1246
1247 //_______________________________________________________________________
1248 void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min,
1249                      Float_t c1max,Int_t nc2,Float_t c2min,Float_t c2max,
1250                      Float_t rmin,Float_t rmax,Float_t zmax, AliLegoGenerator* gener)
1251 {
1252   //
1253   // Generates lego plots of:
1254   //    - radiation length map phi vs theta
1255   //    - radiation length map phi vs eta
1256   //    - interaction length map
1257   //    - g/cm2 length map
1258   //
1259   //  ntheta    bins in theta, eta
1260   //  themin    minimum angle in theta (degrees)
1261   //  themax    maximum angle in theta (degrees)
1262   //  nphi      bins in phi
1263   //  phimin    minimum angle in phi (degrees)
1264   //  phimax    maximum angle in phi (degrees)
1265   //  rmin      minimum radius
1266   //  rmax      maximum radius
1267   //  
1268   //
1269   //  The number of events generated = ntheta*nphi
1270   //  run input parameters in macro setup (default="Config.C")
1271   //
1272   //  Use macro "lego.C" to visualize the 3 lego plots in spherical coordinates
1273   //Begin_Html
1274   /*
1275     <img src="picts/AliRunLego1.gif">
1276   */
1277   //End_Html
1278   //Begin_Html
1279   /*
1280     <img src="picts/AliRunLego2.gif">
1281   */
1282   //End_Html
1283   //Begin_Html
1284   /*
1285     <img src="picts/AliRunLego3.gif">
1286   */
1287   //End_Html
1288   //
1289
1290   // check if initialisation has been done
1291   if (!fInitDone) InitMC(setup);
1292   //Save current generator
1293   AliGenerator *gen=Generator();
1294
1295   // Set new generator
1296   if (!gener) gener  = new AliLegoGenerator();
1297   ResetGenerator(gener);
1298   //
1299   // Configure Generator
1300   gener->SetRadiusRange(rmin, rmax);
1301   gener->SetZMax(zmax);
1302   gener->SetCoor1Range(nc1, c1min, c1max);
1303   gener->SetCoor2Range(nc2, c2min, c2max);
1304   
1305   
1306   //Create Lego object  
1307   fLego = new AliLego("lego",gener);
1308
1309   //Prepare MC for Lego Run
1310   gMC->InitLego();
1311   
1312   //Run Lego Object
1313
1314   //gMC->ProcessRun(nc1*nc2+1);
1315   gMC->ProcessRun(nc1*nc2);
1316   
1317   // Create only the Root event Tree
1318   fRunLoader->MakeTree("E");
1319   
1320   // End of this run, close files
1321   FinishRun();
1322   // Restore current generator
1323   ResetGenerator(gen);
1324   // Delete Lego Object
1325   delete fLego; fLego=0;
1326 }
1327
1328 //_______________________________________________________________________
1329 void AliRun::SetConfigFunction(const char * config) 
1330 {
1331   //
1332   // Set the signature of the function contained in Config.C to configure
1333   // the run
1334   //
1335   fConfigFunction=config;
1336 }
1337
1338 //_______________________________________________________________________
1339 void AliRun::SetCurrentTrack(Int_t track)
1340
1341   //
1342   // Set current track number
1343   //
1344     fRunLoader->Stack()->SetCurrentTrack(track); 
1345 }
1346  
1347 //_______________________________________________________________________
1348 void AliRun::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
1349                       Float_t *vpos, Float_t *polar, Float_t tof,
1350                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
1351
1352 // Delegate to stack
1353 //
1354     fRunLoader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
1355                      mech, ntr, weight, is);
1356 }
1357
1358 //_______________________________________________________________________
1359 void AliRun::PushTrack(Int_t done, Int_t parent, Int_t pdg,
1360                       Double_t px, Double_t py, Double_t pz, Double_t e,
1361                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
1362                       Double_t polx, Double_t poly, Double_t polz,
1363                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
1364
1365   // Delegate to stack
1366   //
1367   fRunLoader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
1368                                 polx, poly, polz, mech, ntr, weight, is);
1369 }
1370
1371 //_______________________________________________________________________
1372 void AliRun::SetHighWaterMark(const Int_t nt)
1373 {
1374     //
1375     // Set high water mark for last track in event
1376     fRunLoader->Stack()->SetHighWaterMark(nt);
1377 }
1378
1379 //_______________________________________________________________________
1380 void AliRun::KeepTrack(const Int_t track)
1381
1382   //
1383   // Delegate to stack
1384   //
1385     fRunLoader->Stack()->KeepTrack(track);
1386 }
1387  
1388 // 
1389 // MC Application
1390 // 
1391
1392 //_______________________________________________________________________
1393 void  AliRun::ConstructGeometry() 
1394 {
1395   //
1396   // Create modules, materials, geometry
1397   //
1398
1399     TStopwatch stw;
1400     TIter next(fModules);
1401     AliModule *detector;
1402     if (GetDebug()) Info("ConstructGeometry","Geometry creation:");
1403     while((detector = dynamic_cast<AliModule*>(next()))) {
1404       stw.Start();
1405       // Initialise detector materials and geometry
1406       detector->CreateMaterials();
1407       detector->CreateGeometry();
1408       printf("%10s R:%.2fs C:%.2fs\n",
1409              detector->GetName(),stw.RealTime(),stw.CpuTime());
1410     }
1411 }
1412
1413 //_______________________________________________________________________
1414 void  AliRun::InitGeometry()
1415
1416   //
1417   // Initialize detectors and display geometry
1418   //
1419
1420    printf("Initialisation:\n");
1421     TStopwatch stw;
1422     TIter next(fModules);
1423     AliModule *detector;
1424     while((detector = dynamic_cast<AliModule*>(next()))) {
1425       stw.Start();
1426       // Initialise detector and display geometry
1427       detector->Init();
1428       detector->BuildGeometry();
1429       printf("%10s R:%.2fs C:%.2fs\n",
1430              detector->GetName(),stw.RealTime(),stw.CpuTime());
1431     }
1432  
1433 }
1434 //_______________________________________________________________________
1435
1436 void  AliRun::GeneratePrimaries()
1437
1438   //
1439   // Generate primary particles and fill them in the stack.
1440   //
1441
1442   Generator()->Generate();
1443 }
1444 //_______________________________________________________________________
1445
1446 void AliRun::BeginPrimary()
1447 {
1448   //
1449   // Called  at the beginning of each primary track
1450   //
1451   
1452   // Reset Hits info
1453   gAlice->ResetHits();
1454   gAlice->ResetTrackReferences();
1455
1456 }
1457
1458 //_______________________________________________________________________
1459 void AliRun::PreTrack()
1460 {
1461      TObjArray &dets = *fModules;
1462      AliModule *module;
1463
1464      for(Int_t i=0; i<=fNdets; i++)
1465        if((module = dynamic_cast<AliModule*>(dets[i])))
1466          module->PreTrack();
1467
1468      fMCQA->PreTrack();
1469 }
1470
1471 //_______________________________________________________________________
1472 void AliRun::Stepping() 
1473 {
1474   //
1475   // Called at every step during transport
1476   //
1477
1478   Int_t id = DetFromMate(gMC->GetMedium());
1479   if (id < 0) return;
1480
1481   //
1482   // --- If lego option, do it and leave 
1483   if (fLego)
1484     fLego->StepManager();
1485   else {
1486     Int_t copy;
1487     //Update energy deposition tables
1488     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
1489   
1490     //Call the appropriate stepping routine;
1491     AliModule *det = dynamic_cast<AliModule*>(fModules->At(id));
1492     if(det && det->StepManagerIsEnabled()) {
1493       fMCQA->StepManager(id);
1494       det->StepManager();
1495     }
1496   }
1497 }
1498
1499 //_______________________________________________________________________
1500 void AliRun::PostTrack()
1501 {
1502      TObjArray &dets = *fModules;
1503      AliModule *module;
1504
1505      for(Int_t i=0; i<=fNdets; i++)
1506        if((module = dynamic_cast<AliModule*>(dets[i])))
1507          module->PostTrack();
1508 }
1509
1510 //_______________________________________________________________________
1511 void AliRun::FinishPrimary()
1512 {
1513   //
1514   // Called  at the end of each primary track
1515   //
1516   
1517   //  static Int_t count=0;
1518   //  const Int_t times=10;
1519   // This primary is finished, purify stack
1520   fRunLoader->Stack()->PurifyKine();
1521
1522   TIter next(fModules);
1523   AliModule *detector;
1524   while((detector = dynamic_cast<AliModule*>(next()))) {
1525     detector->FinishPrimary();
1526     if(detector->GetLoader())
1527      {
1528        detector->GetLoader()->TreeH()->Fill();
1529      }
1530   }
1531
1532   // Write out track references if any
1533   if (fRunLoader->TreeTR()) 
1534    {
1535     fRunLoader->TreeTR()->Fill();
1536    }
1537 }
1538
1539 //_______________________________________________________________________
1540 void AliRun::FinishEvent()
1541 {
1542   //
1543   // Called at the end of the event.
1544   //
1545   
1546   //
1547   if(fLego) fLego->FinishEvent();
1548
1549   TIter next(fModules);
1550   AliModule *detector;
1551   while((detector = dynamic_cast<AliModule*>(next()))) {
1552     detector->FinishEvent();
1553   }
1554
1555   //Update the energy deposit tables
1556   Int_t i;
1557   for(i=0;i<fEventEnergy.GetSize();i++) 
1558    {
1559     fSummEnergy[i]+=fEventEnergy[i];
1560     fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
1561    }
1562
1563   AliHeader* header = fRunLoader->GetHeader();
1564   AliStack* stack = fRunLoader->Stack();
1565   if ( (header == 0x0) || (stack == 0x0) )
1566    {//check if we got header and stack. If not cry and exit aliroot
1567     Fatal("AliRun","Can not get the stack or header from LOADER");
1568     return;//never reached
1569    }  
1570   // Update Header information 
1571   header->SetNprimary(stack->GetNprimary());
1572   header->SetNtrack(stack->GetNtrack());  
1573
1574   
1575   // Write out the kinematics
1576   stack->FinishEvent();
1577    
1578   // Write out the event Header information
1579   TTree* treeE = fRunLoader->TreeE();
1580   if (treeE) 
1581    {
1582       header->SetStack(stack);
1583       treeE->Fill();
1584    }
1585   else
1586    {
1587     Error("FinishEvent","Can not get TreeE from RL");
1588    }
1589   
1590   fRunLoader->WriteKinematics("OVERWRITE");
1591   fRunLoader->WriteTrackRefs("OVERWRITE");
1592   fRunLoader->WriteHits("OVERWRITE");
1593
1594   if (GetDebug()) 
1595    { 
1596      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1597      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1598      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1599      Info("FinishEvent","          FINISHING EVENT               ");
1600      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1601      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1602      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
1603    }
1604 }
1605
1606 //_______________________________________________________________________
1607 void AliRun::Field(const Double_t* x, Double_t *b) const
1608 {
1609    Float_t xfloat[3];
1610    for (Int_t i=0; i<3; i++) xfloat[i] = x[i]; 
1611
1612    if (Field()) {
1613          Float_t bfloat[3];
1614          Field()->Field(xfloat,bfloat);
1615          for (Int_t j=0; j<3; j++) b[j] = bfloat[j]; 
1616    } 
1617    else {
1618          printf("No mag field defined!\n");
1619          b[0]=b[1]=b[2]=0.;
1620    }
1621
1622 }      
1623
1624 // 
1625 // End of MC Application
1626 // 
1627
1628 //_______________________________________________________________________
1629 void AliRun::Streamer(TBuffer &R__b)
1630 {
1631   // Stream an object of class AliRun.
1632
1633   if (R__b.IsReading()) {
1634     if (!gAlice) gAlice = this;
1635     AliRun::Class()->ReadBuffer(R__b, this);
1636     gROOT->GetListOfBrowsables()->Add(this,"Run");
1637
1638     gRandom = fRandom;
1639   } else {
1640     AliRun::Class()->WriteBuffer(R__b, this);
1641   }
1642 }
1643
1644
1645 //_______________________________________________________________________
1646 Int_t AliRun::GetCurrentTrackNumber() const {
1647   //
1648   // Returns current track
1649   //
1650   return fRunLoader->Stack()->GetCurrentTrackNumber();
1651 }
1652
1653 //_______________________________________________________________________
1654 Int_t AliRun::GetNtrack() const {
1655   //
1656   // Returns number of tracks in stack
1657   //
1658   return fRunLoader->Stack()->GetNtrack();
1659 }
1660 //_______________________________________________________________________
1661
1662 //_______________________________________________________________________
1663 TObjArray* AliRun::Particles() const {
1664   //
1665   // Returns pointer to Particles array
1666   //
1667   if (fRunLoader)
1668    if (fRunLoader->Stack())
1669     return fRunLoader->Stack()->Particles();
1670   return 0x0;
1671 }
1672
1673 //___________________________________________________________________________
1674
1675 //_______________________________________________________________________
1676 void AliRun::SetGenEventHeader(AliGenEventHeader* header)
1677 {
1678   fRunLoader->GetHeader()->SetGenEventHeader(header);
1679 }
1680
1681 //___________________________________________________________________________
1682
1683 Int_t AliRun::GetEvNumber() const
1684
1685 //Returns number of current event  
1686   if (fRunLoader == 0x0)
1687    {
1688      Error("GetEvent","RunLoader is not set. Can not load data.");
1689      return -1;
1690    }
1691
1692   return fRunLoader->GetEventNumber();
1693 }
1694
1695 void AliRun::SetRunLoader(AliRunLoader* rloader)
1696 {
1697   fRunLoader = rloader;
1698   if (fRunLoader == 0x0) return;
1699   
1700   TString evfoldname;
1701   TFolder* evfold = fRunLoader->GetEventFolder();
1702   if (evfold) evfoldname = evfold->GetName();
1703   else Warning("SetRunLoader","Did not get Event Folder from Run Loader");
1704   
1705   if ( fRunLoader->GetAliRun() )
1706    {//if alrun already exists in folder
1707     if (fRunLoader->GetAliRun() != this )
1708      {//and is different than this - crash
1709        Fatal("AliRun","AliRun is already in Folder and it is not this object");
1710        return;//pro forma
1711      }//else do nothing
1712    }
1713   else
1714    {
1715      evfold->Add(this);//Post this AliRun to Folder
1716    }
1717   
1718   TIter next(fModules);
1719   AliModule *module;
1720   while((module = (AliModule*)next())) 
1721    {
1722      if (evfold) AliConfig::Instance()->Add(module,evfoldname);
1723      AliDetector* detector = dynamic_cast<AliDetector*>(module);
1724      if (detector)
1725       {
1726         AliLoader* loader = fRunLoader->GetLoader(detector);
1727         if (loader == 0x0)
1728          {
1729            Error("SetRunLoader","Can not get loader for detector %s",detector->GetName());
1730          }
1731         else
1732          {
1733            if (GetDebug()) Info("SetRunLoader","Setting loader for detector %s",detector->GetName());
1734            detector->SetLoader(loader);
1735          }
1736       }
1737    }
1738 }
1739
1740 void AliRun::AddModule(AliModule* mod)
1741 {
1742   if (mod == 0x0) return;
1743   if (strlen(mod->GetName()) == 0) return;
1744   if (GetModuleID(mod->GetName()) >= 0) return;
1745   
1746   if (GetDebug()) Info("AddModule","%s",mod->GetName());
1747   if (fRunLoader == 0x0) AliConfig::Instance()->Add(mod);
1748   else AliConfig::Instance()->Add(mod,fRunLoader->GetEventFolder()->GetName());
1749
1750   Modules()->Add(mod);
1751 }