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