STEER/AliRun.cxx

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