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