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.31  2000/04/26 10:17:32  fca
  19 Changes in Lego for G4 compatibility
  20
  21 Revision 1.30  2000/04/18 19:11:40  fca
  22 Introduce variable Config.C function signature
  23
  24 Revision 1.29  2000/04/07 11:12:34  fca
  25 G4 compatibility changes
  26
  27 Revision 1.28  2000/04/05 06:51:06  fca
  28 Workaround for an HP compiler problem
  29
  30 Revision 1.27  2000/03/22 18:08:07  fca
  31 Rationalisation of the virtual MC interfaces
  32
  33 Revision 1.26  2000/03/22 13:42:26  fca
  34 SetGenerator does not replace an existing generator, ResetGenerator does
  35
  36 Revision 1.25  2000/02/23 16:25:22  fca
  37 AliVMC and AliGeant3 classes introduced
  38 ReadEuclid moved from AliRun to AliModule
  39
  40 Revision 1.24  2000/01/19 17:17:20  fca
  41 Introducing a list of lists of hits -- more hits allowed for detector now
  42
  43 Revision 1.23  1999/12/03 11:14:31  fca
  44 Fixing previous wrong checking
  45
  46 Revision 1.21  1999/11/25 10:40:08  fca
  47 Fixing daughters information also in primary tracks
  48
  49 Revision 1.20  1999/10/04 18:08:49  fca
  50 Adding protection against inconsistent Euclid files
  51
  52 Revision 1.19  1999/09/29 07:50:40  fca
  53 Introduction of the Copyright and cvs Log
  54
  55 */
  56
  57 ///////////////////////////////////////////////////////////////////////////////
  58 //                                                                           //
  59 //  Control class for Alice C++                                              //
  60 //  Only one single instance of this class exists.                           //
  61 //  The object is created in main program aliroot                            //
  62 //  and is pointed by the global gAlice.                                     //
  63 //                                                                           //
  64 //   -Supports the list of all Alice Detectors (fModules).                 //
  65 //   -Supports the list of particles (fParticles).                           //
  66 //   -Supports the Trees.                                                    //
  67 //   -Supports the geometry.                                                 //
  68 //   -Supports the event display.                                            //
  69 //Begin_Html
  70 /*
  71 <img src="picts/AliRunClass.gif">
  72 */
  73 //End_Html
  74 //Begin_Html
  75 /*
  76 <img src="picts/alirun.gif">
  77 */
  78 //End_Html
  79 //                                                                           //
  80 ///////////////////////////////////////////////////////////////////////////////
  81
  82 #include <TFile.h>
  83 #include <TRandom.h>
  84 #include <TBRIK.h>
  85 #include <TNode.h>
  86 #include <TCint.h>
  87 #include <TSystem.h>
  88 #include <TObjectTable.h>
  89
  90 #include "TParticle.h"
  91 #include "AliRun.h"
  92 #include "AliDisplay.h"
  93 #include "AliMC.h"
  94 #include "AliLego.h"
  95
  96 #include <stdlib.h>
  97 #include <stdio.h>
  98 #include <string.h>
  99
 100 AliRun *gAlice;
 101
 102 static AliHeader *header;
 103
 104 ClassImp(AliRun)
 105
 106 //_____________________________________________________________________________
 107 AliRun::AliRun()
 108 {
 109   //
 110   // Default constructor for AliRun
 111   //
 112   header=&fHeader;
 113   fRun       = 0;
 114   fEvent     = 0;
 115   fCurrent   = -1;
 116   fModules = 0;
 117   fGenerator = 0;
 118   fTreeD     = 0;
 119   fTreeK     = 0;
 120   fTreeH     = 0;
 121   fTreeE     = 0;
 122   fTreeR     = 0;
 123   fParticles = 0;
 124   fGeometry  = 0;
 125   fDisplay   = 0;
 126   fField     = 0;
 127   fMC       = 0;
 128   fNdets     = 0;
 129   fImedia    = 0;
 130   fTrRmax    = 1.e10;
 131   fTrZmax    = 1.e10;
 132   fInitDone  = kFALSE;
 133   fLego      = 0;
 134   fPDGDB     = 0;        //Particle factory object!
 135   fHitLists  = 0;
 136   fConfigFunction    = "\0";
 137 }
 138
 139 //_____________________________________________________________________________
 140 AliRun::AliRun(const char *name, const char *title)
 141   : TNamed(name,title)
 142 {
 143   //
 144   //  Constructor for the main processor.
 145   //  Creates the geometry
 146   //  Creates the list of Detectors.
 147   //  Creates the list of particles.
 148   //
 149   Int_t i;
 150
 151   gAlice     = this;
 152   fTreeD     = 0;
 153   fTreeK     = 0;
 154   fTreeH     = 0;
 155   fTreeE     = 0;
 156   fTreeR     = 0;
 157   fTrRmax    = 1.e10;
 158   fTrZmax    = 1.e10;
 159   fGenerator = 0;
 160   fInitDone  = kFALSE;
 161   fLego      = 0;
 162   fField     = 0;
 163   fConfigFunction    = "Config();";
 164
 165   gROOT->GetListOfBrowsables()->Add(this,name);
 166   //
 167   // create the support list for the various Detectors
 168   fModules = new TObjArray(77);
 169   //
 170   // Create the TNode geometry for the event display
 171
 172   BuildSimpleGeometry();
 173
 174
 175   fNtrack=0;
 176   fHgwmk=0;
 177   fCurrent=-1;
 178   header=&fHeader;
 179   fRun       = 0;
 180   fEvent     = 0;
 181   //
 182   // Create the particle stack
 183   fParticles = new TClonesArray("TParticle",100);
 184
 185   fDisplay = 0;
 186   //
 187   // Create default mag field
 188   SetField();
 189   //
 190   fMC      = gMC;
 191   //
 192   // Prepare the tracking medium lists
 193   fImedia = new TArrayI(1000);
 194   for(i=0;i<1000;i++) (*fImedia)[i]=-99;
 195   //
 196   // Make particles
 197   fPDGDB     = TDatabasePDG::Instance();        //Particle factory object!
 198   //
 199   // Create HitLists list
 200   fHitLists  = new TList();
 201 }
 202
 203 //_____________________________________________________________________________
 204 AliRun::~AliRun()
 205 {
 206   //
 207   // Defaullt AliRun destructor
 208   //
 209   delete fImedia;
 210   delete fField;
 211   delete fMC;
 212   delete fGeometry;
 213   delete fDisplay;
 214   delete fGenerator;
 215   delete fLego;
 216   delete fTreeD;
 217   delete fTreeK;
 218   delete fTreeH;
 219   delete fTreeE;
 220   delete fTreeR;
 221   if (fModules) {
 222     fModules->Delete();
 223     delete fModules;
 224   }
 225   if (fParticles) {
 226     fParticles->Delete();
 227     delete fParticles;
 228   }
 229   delete fHitLists;
 230 }
 231
 232 //_____________________________________________________________________________
 233 void AliRun::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
 234 {
 235   //
 236   //  Add a hit to detector id
 237   //
 238   TObjArray &dets = *fModules;
 239   if(dets[id]) ((AliModule*) dets[id])->AddHit(track,vol,hits);
 240 }
 241
 242 //_____________________________________________________________________________
 243 void AliRun::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
 244 {
 245   //
 246   // Add digit to detector id
 247   //
 248   TObjArray &dets = *fModules;
 249   if(dets[id]) ((AliModule*) dets[id])->AddDigit(tracks,digits);
 250 }
 251
 252 //_____________________________________________________________________________
 253 void AliRun::Browse(TBrowser *b)
 254 {
 255   //
 256   // Called when the item "Run" is clicked on the left pane
 257   // of the Root browser.
 258   // It displays the Root Trees and all detectors.
 259   //
 260   if (fTreeK) b->Add(fTreeK,fTreeK->GetName());
 261   if (fTreeH) b->Add(fTreeH,fTreeH->GetName());
 262   if (fTreeD) b->Add(fTreeD,fTreeD->GetName());
 263   if (fTreeE) b->Add(fTreeE,fTreeE->GetName());
 264   if (fTreeR) b->Add(fTreeR,fTreeR->GetName());
 265
 266   TIter next(fModules);
 267   AliModule *detector;
 268   while((detector = (AliModule*)next())) {
 269     b->Add(detector,detector->GetName());
 270   }
 271 }
 272
 273 //_____________________________________________________________________________
 274 void AliRun::Build()
 275 {
 276   //
 277   // Initialize Alice geometry
 278   // Dummy routine
 279   //
 280 }
 281
 282 //_____________________________________________________________________________
 283 void AliRun::BuildSimpleGeometry()
 284 {
 285   //
 286   // Create a simple TNode geometry used by Root display engine
 287   //
 288   // Initialise geometry
 289   //
 290   fGeometry = new TGeometry("AliceGeom","Galice Geometry for Hits");
 291   new TMaterial("void","Vacuum",0,0,0);  //Everything is void
 292   TBRIK *brik = new TBRIK("S_alice","alice volume","void",2000,2000,3000);
 293   brik->SetVisibility(0);
 294   new TNode("alice","alice","S_alice");
 295 }
 296
 297 //_____________________________________________________________________________
 298 void AliRun::CleanDetectors()
 299 {
 300   //
 301   // Clean Detectors at the end of event
 302   //
 303   TIter next(fModules);
 304   AliModule *detector;
 305   while((detector = (AliModule*)next())) {
 306     detector->FinishEvent();
 307   }
 308 }
 309
 310 //_____________________________________________________________________________
 311 void AliRun::CleanParents()
 312 {
 313   //
 314   // Clean Particles stack.
 315   // Set parent/daughter relations
 316   //
 317   TClonesArray &particles = *(gAlice->Particles());
 318   TParticle *part;
 319   int i;
 320   for(i=0; i<fNtrack; i++) {
 321     part = (TParticle *)particles.UncheckedAt(i);
 322     if(!part->TestBit(Daughters_Bit)) {
 323       part->SetFirstDaughter(-1);
 324       part->SetLastDaughter(-1);
 325     }
 326   }
 327 }
 328
 329 //_____________________________________________________________________________
 330 Int_t AliRun::DistancetoPrimitive(Int_t, Int_t)
 331 {
 332   //
 333   // Return the distance from the mouse to the AliRun object
 334   // Dummy routine
 335   //
 336   return 9999;
 337 }
 338
 339 //_____________________________________________________________________________
 340 void AliRun::DumpPart (Int_t i)
 341 {
 342   //
 343   // Dumps particle i in the stack
 344   //
 345   TClonesArray &particles = *fParticles;
 346   ((TParticle*) particles[i])->Print();
 347 }
 348
 349 //_____________________________________________________________________________
 350 void AliRun::DumpPStack ()
 351 {
 352   //
 353   // Dumps the particle stack
 354   //
 355   TClonesArray &particles = *fParticles;
 356   printf(
 357          "\n\n=======================================================================\n");
 358   for (Int_t i=0;i<fNtrack;i++)
 359     {
 360       printf("-> %d ",i); ((TParticle*) particles[i])->Print();
 361       printf("--------------------------------------------------------------\n");
 362     }
 363   printf(
 364          "\n=======================================================================\n\n");
 365 }
 366
 367 //_____________________________________________________________________________
 368 void AliRun::SetField(Int_t type, Int_t version, Float_t scale,
 369                       Float_t maxField, char* filename)
 370 {
 371   //
 372   //  Set magnetic field parameters
 373   //  type      Magnetic field transport flag 0=no field, 2=helix, 3=Runge Kutta
 374   //  version   Magnetic field map version (only 1 active now)
 375   //  scale     Scale factor for the magnetic field
 376   //  maxField  Maximum value for the magnetic field
 377
 378   //
 379   // --- Sanity check on mag field flags
 380   if(type<0 || type > 2) {
 381     Warning("SetField",
 382             "Invalid magnetic field flag: %5d; Helix tracking chosen instead\n"
 383            ,type);
 384     type=2;
 385   }
 386   if(fField) delete fField;
 387   if(version==1) {
 388     fField = new AliMagFC("Map1"," ",type,version,scale,maxField);
 389   } else if(version<=3) {
 390     fField = new AliMagFCM("Map2-3",filename,type,version,scale,maxField);
 391     fField->ReadField();
 392   } else {
 393     Warning("SetField","Invalid map %d\n",version);
 394   }
 395 }
 396
 397 //_____________________________________________________________________________
 398 void AliRun::FillTree()
 399 {
 400   //
 401   // Fills all AliRun TTrees
 402   //
 403   if (fTreeK) fTreeK->Fill();
 404   if (fTreeH) fTreeH->Fill();
 405   if (fTreeD) fTreeD->Fill();
 406   if (fTreeR) fTreeR->Fill();
 407 }
 408
 409 //_____________________________________________________________________________
 410 void AliRun::FinishPrimary()
 411 {
 412   //
 413   // Called  at the end of each primary track
 414   //
 415
 416   //  static Int_t count=0;
 417   //  const Int_t times=10;
 418   // This primary is finished, purify stack
 419   PurifyKine();
 420
 421   // Write out hits if any
 422   if (gAlice->TreeH()) {
 423     gAlice->TreeH()->Fill();
 424   }
 425
 426   // Reset Hits info
 427   gAlice->ResetHits();
 428
 429   //
 430   //  if(++count%times==1) gObjectTable->Print();
 431 }
 432
 433 //_____________________________________________________________________________
 434 void AliRun::FinishEvent()
 435 {
 436   //
 437   // Called at the end of the event.
 438   //
 439
 440   //
 441   if(fLego) fLego->FinishEvent();
 442
 443   //Update the energy deposit tables
 444   Int_t i;
 445   for(i=0;i<fEventEnergy.GetSize();i++) {
 446     fSummEnergy[i]+=fEventEnergy[i];
 447     fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
 448   }
 449   fEventEnergy.Reset();
 450
 451   // Clean detector information
 452   CleanDetectors();
 453
 454   // Write out the kinematics
 455   if (fTreeK) {
 456     CleanParents();
 457     fTreeK->Fill();
 458   }
 459
 460   // Write out the digits
 461   if (fTreeD) {
 462     fTreeD->Fill();
 463     ResetDigits();
 464   }
 465
 466   // Write out reconstructed clusters
 467   if (fTreeR) {
 468     fTreeR->Fill();
 469   }
 470
 471   // Write out the event Header information
 472   if (fTreeE) fTreeE->Fill();
 473
 474   // Reset stack info
 475   ResetStack();
 476
 477   // Write Tree headers
 478   //  Int_t ievent = fHeader.GetEvent();
 479   //  char hname[30];
 480   //  sprintf(hname,"TreeK%d",ievent);
 481   if (fTreeK) fTreeK->Write();
 482   //  sprintf(hname,"TreeH%d",ievent);
 483   if (fTreeH) fTreeH->Write();
 484   //  sprintf(hname,"TreeD%d",ievent);
 485   if (fTreeD) fTreeD->Write();
 486   //  sprintf(hname,"TreeR%d",ievent);
 487   if (fTreeR) fTreeR->Write();
 488
 489   ++fEvent;
 490 }
 491
 492 //_____________________________________________________________________________
 493 void AliRun::FinishRun()
 494 {
 495   //
 496   // Called at the end of the run.
 497   //
 498
 499   //
 500   if(fLego) fLego->FinishRun();
 501
 502   // Clean detector information
 503   TIter next(fModules);
 504   AliModule *detector;
 505   while((detector = (AliModule*)next())) {
 506     detector->FinishRun();
 507   }
 508
 509   //Output energy summary tables
 510   EnergySummary();
 511
 512   // file is retrieved from whatever tree
 513   TFile *File = 0;
 514   if (fTreeK) File = fTreeK->GetCurrentFile();
 515   if ((!File) && (fTreeH)) File = fTreeH->GetCurrentFile();
 516   if ((!File) && (fTreeD)) File = fTreeD->GetCurrentFile();
 517   if ((!File) && (fTreeE)) File = fTreeE->GetCurrentFile();
 518   if( NULL==File ) {
 519     Error("FinishRun","There isn't root file!");
 520     exit(1);
 521   }
 522   File->cd();
 523   fTreeE->Write();
 524
 525   // Clean tree information
 526   delete fTreeK; fTreeK = 0;
 527   delete fTreeH; fTreeH = 0;
 528   delete fTreeD; fTreeD = 0;
 529   delete fTreeR; fTreeR = 0;
 530   delete fTreeE; fTreeE = 0;
 531
 532   // Write AliRun info and all detectors parameters
 533   Write();
 534
 535   // Close output file
 536   File->Write();
 537 }
 538
 539 //_____________________________________________________________________________
 540 void AliRun::FlagTrack(Int_t track)
 541 {
 542   //
 543   // Flags a track and all its family tree to be kept
 544   //
 545   int curr;
 546   TParticle *particle;
 547
 548   curr=track;
 549   while(1) {
 550     particle=(TParticle*)fParticles->UncheckedAt(curr);
 551
 552     // If the particle is flagged the three from here upward is saved already
 553     if(particle->TestBit(Keep_Bit)) return;
 554
 555     // Save this particle
 556     particle->SetBit(Keep_Bit);
 557
 558     // Move to father if any
 559     if((curr=particle->GetFirstMother())==-1) return;
 560   }
 561 }
 562
 563 //_____________________________________________________________________________
 564 void AliRun::EnergySummary()
 565 {
 566   //
 567   // Print summary of deposited energy
 568   //
 569
 570   Int_t ndep=0;
 571   Float_t edtot=0;
 572   Float_t ed, ed2;
 573   Int_t kn, i, left, j, id;
 574   const Float_t zero=0;
 575   Int_t ievent=fHeader.GetEvent()+1;
 576   //
 577   // Energy loss information
 578   if(ievent) {
 579     printf("***************** Energy Loss Information per event (GEV) *****************\n");
 580     for(kn=1;kn<fEventEnergy.GetSize();kn++) {
 581       ed=fSummEnergy[kn];
 582       if(ed>0) {
 583         fEventEnergy[ndep]=kn;
 584         if(ievent>1) {
 585           ed=ed/ievent;
 586           ed2=fSum2Energy[kn];
 587           ed2=ed2/ievent;
 588           ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,zero))/ed;
 589         } else
 590           ed2=99;
 591         fSummEnergy[ndep]=ed;
 592         fSum2Energy[ndep]=TMath::Min((Float_t) 99.,TMath::Max(ed2,zero));
 593         edtot+=ed;
 594         ndep++;
 595       }
 596     }
 597     for(kn=0;kn<(ndep-1)/3+1;kn++) {
 598       left=ndep-kn*3;
 599       for(i=0;i<(3<left?3:left);i++) {
 600         j=kn*3+i;
 601         id=Int_t (fEventEnergy[j]+0.1);
 602         printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
 603       }
 604       printf("\n");
 605     }
 606     //
 607     // Relative energy loss in different detectors
 608     printf("******************** Relative Energy Loss per event ********************\n");
 609     printf("Total energy loss per event %10.3f GeV\n",edtot);
 610     for(kn=0;kn<(ndep-1)/5+1;kn++) {
 611       left=ndep-kn*5;
 612       for(i=0;i<(5<left?5:left);i++) {
 613         j=kn*5+i;
 614         id=Int_t (fEventEnergy[j]+0.1);
 615         printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
 616       }
 617       printf("\n");
 618     }
 619     for(kn=0;kn<75;kn++) printf("*");
 620     printf("\n");
 621   }
 622   //
 623   // Reset the TArray's
 624   //  fEventEnergy.Set(0);
 625   //  fSummEnergy.Set(0);
 626   //  fSum2Energy.Set(0);
 627 }
 628
 629 //_____________________________________________________________________________
 630 AliModule *AliRun::GetModule(const char *name)
 631 {
 632   //
 633   // Return pointer to detector from name
 634   //
 635   return (AliModule*)fModules->FindObject(name);
 636 }
 637
 638 //_____________________________________________________________________________
 639 AliDetector *AliRun::GetDetector(const char *name)
 640 {
 641   //
 642   // Return pointer to detector from name
 643   //
 644   return (AliDetector*)fModules->FindObject(name);
 645 }
 646
 647 //_____________________________________________________________________________
 648 Int_t AliRun::GetModuleID(const char *name)
 649 {
 650   //
 651   // Return galice internal detector identifier from name
 652   //
 653   Int_t i=-1;
 654   TObject *mod=fModules->FindObject(name);
 655   if(mod) i=fModules->IndexOf(mod);
 656   return i;
 657 }
 658
 659 //_____________________________________________________________________________
 660 Int_t AliRun::GetEvent(Int_t event)
 661 {
 662   //
 663   // Connect the Trees Kinematics and Hits for event # event
 664   // Set branch addresses
 665   //
 666
 667   // Reset existing structures
 668   ResetStack();
 669   ResetHits();
 670   ResetDigits();
 671
 672   // Delete Trees already connected
 673   if (fTreeK) delete fTreeK;
 674   if (fTreeH) delete fTreeH;
 675   if (fTreeD) delete fTreeD;
 676   if (fTreeR) delete fTreeR;
 677
 678   // Get header from file
 679   if(fTreeE) fTreeE->GetEntry(event);
 680   else Error("GetEvent","Cannot file Header Tree\n");
 681
 682   // Get Kine Tree from file
 683   char treeName[20];
 684   sprintf(treeName,"TreeK%d",event);
 685   fTreeK = (TTree*)gDirectory->Get(treeName);
 686   if (fTreeK) fTreeK->SetBranchAddress("Particles", &fParticles);
 687   else    Error("GetEvent","cannot find Kine Tree for event:%d\n",event);
 688
 689   // Get Hits Tree header from file
 690   sprintf(treeName,"TreeH%d",event);
 691   fTreeH = (TTree*)gDirectory->Get(treeName);
 692   if (!fTreeH) {
 693     Error("GetEvent","cannot find Hits Tree for event:%d\n",event);
 694   }
 695
 696   // Get Digits Tree header from file
 697   sprintf(treeName,"TreeD%d",event);
 698   fTreeD = (TTree*)gDirectory->Get(treeName);
 699   if (!fTreeD) {
 700     Warning("GetEvent","cannot find Digits Tree for event:%d\n",event);
 701   }
 702
 703
 704   // Get Reconstruct Tree header from file
 705   sprintf(treeName,"TreeR%d",event);
 706   fTreeR = (TTree*)gDirectory->Get(treeName);
 707   if (!fTreeR) {
 708     //    printf("WARNING: cannot find Reconstructed Tree for event:%d\n",event);
 709   }
 710
 711   // Set Trees branch addresses
 712   TIter next(fModules);
 713   AliModule *detector;
 714   while((detector = (AliModule*)next())) {
 715     detector->SetTreeAddress();
 716   }
 717
 718   if (fTreeK) fTreeK->GetEvent(0);
 719   fNtrack = Int_t (fParticles->GetEntries());
 720   return fNtrack;
 721 }
 722
 723 //_____________________________________________________________________________
 724 TGeometry *AliRun::GetGeometry()
 725 {
 726   //
 727   // Import Alice geometry from current file
 728   // Return pointer to geometry object
 729   //
 730   if (!fGeometry) fGeometry = (TGeometry*)gDirectory->Get("AliceGeom");
 731   //
 732   // Unlink and relink nodes in detectors
 733   // This is bad and there must be a better way...
 734   //
 735
 736   TIter next(fModules);
 737   AliModule *detector;
 738   while((detector = (AliModule*)next())) {
 739     detector->SetTreeAddress();
 740     TList *dnodes=detector->Nodes();
 741     Int_t j;
 742     TNode *node, *node1;
 743     for ( j=0; j<dnodes->GetSize(); j++) {
 744       node = (TNode*) dnodes->At(j);
 745       node1 = fGeometry->GetNode(node->GetName());
 746       dnodes->Remove(node);
 747       dnodes->AddAt(node1,j);
 748     }
 749   }
 750   return fGeometry;
 751 }
 752
 753 //_____________________________________________________________________________
 754 void AliRun::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom,
 755                           Float_t &e, Float_t *vpos, Float_t *polar,
 756                           Float_t &tof)
 757 {
 758   //
 759   // Return next track from stack of particles
 760   //
 761   TVector3 pol;
 762   fCurrent=-1;
 763   TParticle *track;
 764   for(Int_t i=fNtrack-1; i>=0; i--) {
 765     track=(TParticle*) fParticles->UncheckedAt(i);
 766     if(!track->TestBit(Done_Bit)) {
 767       //
 768       // The track has not yet been processed
 769       fCurrent=i;
 770       ipart=track->GetPdgCode();
 771       pmom[0]=track->Px();
 772       pmom[1]=track->Py();
 773       pmom[2]=track->Pz();
 774       e     =track->Energy();
 775       vpos[0]=track->Vx();
 776       vpos[1]=track->Vy();
 777       vpos[2]=track->Vz();
 778       track->GetPolarisation(pol);
 779       polar[0]=pol.X();
 780       polar[1]=pol.Y();
 781       polar[2]=pol.Z();
 782       tof=track->T();
 783       track->SetBit(Done_Bit);
 784       break;
 785     }
 786   }
 787   mtrack=fCurrent;
 788   //
 789   // stop and start timer when we start a primary track
 790   Int_t nprimaries = fHeader.GetNprimary();
 791   if (fCurrent >= nprimaries) return;
 792   if (fCurrent < nprimaries-1) {
 793     fTimer.Stop();
 794     track=(TParticle*) fParticles->UncheckedAt(fCurrent+1);
 795     //    track->SetProcessTime(fTimer.CpuTime());
 796   }
 797   fTimer.Start();
 798 }
 799
 800 //_____________________________________________________________________________
 801 Int_t AliRun::GetPrimary(Int_t track)
 802 {
 803   //
 804   // return number of primary that has generated track
 805   //
 806   int current, parent;
 807   TParticle *part;
 808   //
 809   parent=track;
 810   while (1) {
 811     current=parent;
 812     part = (TParticle *)fParticles->UncheckedAt(current);
 813     parent=part->GetFirstMother();
 814     if(parent<0) return current;
 815   }
 816 }
 817
 818 //_____________________________________________________________________________
 819 void AliRun::InitMC(const char *setup)
 820 {
 821   //
 822   // Initialize the Alice setup
 823   //
 824
 825   gROOT->LoadMacro(setup);
 826   gInterpreter->ProcessLine(fConfigFunction.Data());
 827
 828   gMC->DefineParticles();  //Create standard MC particles
 829
 830   TObject *objfirst, *objlast;
 831
 832   fNdets = fModules->GetLast()+1;
 833
 834   //
 835   //=================Create Materials and geometry
 836   gMC->Init();
 837
 838    TIter next(fModules);
 839    AliModule *detector;
 840    while((detector = (AliModule*)next())) {
 841       detector->SetTreeAddress();
 842       objlast = gDirectory->GetList()->Last();
 843
 844       // Add Detector histograms in Detector list of histograms
 845       if (objlast) objfirst = gDirectory->GetList()->After(objlast);
 846       else         objfirst = gDirectory->GetList()->First();
 847       while (objfirst) {
 848         detector->Histograms()->Add(objfirst);
 849         objfirst = gDirectory->GetList()->After(objfirst);
 850       }
 851    }
 852    SetTransPar(); //Read the cuts for all materials
 853
 854    MediaTable(); //Build the special IMEDIA table
 855
 856    //Initialise geometry deposition table
 857    fEventEnergy.Set(gMC->NofVolumes()+1);
 858    fSummEnergy.Set(gMC->NofVolumes()+1);
 859    fSum2Energy.Set(gMC->NofVolumes()+1);
 860
 861    //Compute cross-sections
 862    gMC->BuildPhysics();
 863
 864    //Write Geometry object to current file.
 865    fGeometry->Write();
 866
 867    fInitDone = kTRUE;
 868 }
 869
 870 //_____________________________________________________________________________
 871 void AliRun::MediaTable()
 872 {
 873   //
 874   // Built media table to get from the media number to
 875   // the detector id
 876   //
 877   Int_t kz, nz, idt, lz, i, k, ind;
 878   //  Int_t ibeg;
 879   TObjArray &dets = *gAlice->Detectors();
 880   AliModule *det;
 881   //
 882   // For all detectors
 883   for (kz=0;kz<fNdets;kz++) {
 884     // If detector is defined
 885     if((det=(AliModule*) dets[kz])) {
 886         TArrayI &idtmed = *(det->GetIdtmed());
 887         for(nz=0;nz<100;nz++) {
 888         // Find max and min material number
 889         if((idt=idtmed[nz])) {
 890           det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
 891           det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
 892         }
 893       }
 894       if(det->LoMedium() > det->HiMedium()) {
 895         det->LoMedium() = 0;
 896         det->HiMedium() = 0;
 897       } else {
 898         if(det->HiMedium() > fImedia->GetSize()) {
 899           Error("MediaTable","Increase fImedia from %d to %d",
 900                 fImedia->GetSize(),det->HiMedium());
 901           return;
 902         }
 903         // Tag all materials in rage as belonging to detector kz
 904         for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
 905           (*fImedia)[lz]=kz;
 906         }
 907       }
 908     }
 909   }
 910   //
 911   // Print summary table
 912   printf(" Traking media ranges:\n");
 913   for(i=0;i<(fNdets-1)/6+1;i++) {
 914     for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
 915       ind=i*6+k;
 916       det=(AliModule*)dets[ind];
 917       if(det)
 918         printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
 919                det->HiMedium());
 920       else
 921         printf(" %6s: %3d -> %3d;","NULL",0,0);
 922     }
 923     printf("\n");
 924   }
 925 }
 926
 927 //____________________________________________________________________________
 928 void AliRun::SetGenerator(AliGenerator *generator)
 929 {
 930   //
 931   // Load the event generator
 932   //
 933   if(!fGenerator) fGenerator = generator;
 934 }
 935
 936 //____________________________________________________________________________
 937 void AliRun::ResetGenerator(AliGenerator *generator)
 938 {
 939   //
 940   // Load the event generator
 941   //
 942   if(fGenerator)
 943     if(generator)
 944       Warning("ResetGenerator","Replacing generator %s with %s\n",
 945               fGenerator->GetName(),generator->GetName());
 946     else
 947       Warning("ResetGenerator","Replacing generator %s with NULL\n",
 948               fGenerator->GetName());
 949   fGenerator = generator;
 950 }
 951
 952 //____________________________________________________________________________
 953 void AliRun::SetTransPar(char* filename)
 954 {
 955   //
 956   // Read filename to set the transport parameters
 957   //
 958
 959
 960   const Int_t ncuts=10;
 961   const Int_t nflags=11;
 962   const Int_t npars=ncuts+nflags;
 963   const char pars[npars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
 964                                "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
 965                                "BREM","COMP","DCAY","DRAY","HADR","LOSS",
 966                                "MULS","PAIR","PHOT","RAYL"};
 967   char line[256];
 968   char detName[7];
 969   char* filtmp;
 970   Float_t cut[ncuts];
 971   Int_t flag[nflags];
 972   Int_t i, itmed, iret, ktmed, kz;
 973   FILE *lun;
 974   //
 975   // See whether the file is there
 976   filtmp=gSystem->ExpandPathName(filename);
 977   lun=fopen(filtmp,"r");
 978   delete [] filtmp;
 979   if(!lun) {
 980     Warning("SetTransPar","File %s does not exist!\n",filename);
 981     return;
 982   }
 983   //
 984   printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
 985   printf(" *%59s\n","*");
 986   printf(" *       Please check carefully what you are doing!%10s\n","*");
 987   printf(" *%59s\n","*");
 988   //
 989   while(1) {
 990     // Initialise cuts and flags
 991     for(i=0;i<ncuts;i++) cut[i]=-99;
 992     for(i=0;i<nflags;i++) flag[i]=-99;
 993     itmed=0;
 994     for(i=0;i<256;i++) line[i]='\0';
 995     // Read up to the end of line excluded
 996     iret=fscanf(lun,"%[^\n]",line);
 997     if(iret<0) {
 998       //End of file
 999       fclose(lun);
1000       printf(" *%59s\n","*");
1001       printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1002       return;
1003     }
1004     // Read the end of line
1005     fscanf(lun,"%*c");
1006     if(!iret) continue;
1007     if(line[0]=='*') continue;
1008     // Read the numbers
1009     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",
1010                 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
1011                 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
1012                 &flag[8],&flag[9],&flag[10]);
1013     if(!iret) continue;
1014     if(iret<0) {
1015       //reading error
1016       Warning("SetTransPar","Error reading file %s\n",filename);
1017       continue;
1018     }
1019     // Check that the module exist
1020     AliModule *mod = GetModule(detName);
1021     if(mod) {
1022       // Get the array of media numbers
1023       TArrayI &idtmed = *mod->GetIdtmed();
1024       // Check that the tracking medium code is valid
1025       if(0<=itmed && itmed < 100) {
1026         ktmed=idtmed[itmed];
1027         if(!ktmed) {
1028           Warning("SetTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
1029           continue;
1030         }
1031         // Set energy thresholds
1032         for(kz=0;kz<ncuts;kz++) {
1033           if(cut[kz]>=0) {
1034             printf(" *  %-6s set to %10.3E for tracking medium code %4d for %s\n",
1035                    pars[kz],cut[kz],itmed,mod->GetName());
1036             gMC->Gstpar(ktmed,pars[kz],cut[kz]);
1037           }
1038         }
1039         // Set transport mechanisms
1040         for(kz=0;kz<nflags;kz++) {
1041           if(flag[kz]>=0) {
1042             printf(" *  %-6s set to %10d for tracking medium code %4d for %s\n",
1043                    pars[ncuts+kz],flag[kz],itmed,mod->GetName());
1044             gMC->Gstpar(ktmed,pars[ncuts+kz],Float_t(flag[kz]));
1045           }
1046         }
1047       } else {
1048         Warning("SetTransPar","Invalid medium code %d *\n",itmed);
1049         continue;
1050       }
1051     } else {
1052       Warning("SetTransPar","Module %s not present\n",detName);
1053       continue;
1054     }
1055   }
1056 }
1057
1058 //_____________________________________________________________________________
1059 void AliRun::MakeTree(Option_t *option)
1060 {
1061   //
1062   //  Create the ROOT trees
1063   //  Loop on all detectors to create the Root branch (if any)
1064   //
1065
1066   char hname[30];
1067   //
1068   // Analyse options
1069   char *K = strstr(option,"K");
1070   char *H = strstr(option,"H");
1071   char *E = strstr(option,"E");
1072   char *D = strstr(option,"D");
1073   char *R = strstr(option,"R");
1074   //
1075   if (K && !fTreeK) {
1076     sprintf(hname,"TreeK%d",fEvent);
1077     fTreeK = new TTree(hname,"Kinematics");
1078     //  Create a branch for particles
1079     fTreeK->Branch("Particles",&fParticles,4000);
1080   }
1081   if (H && !fTreeH) {
1082     sprintf(hname,"TreeH%d",fEvent);
1083     fTreeH = new TTree(hname,"Hits");
1084     fTreeH->SetAutoSave(1000000000); //no autosave
1085   }
1086   if (D && !fTreeD) {
1087     sprintf(hname,"TreeD%d",fEvent);
1088     fTreeD = new TTree(hname,"Digits");
1089   }
1090   if (R && !fTreeR) {
1091     sprintf(hname,"TreeR%d",fEvent);
1092     fTreeR = new TTree(hname,"Reconstruction");
1093   }
1094   if (E && !fTreeE) {
1095     fTreeE = new TTree("TE","Header");
1096     //  Create a branch for Header
1097     fTreeE->Branch("Header","AliHeader",&header,4000);
1098   }
1099   //
1100   // Create a branch for hits/digits for each detector
1101   // Each branch is a TClonesArray. Each data member of the Hits classes
1102   // will be in turn a subbranch of the detector master branch
1103   TIter next(fModules);
1104   AliModule *detector;
1105   while((detector = (AliModule*)next())) {
1106      if (H || D || R) detector->MakeBranch(option);
1107   }
1108 }
1109
1110 //_____________________________________________________________________________
1111 Int_t AliRun::PurifyKine(Int_t lastSavedTrack, Int_t nofTracks)
1112 {
1113   //
1114   // PurifyKine with external parameters
1115   //
1116   fHgwmk = lastSavedTrack;
1117   fNtrack = nofTracks;
1118   PurifyKine();
1119   return fHgwmk;
1120 }
1121
1122 //_____________________________________________________________________________
1123 void AliRun::PurifyKine()
1124 {
1125   //
1126   // Compress kinematic tree keeping only flagged particles
1127   // and renaming the particle id's in all the hits
1128   //
1129   TClonesArray &particles = *fParticles;
1130   int nkeep=fHgwmk+1, parent, i;
1131   TParticle *part, *partnew, *father;
1132   int *map = new int[particles.GetEntries()];
1133
1134   // Save in Header total number of tracks before compression
1135   fHeader.SetNtrack(fHeader.GetNtrack()+fNtrack-fHgwmk);
1136
1137   // Preset map, to be removed later
1138   for(i=0; i<fNtrack; i++) {
1139     if(i<=fHgwmk) map[i]=i ; else map[i] = -99 ;}
1140   // Second pass, build map between old and new numbering
1141   for(i=fHgwmk+1; i<fNtrack; i++) {
1142     part = (TParticle *)particles.UncheckedAt(i);
1143     if(part->TestBit(Keep_Bit)) {
1144
1145       // This particle has to be kept
1146       map[i]=nkeep;
1147       if(i!=nkeep) {
1148
1149         // Old and new are different, have to copy
1150         partnew = (TParticle *)particles.UncheckedAt(nkeep);
1151         // Change due to a bug in the HP compiler
1152         //      *partnew = *part;
1153         memcpy(partnew,part,sizeof(TParticle));
1154       } else partnew = part;
1155
1156       // as the parent is always *before*, it must be already
1157       // in place. This is what we are checking anyway!
1158       if((parent=partnew->GetFirstMother())>fHgwmk) {
1159         if(map[parent]==-99) printf("map[%d] = -99!\n",parent);
1160         partnew->SetFirstMother(map[parent]);
1161       }
1162       nkeep++;
1163     }
1164   }
1165   fNtrack=nkeep;
1166
1167   // Fix daughters information
1168   for (i=0; i<fNtrack; i++) {
1169     part = (TParticle *)particles.UncheckedAt(i);
1170     parent = part->GetFirstMother();
1171     if(parent>=0) {
1172       father = (TParticle *)particles.UncheckedAt(parent);
1173       if(father->TestBit(Daughters_Bit)) {
1174
1175         if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i);
1176         if(i>father->GetLastDaughter())  father->SetLastDaughter(i);
1177       } else {
1178         // Iitialise daughters info for first pass
1179         father->SetFirstDaughter(i);
1180         father->SetLastDaughter(i);
1181         father->SetBit(Daughters_Bit);
1182       }
1183     }
1184   }
1185
1186 #ifdef old
1187   // Now loop on all detectors and reset the hits
1188   AliHit *OneHit;
1189   TIter next(fModules);
1190   AliModule *detector;
1191   while((detector = (AliModule*)next())) {
1192     if (!detector->Hits()) continue;
1193     TClonesArray &vHits=*(detector->Hits());
1194     if(vHits.GetEntries() != detector->GetNhits())
1195       printf("vHits.GetEntries()!=detector->GetNhits(): %d != %d\n",
1196              vHits.GetEntries(),detector->GetNhits());
1197     for (i=0; i<detector->GetNhits(); i++) {
1198       OneHit = (AliHit *)vHits.UncheckedAt(i);
1199       OneHit->SetTrack(map[OneHit->GetTrack()]);
1200     }
1201   }
1202 #else
1203
1204   // Now loop on all registered hit lists
1205   TIter next(fHitLists);
1206   TCollection *hitList;
1207   while((hitList = (TCollection*)next())) {
1208     TIter nexthit(hitList);
1209     AliHit *hit;
1210     while((hit = (AliHit*)nexthit())) {
1211       hit->SetTrack(map[hit->GetTrack()]);
1212     }
1213   }
1214 #endif
1215
1216   fHgwmk=nkeep-1;
1217   particles.SetLast(fHgwmk);
1218   delete [] map;
1219 }
1220
1221 //_____________________________________________________________________________
1222 void AliRun::BeginEvent()
1223 {
1224   //
1225   //  Reset all Detectors & kinematics & trees
1226   //
1227   char hname[30];
1228   //
1229
1230   //
1231   if(fLego) {
1232     fLego->BeginEvent();
1233     return;
1234   }
1235
1236   //
1237   ResetStack();
1238   ResetHits();
1239   ResetDigits();
1240
1241   // Initialise event header
1242   fHeader.Reset(fRun,fEvent);
1243
1244   if(fTreeK) {
1245     fTreeK->Reset();
1246     sprintf(hname,"TreeK%d",fEvent);
1247     fTreeK->SetName(hname);
1248   }
1249   if(fTreeH) {
1250     fTreeH->Reset();
1251     sprintf(hname,"TreeH%d",fEvent);
1252     fTreeH->SetName(hname);
1253   }
1254   if(fTreeD) {
1255     fTreeD->Reset();
1256     sprintf(hname,"TreeD%d",fEvent);
1257     fTreeD->SetName(hname);
1258   }
1259   if(fTreeR) {
1260     fTreeR->Reset();
1261     sprintf(hname,"TreeR%d",fEvent);
1262     fTreeR->SetName(hname);
1263   }
1264 }
1265
1266 //_____________________________________________________________________________
1267 void AliRun::ResetDigits()
1268 {
1269   //
1270   //  Reset all Detectors digits
1271   //
1272   TIter next(fModules);
1273   AliModule *detector;
1274   while((detector = (AliModule*)next())) {
1275      detector->ResetDigits();
1276   }
1277 }
1278
1279 //_____________________________________________________________________________
1280 void AliRun::ResetHits()
1281 {
1282   //
1283   //  Reset all Detectors hits
1284   //
1285   TIter next(fModules);
1286   AliModule *detector;
1287   while((detector = (AliModule*)next())) {
1288      detector->ResetHits();
1289   }
1290 }
1291
1292 //_____________________________________________________________________________
1293 void AliRun::ResetPoints()
1294 {
1295   //
1296   // Reset all Detectors points
1297   //
1298   TIter next(fModules);
1299   AliModule *detector;
1300   while((detector = (AliModule*)next())) {
1301      detector->ResetPoints();
1302   }
1303 }
1304
1305 //_____________________________________________________________________________
1306 void AliRun::RunMC(Int_t nevent, const char *setup)
1307 {
1308   //
1309   // Main function to be called to process a galice run
1310   // example
1311   //    Root > gAlice.Run();
1312   // a positive number of events will cause the finish routine
1313   // to be called
1314   //
1315
1316   // check if initialisation has been done
1317   if (!fInitDone) InitMC(setup);
1318
1319   // Create the Root Tree with one branch per detector
1320   MakeTree("KHDER");
1321
1322   gMC->ProcessRun(nevent);
1323
1324   // End of this run, close files
1325   if(nevent>0) FinishRun();
1326 }
1327
1328 //_____________________________________________________________________________
1329 void AliRun::RunLego(const char *setup,Int_t ntheta,Float_t themin,
1330                      Float_t themax,Int_t nphi,Float_t phimin,Float_t phimax,
1331                      Float_t rmin,Float_t rmax,Float_t zmax)
1332 {
1333   //
1334   // Generates lego plots of:
1335   //    - radiation length map phi vs theta
1336   //    - radiation length map phi vs eta
1337   //    - interaction length map
1338   //    - g/cm2 length map
1339   //
1340   //  ntheta    bins in theta, eta
1341   //  themin    minimum angle in theta (degrees)
1342   //  themax    maximum angle in theta (degrees)
1343   //  nphi      bins in phi
1344   //  phimin    minimum angle in phi (degrees)
1345   //  phimax    maximum angle in phi (degrees)
1346   //  rmin      minimum radius
1347   //  rmax      maximum radius
1348   //
1349   //
1350   //  The number of events generated = ntheta*nphi
1351   //  run input parameters in macro setup (default="Config.C")
1352   //
1353   //  Use macro "lego.C" to visualize the 3 lego plots in spherical coordinates
1354   //Begin_Html
1355   /*
1356     <img src="picts/AliRunLego1.gif">
1357   */
1358   //End_Html
1359   //Begin_Html
1360   /*
1361     <img src="picts/AliRunLego2.gif">
1362   */
1363   //End_Html
1364   //Begin_Html
1365   /*
1366     <img src="picts/AliRunLego3.gif">
1367   */
1368   //End_Html
1369   //
1370
1371   // check if initialisation has been done
1372   if (!fInitDone) InitMC(setup);
1373
1374   //Save current generator
1375   AliGenerator *gen=Generator();
1376
1377   //Create Lego object
1378   fLego = new AliLego("lego",ntheta,themin,themax,nphi,phimin,phimax,rmin,rmax,zmax);
1379
1380   //Prepare MC for Lego Run
1381   gMC->InitLego();
1382
1383   //Run Lego Object
1384   gMC->ProcessRun(ntheta*nphi+1);
1385
1386   // Create only the Root event Tree
1387   MakeTree("E");
1388
1389   // End of this run, close files
1390   FinishRun();
1391
1392   // Delete Lego Object
1393   delete fLego; fLego=0;
1394
1395   // Restore current generator
1396   SetGenerator(gen);
1397 }
1398
1399 //_____________________________________________________________________________
1400 void AliRun::SetCurrentTrack(Int_t track)
1401 {
1402   //
1403   // Set current track number
1404   //
1405   fCurrent = track;
1406 }
1407
1408 //_____________________________________________________________________________
1409 void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
1410                       Float_t *vpos, Float_t *polar, Float_t tof,
1411                       const char *mecha, Int_t &ntr, Float_t weight)
1412 {
1413   //
1414   // Load a track on the stack
1415   //
1416   // done     0 if the track has to be transported
1417   //          1 if not
1418   // parent   identifier of the parent track. -1 for a primary
1419   // pdg    particle code
1420   // pmom     momentum GeV/c
1421   // vpos     position
1422   // polar    polarisation
1423   // tof      time of flight in seconds
1424   // mecha    production mechanism
1425   // ntr      on output the number of the track stored
1426   //
1427   TClonesArray &particles = *fParticles;
1428   TParticle *particle;
1429   Float_t mass;
1430   const Int_t firstdaughter=-1;
1431   const Int_t lastdaughter=-1;
1432   const Int_t KS=0;
1433   //  const Float_t tlife=0;
1434
1435   //
1436   // Here we get the static mass
1437   // For MC is ok, but a more sophisticated method could be necessary
1438   // if the calculated mass is required
1439   // also, this method is potentially dangerous if the mass
1440   // used in the MC is not the same of the PDG database
1441   //
1442   mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
1443   Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
1444                         pmom[1]*pmom[1]+pmom[2]*pmom[2]);
1445
1446   //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",
1447   //pname,mass,e,fNtrack,pdg,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],KS,mecha);
1448
1449   particle=new(particles[fNtrack]) TParticle(pdg,KS,parent,-1,firstdaughter,
1450                                              lastdaughter,pmom[0],pmom[1],pmom[2],
1451                                              e,vpos[0],vpos[1],vpos[2],tof);
1452   //                                         polar[0],polar[1],polar[2],tof,
1453   //                                         mecha,weight);
1454   ((TParticle*)particles[fNtrack])->SetPolarisation(TVector3(polar[0],polar[1],polar[2]));
1455   ((TParticle*)particles[fNtrack])->SetWeight(weight);
1456   if(!done) particle->SetBit(Done_Bit);
1457
1458   if(parent>=0) {
1459     particle=(TParticle*) fParticles->UncheckedAt(parent);
1460     particle->SetLastDaughter(fNtrack);
1461     if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
1462   } else {
1463     //
1464     // This is a primary track. Set high water mark for this event
1465     fHgwmk=fNtrack;
1466     //
1467     // Set also number if primary tracks
1468     fHeader.SetNprimary(fHgwmk+1);
1469     fHeader.SetNtrack(fHgwmk+1);
1470   }
1471   ntr = fNtrack++;
1472 }
1473
1474 //_____________________________________________________________________________
1475 void AliRun::KeepTrack(const Int_t track)
1476 {
1477   //
1478   // flags a track to be kept
1479   //
1480   TClonesArray &particles = *fParticles;
1481   ((TParticle*)particles[track])->SetBit(Keep_Bit);
1482 }
1483
1484 //_____________________________________________________________________________
1485 void AliRun::StepManager(Int_t id)
1486 {
1487   //
1488   // Called at every step during transport
1489   //
1490
1491   //
1492   // --- If lego option, do it and leave
1493   if (fLego)
1494     fLego->StepManager();
1495   else {
1496     Int_t copy;
1497     //Update energy deposition tables
1498     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
1499
1500     //Call the appropriate stepping routine;
1501     AliModule *det = (AliModule*)fModules->At(id);
1502     if(det) det->StepManager();
1503   }
1504 }
1505
1506 //_____________________________________________________________________________
1507 void AliRun::Streamer(TBuffer &R__b)
1508 {
1509   //
1510   // Stream an object of class AliRun.
1511   //
1512   if (R__b.IsReading()) {
1513     Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1514     TNamed::Streamer(R__b);
1515     if (!gAlice) gAlice = this;
1516     gROOT->GetListOfBrowsables()->Add(this,"Run");
1517     fTreeE = (TTree*)gDirectory->Get("TE");
1518     if (fTreeE) fTreeE->SetBranchAddress("Header", &header);
1519     else    Error("Streamer","cannot find Header Tree\n");
1520     R__b >> fNtrack;
1521     R__b >> fHgwmk;
1522     R__b >> fDebug;
1523     fHeader.Streamer(R__b);
1524     R__b >> fModules;
1525     R__b >> fParticles;
1526     R__b >> fField;
1527     //    R__b >> fMC;
1528     R__b >> fNdets;
1529     R__b >> fTrRmax;
1530     R__b >> fTrZmax;
1531     R__b >> fGenerator;
1532     if(R__v>1) {
1533       R__b >> fPDGDB;        //Particle factory object!
1534       fTreeE->GetEntry(0);
1535     } else {
1536       fHeader.SetEvent(0);
1537       fPDGDB     = TDatabasePDG::Instance();        //Particle factory object!
1538     }
1539     if(R__v>2) {
1540       fConfigFunction.Streamer(R__b);
1541     } else {
1542       fConfigFunction="Config();";
1543     }
1544   } else {
1545     R__b.WriteVersion(AliRun::IsA());
1546     TNamed::Streamer(R__b);
1547     R__b << fNtrack;
1548     R__b << fHgwmk;
1549     R__b << fDebug;
1550     fHeader.Streamer(R__b);
1551     R__b << fModules;
1552     R__b << fParticles;
1553     R__b << fField;
1554     //    R__b << fMC;
1555     R__b << fNdets;
1556     R__b << fTrRmax;
1557     R__b << fTrZmax;
1558     R__b << fGenerator;
1559     R__b << fPDGDB;        //Particle factory object!
1560     fConfigFunction.Streamer(R__b);
1561   }
1562 }