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