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