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