STEER/AliRun.cxx

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