STEER/AliRun.cxx

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