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