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