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