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