STEER/STEER/AliMC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // This class is extracted from the AliRun class
  19 // and contains all the MC-related functionality
  20 // The number of dependencies has to be reduced...
  21 // Author: F.Carminati
  22 //         Federico.Carminati@cern.ch
  23
  24 #include <string.h>
  25
  26 #include <RVersion.h>
  27 #include <TArrayI.h>
  28 #include <TClonesArray.h>
  29 #include <TFile.h>
  30 #include <TGeoGlobalMagField.h>
  31 #include <TGeoManager.h>
  32 #include <TParticle.h>
  33 #include <TROOT.h>
  34 #include <TStopwatch.h>
  35 #include <TSystem.h>
  36 #include <TVirtualMC.h>
  37 #include <TTree.h>
  38
  39 #include "AliCDBEntry.h"
  40 #include "AliCDBManager.h"
  41 #include "AliCDBStorage.h"
  42 #include "AliDetector.h"
  43 #include "AliGenerator.h"
  44 #include "AliGeomManager.h"
  45 #include "AliHeader.h"
  46 #include "AliHit.h"
  47 #include "AliLego.h"
  48 #include "AliLog.h"
  49 #include "AliMC.h"
  50 #include "AliMagF.h"
  51 #include "AliRun.h"
  52 #include "AliSimulation.h"
  53 #include "AliStack.h"
  54 #include "AliTrackReference.h"
  55
  56 ClassImp(AliMC)
  57
  58 //_______________________________________________________________________
  59 AliMC::AliMC() :
  60   fGenerator(0),
  61   fSaveRndmStatus(kFALSE),
  62   fSaveRndmEventStatus(kFALSE),
  63   fReadRndmStatus(kFALSE),
  64   fRndmFileName("random.root"),
  65   fEventEnergy(0),
  66   fSummEnergy(0),
  67   fSum2Energy(0),
  68   fTrRmax(1.e10),
  69   fTrZmax(1.e10),
  70   fRDecayMax(1.e10),
  71   fRDecayMin(-1.),
  72   fDecayPdg(0),
  73   fImedia(0),
  74   fTransParName("\0"),
  75   fHitLists(0),
  76   fTmpTreeTR(0),
  77   fTmpFileTR(0),
  78   fTrackReferences(),
  79   fTmpTrackReferences()
  80
  81 {
  82   //default constructor
  83   DecayLimits();
  84 }
  85
  86 //_______________________________________________________________________
  87 AliMC::AliMC(const char *name, const char *title) :
  88   TVirtualMCApplication(name, title),
  89   fGenerator(0),
  90   fSaveRndmStatus(kFALSE),
  91   fSaveRndmEventStatus(kFALSE),
  92   fReadRndmStatus(kFALSE),
  93   fRndmFileName("random.root"),
  94   fEventEnergy(0),
  95   fSummEnergy(0),
  96   fSum2Energy(0),
  97   fTrRmax(1.e10),
  98   fTrZmax(1.e10),
  99   fRDecayMax(1.e10),
 100   fRDecayMin(-1.),
 101   fDecayPdg(0),
 102   fImedia(new TArrayI(1000)),
 103   fTransParName("\0"),
 104   fHitLists(new TList()),
 105   fTmpTreeTR(0),
 106   fTmpFileTR(0),
 107   fTrackReferences("AliTrackReference", 100),
 108   fTmpTrackReferences("AliTrackReference", 100)
 109 {
 110   //constructor
 111   // Set transport parameters
 112   SetTransPar();
 113   DecayLimits();
 114   // Prepare the tracking medium lists
 115   for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
 116 }
 117
 118 //_______________________________________________________________________
 119 AliMC::~AliMC()
 120 {
 121   //destructor
 122   delete fGenerator;
 123   delete fImedia;
 124   delete fHitLists;
 125   // Delete track references
 126 }
 127
 128 //_______________________________________________________________________
 129 void  AliMC::ConstructGeometry()
 130 {
 131   //
 132   // Either load geometry from file or create it through usual
 133   // loop on detectors. In the first case the method
 134   // AliModule::CreateMaterials() only builds fIdtmed and is postponed
 135   // at InitGeometry().
 136   //
 137
 138   if(AliSimulation::Instance()->IsGeometryFromFile()){ //load geometry either from CDB or from file
 139     if(IsGeometryFromCDB()){
 140       AliInfo("Loading geometry from CDB default storage");
 141       AliCDBPath path("GRP","Geometry","Data");
 142       AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
 143       if(!entry) AliFatal("Unable to load geometry from CDB!");
 144       entry->SetOwner(0);
 145       gGeoManager = (TGeoManager*) entry->GetObject();
 146       if (!gGeoManager) AliFatal("TGeoManager object not found in the specified CDB entry!");
 147     }else{
 148       // Load geometry
 149       const char *geomfilename = AliSimulation::Instance()->GetGeometryFile();
 150       if(gSystem->ExpandPathName(geomfilename)){
 151         AliInfo(Form("Loading geometry from file:\n %40s",geomfilename));
 152         TGeoManager::Import(geomfilename);
 153       }else{
 154         AliInfo(Form("Geometry file %40s not found!\n",geomfilename));
 155         return;
 156       }
 157     }
 158     gMC->SetRootGeometry();
 159   }else{
 160     // Create modules, materials, geometry
 161     if (!gGeoManager) new TGeoManager("ALICE", "ALICE geometry");
 162     TStopwatch stw;
 163     TIter next(gAlice->Modules());
 164     AliModule *detector;
 165     AliDebug(1, "Geometry creation:");
 166     while((detector = dynamic_cast<AliModule*>(next()))) {
 167       stw.Start();
 168       // Initialise detector materials and geometry
 169       detector->CreateMaterials();
 170       detector->CreateGeometry();
 171       AliInfo(Form("%10s R:%.2fs C:%.2fs",
 172                    detector->GetName(),stw.RealTime(),stw.CpuTime()));
 173     }
 174   }
 175
 176 }
 177
 178 //_______________________________________________________________________
 179 Bool_t  AliMC::MisalignGeometry()
 180 {
 181   // Call misalignment code if AliSimulation object was defined.
 182
 183   if(!AliSimulation::Instance()->IsGeometryFromFile()){
 184     //Set alignable volumes for the whole geometry
 185     SetAllAlignableVolumes();
 186   }
 187   // Misalign geometry via AliSimulation instance
 188   if (!AliSimulation::Instance()) return kFALSE;
 189   AliGeomManager::SetGeometry(gGeoManager);
 190   if(!AliGeomManager::CheckSymNamesLUT("ALL"))
 191     AliFatal("Current loaded geometry differs in the definition of symbolic names!");
 192
 193   return AliSimulation::Instance()->MisalignGeometry(AliRunLoader::Instance());
 194 }
 195
 196 //_______________________________________________________________________
 197 void  AliMC::ConstructOpGeometry()
 198 {
 199   //
 200   // Loop all detector modules and call DefineOpticalProperties() method
 201   //
 202
 203   TIter next(gAlice->Modules());
 204   AliModule *detector;
 205   AliInfo("Optical properties definition");
 206   while((detector = dynamic_cast<AliModule*>(next()))) {
 207     // Initialise detector geometry
 208     if(AliSimulation::Instance()->IsGeometryFromFile()) detector->CreateMaterials();
 209     // Initialise detector optical properties
 210     detector->DefineOpticalProperties();
 211   }
 212 }
 213
 214 //_______________________________________________________________________
 215 void  AliMC::InitGeometry()
 216 {
 217   //
 218   // Initialize detectors
 219   //
 220
 221   AliInfo("Initialisation:");
 222   TStopwatch stw;
 223   TIter next(gAlice->Modules());
 224   AliModule *detector;
 225   while((detector = dynamic_cast<AliModule*>(next()))) {
 226     stw.Start();
 227     detector->Init();
 228     AliInfo(Form("%10s R:%.2fs C:%.2fs",
 229                  detector->GetName(),stw.RealTime(),stw.CpuTime()));
 230   }
 231 }
 232
 233 //_______________________________________________________________________
 234 void AliMC::SetGeometryFromCDB()
 235 {
 236   // Set the loading of geometry from cdb instead of creating it
 237   // A default CDB storage needs to be set before this method is called
 238   if(AliCDBManager::Instance()->IsDefaultStorageSet() &&
 239      AliCDBManager::Instance()->GetRun() >= 0)
 240     AliSimulation::Instance()->SetGeometryFile("*OCDB*");
 241   else
 242     AliError("Loading of geometry from CDB ignored. First set a default CDB storage!");
 243 }
 244
 245 //_______________________________________________________________________
 246 Bool_t AliMC::IsGeometryFromCDB() const
 247 {
 248   return (strcmp(AliSimulation::Instance()->GetGeometryFile(),"*OCDB*")==0);
 249 }
 250
 251 //_______________________________________________________________________
 252 void  AliMC::SetAllAlignableVolumes()
 253 {
 254   //
 255   // Add alignable volumes (TGeoPNEntries) looping on all
 256   // active modules
 257   //
 258
 259   AliInfo(Form("Setting entries for all alignable volumes of active detectors"));
 260   AliModule *detector;
 261   TIter next(gAlice->Modules());
 262   while((detector = dynamic_cast<AliModule*>(next()))) {
 263     detector->AddAlignableVolumes();
 264   }
 265 }
 266
 267 //_______________________________________________________________________
 268 void  AliMC::GeneratePrimaries()
 269 {
 270   //
 271   // Generate primary particles and fill them in the stack.
 272   //
 273
 274   Generator()->Generate();
 275 }
 276
 277 //_______________________________________________________________________
 278 void AliMC::SetGenerator(AliGenerator *generator)
 279 {
 280   //
 281   // Load the event generator
 282   //
 283   if(!fGenerator) fGenerator = generator;
 284 }
 285
 286 //_______________________________________________________________________
 287 void AliMC::ResetGenerator(AliGenerator *generator)
 288 {
 289   //
 290   // Load the event generator
 291   //
 292   if(fGenerator) {
 293     if(generator) {
 294       AliWarning(Form("Replacing generator %s with %s",
 295                       fGenerator->GetName(),generator->GetName()));
 296     }
 297     else {
 298       AliWarning(Form("Replacing generator %s with NULL",
 299                       fGenerator->GetName()));
 300     }
 301   }
 302   fGenerator = generator;
 303 }
 304
 305 //_______________________________________________________________________
 306 void AliMC::FinishRun()
 307 {
 308   // Clean generator information
 309   AliDebug(1, "fGenerator->FinishRun()");
 310   fGenerator->FinishRun();
 311
 312   //Output energy summary tables
 313   AliDebug(1, "EnergySummary()");
 314   ToAliDebug(1, EnergySummary());
 315 }
 316
 317 //_______________________________________________________________________
 318 void AliMC::BeginPrimary()
 319 {
 320   //
 321   // Called  at the beginning of each primary track
 322   //
 323
 324   // Reset Hits info
 325   ResetHits();
 326   ResetTrackReferences();
 327 }
 328
 329 //_______________________________________________________________________
 330 void AliMC::PreTrack()
 331 {
 332   // Actions before the track's transport
 333
 334      TObjArray &dets = *gAlice->Modules();
 335      AliModule *module;
 336
 337      for(Int_t i=0; i<=gAlice->GetNdets(); i++)
 338        if((module = dynamic_cast<AliModule*>(dets[i])))
 339          module->PreTrack();
 340 }
 341
 342 //_______________________________________________________________________
 343 void AliMC::Stepping()
 344 {
 345   //
 346   // Called at every step during transport
 347   //
 348   Int_t id = DetFromMate(gMC->CurrentMedium());
 349   if (id < 0) return;
 350
 351
 352   if ( gMC->IsNewTrack()            &&
 353        gMC->TrackTime() == 0.       &&
 354        fRDecayMin >= 0.             &&
 355        fRDecayMax > fRDecayMin      &&
 356        gMC->TrackPid() == fDecayPdg )
 357   {
 358       FixParticleDecaytime();
 359   }
 360
 361
 362
 363   //
 364   // --- If lego option, do it and leave
 365   if (AliSimulation::Instance()->Lego())
 366     AliSimulation::Instance()->Lego()->StepManager();
 367   else {
 368     Int_t copy;
 369     //Update energy deposition tables
 370     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
 371     //
 372     // write tracke reference for track which is dissapearing - MI
 373
 374     if (gMC->IsTrackDisappeared() && !(gMC->IsTrackAlive())) {
 375         if (gMC->Etot() > 0.05) AddTrackReference(GetCurrentTrackNumber(),
 376                                                 AliTrackReference::kDisappeared);
 377
 378
 379     }
 380
 381     //Call the appropriate stepping routine;
 382     AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
 383     if(det && det->StepManagerIsEnabled()) {
 384       det->StepManager();
 385     }
 386   }
 387 }
 388
 389 //_______________________________________________________________________
 390 void AliMC::EnergySummary()
 391 {
 392   //e
 393   // Print summary of deposited energy
 394   //
 395
 396   Int_t ndep=0;
 397   Float_t edtot=0;
 398   Float_t ed, ed2;
 399   Int_t kn, i, left, j, id;
 400   const Float_t kzero=0;
 401   Int_t ievent=AliRunLoader::Instance()->GetHeader()->GetEvent()+1;
 402   //
 403   // Energy loss information
 404   if(ievent) {
 405     printf("***************** Energy Loss Information per event (GEV) *****************\n");
 406     for(kn=1;kn<fEventEnergy.GetSize();kn++) {
 407       ed=fSummEnergy[kn];
 408       if(ed>0) {
 409         fEventEnergy[ndep]=kn;
 410         if(ievent>1) {
 411           ed=ed/ievent;
 412           ed2=fSum2Energy[kn];
 413           ed2=ed2/ievent;
 414           ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
 415         } else
 416           ed2=99;
 417         fSummEnergy[ndep]=ed;
 418         fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
 419         edtot+=ed;
 420         ndep++;
 421       }
 422     }
 423     for(kn=0;kn<(ndep-1)/3+1;kn++) {
 424       left=ndep-kn*3;
 425       for(i=0;i<(3<left?3:left);i++) {
 426         j=kn*3+i;
 427         id=Int_t (fEventEnergy[j]+0.1);
 428         printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
 429       }
 430       printf("\n");
 431     }
 432     //
 433     // Relative energy loss in different detectors
 434     printf("******************** Relative Energy Loss per event ********************\n");
 435     printf("Total energy loss per event %10.3f GeV\n",edtot);
 436     for(kn=0;kn<(ndep-1)/5+1;kn++) {
 437       left=ndep-kn*5;
 438       for(i=0;i<(5<left?5:left);i++) {
 439         j=kn*5+i;
 440         id=Int_t (fEventEnergy[j]+0.1);
 441         printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
 442       }
 443       printf("\n");
 444     }
 445     for(kn=0;kn<75;kn++) printf("*");
 446     printf("\n");
 447   }
 448   //
 449   // Reset the TArray's
 450   //  fEventEnergy.Set(0);
 451   //  fSummEnergy.Set(0);
 452   //  fSum2Energy.Set(0);
 453 }
 454 #include <TFile.h>
 455 //_____________________________________________________________________________
 456 void AliMC::BeginEvent()
 457 {
 458   //
 459   // Clean-up previous event
 460   // Energy scores
 461   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 462   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 463   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 464   AliDebug(1, "          BEGINNING EVENT               ");
 465   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 466   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 467   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 468
 469   AliRunLoader *runloader=AliRunLoader::Instance();
 470
 471   /*******************************/
 472   /*   Clean after eventual      */
 473   /*   previous event            */
 474   /*******************************/
 475
 476
 477   //Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors),
 478   gAlice->SetEventNrInRun(gAlice->GetEventNrInRun()+1);
 479   runloader->SetEventNumber(gAlice->GetEventNrInRun());// sets new files, cleans the previous event stuff, if necessary, etc.,
 480   AliDebug(1, Form("EventNr is %d",gAlice->GetEventNrInRun()));
 481
 482   fEventEnergy.Reset();
 483     // Clean detector information
 484
 485   if (runloader->Stack())
 486       runloader->Stack()->Reset();//clean stack -> tree is unloaded
 487   else
 488       runloader->MakeStack();//or make a new one
 489
 490   // Random engine status
 491   //
 492
 493   if ( fSaveRndmStatus || fSaveRndmEventStatus) {
 494     TString fileName="random";
 495     if ( fSaveRndmEventStatus ) {
 496       fileName += "Evt";
 497       fileName += gAlice->GetEventNrInRun();
 498     }
 499     fileName += ".root";
 500
 501     // write ROOT random engine status
 502     cout << "Saving random engine status in " << fileName.Data() << endl;
 503     TFile f(fileName.Data(),"RECREATE");
 504     gRandom->Write(fileName.Data());
 505   }
 506
 507   if ( fReadRndmStatus ) {
 508     //read ROOT random engine status
 509     cout << "Reading random engine status from " << fRndmFileName.Data() << endl;
 510     TFile f(fRndmFileName.Data());
 511     gRandom->Read(fRndmFileName.Data());
 512   }
 513
 514   if(AliSimulation::Instance()->Lego() == 0x0)
 515   {
 516       AliDebug(1, "fRunLoader->MakeTree(K)");
 517       runloader->MakeTree("K");
 518   }
 519
 520   AliDebug(1, "gMC->SetStack(fRunLoader->Stack())");
 521   gMC->SetStack(runloader->Stack());//Was in InitMC - but was moved here
 522                                      //because we don't have guarantee that
 523                                      //stack pointer is not going to change from event to event
 524                          //since it bellobgs to header and is obtained via RunLoader
 525   //
 526   //  Reset all Detectors & kinematics & make/reset trees
 527   //
 528
 529   runloader->GetHeader()->Reset(AliCDBManager::Instance()->GetRun(),gAlice->GetEvNumber(),
 530                                 gAlice->GetEventNrInRun());
 531 //  fRunLoader->WriteKinematics("OVERWRITE");  is there any reason to rewrite here since MakeTree does so
 532
 533   if(AliSimulation::Instance()->Lego())
 534   {
 535       AliSimulation::Instance()->Lego()->BeginEvent();
 536       return;
 537   }
 538
 539
 540   AliDebug(1, "ResetHits()");
 541   ResetHits();
 542
 543   AliDebug(1, "fRunLoader->MakeTree(H)");
 544   runloader->MakeTree("H");
 545
 546
 547
 548   MakeTmpTrackRefsTree();
 549   //create new branches and SetAdresses
 550   TIter next(gAlice->Modules());
 551   AliModule *detector;
 552   while((detector = (AliModule*)next()))
 553    {
 554        AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName()));
 555        detector->MakeBranch("H");
 556    }
 557 }
 558
 559 //_______________________________________________________________________
 560 void AliMC::ResetHits()
 561 {
 562   //
 563   //  Reset all Detectors hits
 564   //
 565   TIter next(gAlice->Modules());
 566   AliModule *detector;
 567   while((detector = dynamic_cast<AliModule*>(next()))) {
 568      detector->ResetHits();
 569   }
 570 }
 571
 572 //_______________________________________________________________________
 573 void AliMC::ResetDigits()
 574 {
 575   //
 576   //  Reset all Detectors digits
 577   //
 578   TIter next(gAlice->Modules());
 579   AliModule *detector;
 580   while((detector = dynamic_cast<AliModule*>(next()))) {
 581      detector->ResetDigits();
 582   }
 583 }
 584
 585 //_______________________________________________________________________
 586 void AliMC::ResetSDigits()
 587 {
 588   //
 589   //  Reset all Detectors digits
 590   //
 591   TIter next(gAlice->Modules());
 592   AliModule *detector;
 593   while((detector = dynamic_cast<AliModule*>(next()))) {
 594      detector->ResetSDigits();
 595   }
 596 }
 597
 598 //_______________________________________________________________________
 599 void AliMC::PostTrack()
 600 {
 601   // Posts tracks for each module
 602
 603   TObjArray &dets = *gAlice->Modules();
 604   AliModule *module;
 605
 606   for(Int_t i=0; i<=gAlice->GetNdets(); i++)
 607     if((module = dynamic_cast<AliModule*>(dets[i])))
 608       module->PostTrack();
 609 }
 610
 611 //_______________________________________________________________________
 612 void AliMC::FinishPrimary()
 613 {
 614   //
 615   // Called  at the end of each primary track
 616   //
 617
 618   AliRunLoader *runloader=AliRunLoader::Instance();
 619   //  static Int_t count=0;
 620   //  const Int_t times=10;
 621   // This primary is finished, purify stack
 622 #if ROOT_VERSION_CODE > 262152
 623   if (!(gMC->SecondariesAreOrdered())) {
 624       if (runloader->Stack()->ReorderKine()) RemapHits();
 625   }
 626 #endif
 627   if (runloader->Stack()->PurifyKine()) RemapHits();
 628
 629   TIter next(gAlice->Modules());
 630   AliModule *detector;
 631   while((detector = dynamic_cast<AliModule*>(next()))) {
 632     detector->FinishPrimary();
 633     AliLoader* loader = detector->GetLoader();
 634     if(loader)
 635      {
 636        TTree* treeH = loader->TreeH();
 637        if (treeH) treeH->Fill(); //can be Lego run and treeH can not exist
 638      }
 639   }
 640
 641   // Write out track references if any
 642   if (fTmpTreeTR) fTmpTreeTR->Fill();
 643 }
 644
 645 void AliMC::RemapHits()
 646 {
 647 //
 648 // Remaps the track labels of the hits
 649     AliRunLoader *runloader=AliRunLoader::Instance();
 650     AliStack* stack = runloader->Stack();
 651     TList* hitLists = GetHitLists();
 652     TIter next(hitLists);
 653     TCollection *hitList;
 654
 655     while((hitList = dynamic_cast<TCollection*>(next()))) {
 656         TIter nexthit(hitList);
 657         AliHit *hit;
 658         while((hit = dynamic_cast<AliHit*>(nexthit()))) {
 659             hit->SetTrack(stack->TrackLabel(hit->GetTrack()));
 660         }
 661     }
 662
 663     //
 664     // This for detectors which have a special mapping mechanism
 665     // for hits, such as TPC and TRD
 666     //
 667
 668
 669     TObjArray* modules = gAlice->Modules();
 670     TIter nextmod(modules);
 671     AliModule *module;
 672     while((module = (AliModule*) nextmod())) {
 673         AliDetector* det = dynamic_cast<AliDetector*> (module);
 674         if (det) det->RemapTrackHitIDs(stack->TrackLabelMap());
 675     }
 676     //
 677     RemapTrackReferencesIDs(stack->TrackLabelMap());
 678 }
 679
 680 //_______________________________________________________________________
 681 void AliMC::FinishEvent()
 682 {
 683   //
 684   // Called at the end of the event.
 685   //
 686
 687   if(AliSimulation::Instance()->Lego()) AliSimulation::Instance()->Lego()->FinishEvent();
 688
 689   TIter next(gAlice->Modules());
 690   AliModule *detector;
 691   while((detector = dynamic_cast<AliModule*>(next()))) {
 692     detector->FinishEvent();
 693   }
 694
 695   //Update the energy deposit tables
 696   Int_t i;
 697   for(i=0;i<fEventEnergy.GetSize();i++)
 698    {
 699     fSummEnergy[i]+=fEventEnergy[i];
 700     fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
 701    }
 702
 703   AliRunLoader *runloader=AliRunLoader::Instance();
 704
 705   AliHeader* header = runloader->GetHeader();
 706   AliStack* stack = runloader->Stack();
 707   if ( (header == 0x0) || (stack == 0x0) )
 708    {//check if we got header and stack. If not cry and exit aliroot
 709     AliFatal("Can not get the stack or header from LOADER");
 710     return;//never reached
 711    }
 712   // Update Header information
 713   header->SetNprimary(stack->GetNprimary());
 714   header->SetNtrack(stack->GetNtrack());
 715   header->SetTimeStamp(AliSimulation::Instance()->GenerateTimeStamp());
 716
 717   // Write out the kinematics
 718   if (!AliSimulation::Instance()->Lego()) stack->FinishEvent();
 719
 720   // Synchronize the TreeTR with TreeK
 721   if (fTmpTreeTR) ReorderAndExpandTreeTR();
 722
 723   // Write out the event Header information
 724   TTree* treeE = runloader->TreeE();
 725   if (treeE)
 726    {
 727       header->SetStack(stack);
 728       treeE->Fill();
 729    }
 730   else
 731    {
 732     AliError("Can not get TreeE from RL");
 733    }
 734
 735   if(AliSimulation::Instance()->Lego() == 0x0)
 736    {
 737      runloader->WriteKinematics("OVERWRITE");
 738      runloader->WriteTrackRefs("OVERWRITE");
 739      runloader->WriteHits("OVERWRITE");
 740    }
 741
 742   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 743   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 744   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 745   AliDebug(1, "          FINISHING EVENT               ");
 746   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 747   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 748   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 749 }
 750
 751 //_______________________________________________________________________
 752 void AliMC::Init()
 753 {
 754   // MC initialization
 755
 756
 757    //=================Create Materials and geometry
 758    gMC->Init();
 759   // Set alignable volumes for the whole geometry (with old root)
 760 #if ROOT_VERSION_CODE < 331527
 761   SetAllAlignableVolumes();
 762 #endif
 763    //Read the cuts for all materials
 764    ReadTransPar();
 765    //Build the special IMEDIA table
 766    MediaTable();
 767
 768    //Compute cross-sections
 769    gMC->BuildPhysics();
 770
 771    //Initialise geometry deposition table
 772    fEventEnergy.Set(gMC->NofVolumes()+1);
 773    fSummEnergy.Set(gMC->NofVolumes()+1);
 774    fSum2Energy.Set(gMC->NofVolumes()+1);
 775
 776    // Register MC in configuration
 777    AliConfig::Instance()->Add(gMC);
 778 }
 779
 780 //_______________________________________________________________________
 781 void AliMC::MediaTable()
 782 {
 783   //
 784   // Built media table to get from the media number to
 785   // the detector id
 786   //
 787
 788   Int_t kz, nz, idt, lz, i, k, ind;
 789   //  Int_t ibeg;
 790   TObjArray &dets = *gAlice->Detectors();
 791   AliModule *det;
 792   Int_t ndets=gAlice->GetNdets();
 793   //
 794   // For all detectors
 795   for (kz=0;kz<ndets;kz++) {
 796     // If detector is defined
 797     if((det=dynamic_cast<AliModule*>(dets[kz]))) {
 798         TArrayI &idtmed = *(det->GetIdtmed());
 799         for(nz=0;nz<100;nz++) {
 800
 801         // Find max and min material number
 802         if((idt=idtmed[nz])) {
 803           det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
 804           det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
 805         }
 806       }
 807       if(det->LoMedium() > det->HiMedium()) {
 808         det->LoMedium() = 0;
 809         det->HiMedium() = 0;
 810       } else {
 811         if(det->HiMedium() > fImedia->GetSize()) {
 812           AliError(Form("Increase fImedia from %d to %d",
 813                         fImedia->GetSize(),det->HiMedium()));
 814           return;
 815         }
 816         // Tag all materials in rage as belonging to detector kz
 817         for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
 818           (*fImedia)[lz]=kz;
 819         }
 820       }
 821     }
 822   }
 823   //
 824   // Print summary table
 825   AliInfo("Tracking media ranges:");
 826   ToAliInfo(
 827   for(i=0;i<(ndets-1)/6+1;i++) {
 828     for(k=0;k< (6<ndets-i*6?6:ndets-i*6);k++) {
 829       ind=i*6+k;
 830       det=dynamic_cast<AliModule*>(dets[ind]);
 831       if(det)
 832         printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
 833                det->HiMedium());
 834       else
 835         printf(" %6s: %3d -> %3d;","NULL",0,0);
 836     }
 837     printf("\n");
 838   }
 839             );
 840 }
 841
 842 //_______________________________________________________________________
 843 void AliMC::ReadTransPar()
 844 {
 845   //
 846   // Read filename to set the transport parameters
 847   //
 848
 849
 850   const Int_t kncuts=10;
 851   const Int_t knflags=12;
 852   const Int_t knpars=kncuts+knflags;
 853   const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
 854                                "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
 855                                "BREM","COMP","DCAY","DRAY","HADR","LOSS",
 856                                "MULS","PAIR","PHOT","RAYL","STRA"};
 857   char line[256];
 858   char detName[7];
 859   char* filtmp;
 860   Float_t cut[kncuts];
 861   Int_t flag[knflags];
 862   Int_t i, itmed, iret, ktmed, kz;
 863   FILE *lun;
 864   //
 865   // See whether the file is there
 866   filtmp=gSystem->ExpandPathName(fTransParName.Data());
 867   lun=fopen(filtmp,"r");
 868   delete [] filtmp;
 869   if(!lun) {
 870     AliWarning(Form("File %s does not exist!",fTransParName.Data()));
 871     return;
 872   }
 873   //
 874   while(1) {
 875     // Initialise cuts and flags
 876     for(i=0;i<kncuts;i++) cut[i]=-99;
 877     for(i=0;i<knflags;i++) flag[i]=-99;
 878     itmed=0;
 879     memset(line,0,256);
 880     // Read up to the end of line excluded
 881     iret=fscanf(lun,"%255[^\n]",line);
 882     if(iret<0) {
 883       //End of file
 884       fclose(lun);
 885       return;
 886     }
 887     // Read the end of line
 888     fscanf(lun,"%*c");
 889     if(!iret) continue;
 890     if(line[0]=='*') continue;
 891     // Read the numbers
 892     iret=sscanf(line,"%6s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d %d",
 893                 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
 894                 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
 895                 &flag[8],&flag[9],&flag[10],&flag[11]);
 896     if(!iret) continue;
 897     if(iret<0) {
 898       //reading error
 899       AliWarning(Form("Error reading file %s",fTransParName.Data()));
 900       continue;
 901     }
 902     // Check that the module exist
 903     AliModule *mod = gAlice->GetModule(detName);
 904     if(mod) {
 905       // Get the array of media numbers
 906       TArrayI &idtmed = *mod->GetIdtmed();
 907       // Check that the tracking medium code is valid
 908       if(0<=itmed && itmed < 100) {
 909         ktmed=idtmed[itmed];
 910         if(!ktmed) {
 911           AliWarning(Form("Invalid tracking medium code %d for %s",itmed,mod->GetName()));
 912           continue;
 913         }
 914         // Set energy thresholds
 915         for(kz=0;kz<kncuts;kz++) {
 916           if(cut[kz]>=0) {
 917             AliDebug(2, Form("%-6s set to %10.3E for tracking medium code %4d for %s",
 918                              kpars[kz],cut[kz],itmed,mod->GetName()));
 919             gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
 920           }
 921         }
 922         // Set transport mechanisms
 923         for(kz=0;kz<knflags;kz++) {
 924           if(flag[kz]>=0) {
 925             AliDebug(2, Form("%-6s set to %10d for tracking medium code %4d for %s",
 926                              kpars[kncuts+kz],flag[kz],itmed,mod->GetName()));
 927             gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
 928           }
 929         }
 930       } else {
 931         AliWarning(Form("Invalid medium code %d",itmed));
 932         continue;
 933       }
 934     } else {
 935       AliDebug(1, Form("%s not present",detName));
 936       continue;
 937     }
 938   }
 939 }
 940
 941 //_______________________________________________________________________
 942 void AliMC::SetTransPar(const char *filename)
 943 {
 944   //
 945   // Sets the file name for transport parameters
 946   //
 947   fTransParName = filename;
 948 }
 949
 950 //_______________________________________________________________________
 951 void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
 952 {
 953   //
 954   //  Add a hit to detector id
 955   //
 956   TObjArray &dets = *gAlice->Modules();
 957   if(dets[id]) static_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
 958 }
 959
 960 //_______________________________________________________________________
 961 void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
 962 {
 963   //
 964   // Add digit to detector id
 965   //
 966   TObjArray &dets = *gAlice->Modules();
 967   if(dets[id]) static_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
 968 }
 969
 970 //_______________________________________________________________________
 971 Int_t AliMC::GetCurrentTrackNumber() const {
 972   //
 973   // Returns current track
 974   //
 975   return AliRunLoader::Instance()->Stack()->GetCurrentTrackNumber();
 976 }
 977
 978 //_______________________________________________________________________
 979 void AliMC::DumpPart (Int_t i) const
 980 {
 981   //
 982   // Dumps particle i in the stack
 983   //
 984   AliRunLoader * runloader = AliRunLoader::Instance();
 985    if (runloader->Stack())
 986     runloader->Stack()->DumpPart(i);
 987 }
 988
 989 //_______________________________________________________________________
 990 void AliMC::DumpPStack () const
 991 {
 992   //
 993   // Dumps the particle stack
 994   //
 995   AliRunLoader * runloader = AliRunLoader::Instance();
 996    if (runloader->Stack())
 997     runloader->Stack()->DumpPStack();
 998 }
 999
1000 //_______________________________________________________________________
1001 Int_t AliMC::GetNtrack() const {
1002   //
1003   // Returns number of tracks in stack
1004   //
1005   Int_t ntracks = -1;
1006   AliRunLoader * runloader = AliRunLoader::Instance();
1007    if (runloader->Stack())
1008      ntracks = runloader->Stack()->GetNtrack();
1009    return ntracks;
1010 }
1011
1012 //_______________________________________________________________________
1013 Int_t AliMC::GetPrimary(Int_t track) const
1014 {
1015   //
1016   // return number of primary that has generated track
1017   //
1018   Int_t nprimary = -999;
1019   AliRunLoader * runloader = AliRunLoader::Instance();
1020   if (runloader->Stack())
1021     nprimary = runloader->Stack()->GetPrimary(track);
1022   return nprimary;
1023 }
1024
1025 //_______________________________________________________________________
1026 TParticle* AliMC::Particle(Int_t i) const
1027 {
1028   // Returns the i-th particle from the stack taking into account
1029   // the remaping done by PurifyKine
1030   AliRunLoader * runloader = AliRunLoader::Instance();
1031   if (runloader)
1032    if (runloader->Stack())
1033     return runloader->Stack()->Particle(i);
1034   return 0x0;
1035 }
1036
1037 //_______________________________________________________________________
1038 const TObjArray* AliMC::Particles() const {
1039   //
1040   // Returns pointer to Particles array
1041   //
1042   AliRunLoader * runloader = AliRunLoader::Instance();
1043   if (runloader)
1044    if (runloader->Stack())
1045     return runloader->Stack()->Particles();
1046   return 0x0;
1047 }
1048
1049 //_______________________________________________________________________
1050 void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom,
1051                       const Float_t *vpos, const Float_t *polar, Float_t tof,
1052                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
1053 {
1054 // Delegate to stack
1055 //
1056   AliRunLoader * runloader = AliRunLoader::Instance();
1057   if (runloader)
1058     if (runloader->Stack())
1059       runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
1060                                     mech, ntr, weight, is);
1061 }
1062
1063 //_______________________________________________________________________
1064 void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg,
1065                       Double_t px, Double_t py, Double_t pz, Double_t e,
1066                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
1067                       Double_t polx, Double_t poly, Double_t polz,
1068                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
1069 {
1070   // Delegate to stack
1071   //
1072   AliRunLoader * runloader = AliRunLoader::Instance();
1073   if (runloader)
1074     if (runloader->Stack())
1075       runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
1076                                     polx, poly, polz, mech, ntr, weight, is);
1077 }
1078
1079 //_______________________________________________________________________
1080 void AliMC::SetHighWaterMark(Int_t nt) const
1081 {
1082     //
1083     // Set high water mark for last track in event
1084   AliRunLoader * runloader = AliRunLoader::Instance();
1085   if (runloader)
1086     if (runloader->Stack())
1087       runloader->Stack()->SetHighWaterMark(nt);
1088 }
1089
1090 //_______________________________________________________________________
1091 void AliMC::KeepTrack(Int_t track) const
1092 {
1093   //
1094   // Delegate to stack
1095   //
1096   AliRunLoader * runloader = AliRunLoader::Instance();
1097   if (runloader)
1098     if (runloader->Stack())
1099       runloader->Stack()->KeepTrack(track);
1100 }
1101
1102 //_______________________________________________________________________
1103 void AliMC::FlagTrack(Int_t track) const
1104 {
1105   // Delegate to stack
1106   //
1107   AliRunLoader * runloader = AliRunLoader::Instance();
1108   if (runloader)
1109     if (runloader->Stack())
1110       runloader->Stack()->FlagTrack(track);
1111 }
1112
1113 //_______________________________________________________________________
1114 void AliMC::SetCurrentTrack(Int_t track) const
1115 {
1116   //
1117   // Set current track number
1118   //
1119   AliRunLoader * runloader = AliRunLoader::Instance();
1120   if (runloader)
1121     if (runloader->Stack())
1122       runloader->Stack()->SetCurrentTrack(track);
1123 }
1124
1125 //_______________________________________________________________________
1126 AliTrackReference*  AliMC::AddTrackReference(Int_t label, Int_t id)
1127 {
1128   //
1129   // add a trackrefernce to the list
1130   Int_t primary = GetPrimary(label);
1131   Particle(primary)->SetBit(kKeepBit);
1132
1133   Int_t nref = fTmpTrackReferences.GetEntriesFast();
1134   return new(fTmpTrackReferences[nref]) AliTrackReference(label, id);
1135 }
1136
1137
1138
1139 //_______________________________________________________________________
1140 void AliMC::ResetTrackReferences()
1141 {
1142   //
1143   //  Reset all  references
1144   //
1145     fTmpTrackReferences.Clear();
1146 }
1147
1148 //_______________________________________________________________________
1149 void AliMC::RemapTrackReferencesIDs(const Int_t *map)
1150 {
1151   //
1152   // Remapping track reference
1153   // Called at finish primary
1154   //
1155
1156   Int_t nEntries = fTmpTrackReferences.GetEntries();
1157   for (Int_t i=0; i < nEntries; i++){
1158       AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTmpTrackReferences.UncheckedAt(i));
1159       if (ref) {
1160           Int_t newID = map[ref->GetTrack()];
1161           if (newID>=0) ref->SetTrack(newID);
1162           else {
1163               ref->SetBit(kNotDeleted,kFALSE);
1164               fTmpTrackReferences.RemoveAt(i);
1165           }
1166       } // if ref
1167   }
1168   fTmpTrackReferences.Compress();
1169 }
1170
1171 //_______________________________________________________________________
1172 void AliMC::FixParticleDecaytime()
1173 {
1174     //
1175     // Fix the particle decay time according to rmin and rmax for decays
1176     //
1177
1178     TLorentzVector p;
1179     gMC->TrackMomentum(p);
1180     Double_t tmin, tmax;
1181     Double_t b;
1182
1183     // Transverse velocity
1184     Double_t vt    = p.Pt() / p.E();
1185
1186     if ((b = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField()) > 0.) {     // [kG]
1187
1188         // Radius of helix
1189
1190         Double_t rho   = p.Pt() / 0.0003 / b; // [cm]
1191
1192         // Revolution frequency
1193
1194         Double_t omega = vt / rho;
1195
1196         // Maximum and minimum decay time
1197         //
1198         // Check for curlers first
1199         const Double_t kOvRhoSqr2 = 1./(rho*TMath::Sqrt(2.));
1200         if (fRDecayMax * kOvRhoSqr2 > 1.) return;
1201
1202         //
1203
1204         tmax  = TMath::ACos((1.-fRDecayMax*kOvRhoSqr2)*(1.+fRDecayMax*kOvRhoSqr2)) / omega;   // [ct]
1205         tmin  = TMath::ACos((1.-fRDecayMin*kOvRhoSqr2)*(1.+fRDecayMin*kOvRhoSqr2)) / omega;   // [ct]
1206     } else {
1207         tmax =  fRDecayMax / vt;                                                      // [ct]
1208         tmin =  fRDecayMin / vt;                                                      // [ct]
1209     }
1210
1211     //
1212     // Dial t using the two limits
1213     Double_t t = tmin + (tmax - tmin) * gRandom->Rndm();                              // [ct]
1214     //
1215     //
1216     // Force decay time in transport code
1217     //
1218     gMC->ForceDecayTime(t / 2.99792458e10);
1219 }
1220
1221 void AliMC::MakeTmpTrackRefsTree()
1222 {
1223     // Make the temporary track reference tree
1224     fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
1225     fTmpTreeTR = new TTree("TreeTR", "Track References");
1226     TClonesArray* pRef = &fTmpTrackReferences;
1227     fTmpTreeTR->Branch("TrackReferences", &pRef, 4000);
1228 }
1229
1230 //_______________________________________________________________________
1231 void AliMC::ReorderAndExpandTreeTR()
1232 {
1233 //
1234 //  Reorder and expand the temporary track reference tree in order to match the kinematics tree
1235 //
1236
1237     AliRunLoader *rl = AliRunLoader::Instance();
1238 //
1239 //  TreeTR
1240     AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
1241     rl->MakeTrackRefsContainer();
1242     TTree * treeTR = rl->TreeTR();
1243         // make branch for central track references
1244         TClonesArray* pRef = &fTrackReferences;
1245         treeTR->Branch("TrackReferences", &pRef);
1246
1247     AliStack* stack  = rl->Stack();
1248     Int_t np = stack->GetNprimary();
1249     Int_t nt = fTmpTreeTR->GetEntries();
1250     //
1251     // Loop over tracks and find the secondaries with the help of the kine tree
1252     Int_t ifills = 0;
1253     Int_t it = 0;
1254     for (Int_t ip = np - 1; ip > -1; ip--) {
1255         TParticle *part = stack->Particle(ip);
1256         //printf("Particle %5d %5d %5d %5d %5d \n", ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), part->GetLastDaughter());
1257
1258         // Skip primaries that have not been transported
1259         Int_t dau1  = part->GetFirstDaughter();
1260         Int_t dau2  = -1;
1261         if (!part->TestBit(kTransportBit)) continue;
1262         //
1263         fTmpTreeTR->GetEntry(it++);
1264         Int_t nh = fTmpTrackReferences.GetEntries();
1265         // Determine range of secondaries produced by this primary
1266         if (dau1 > -1) {
1267             Int_t inext = ip - 1;
1268             while (dau2 < 0) {
1269                 if (inext >= 0) {
1270                     part = stack->Particle(inext);
1271                     dau2 =  part->GetFirstDaughter();
1272                     if (!(part->TestBit(kTransportBit)) || dau2 == -1 || dau2 < np) {
1273 //                  if (dau2 == -1 || dau2 < np) {
1274                         dau2 = -1;
1275                     } else {
1276                         dau2--;
1277                     }
1278                 } else {
1279                     dau2 = stack->GetNtrack() - 1;
1280                 }
1281                 inext--;
1282             } // find upper bound
1283         }  // dau2 < 0
1284 //      printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
1285         //
1286         // Loop over reference hits and find secondary label
1287         for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
1288             for (Int_t ih = 0; ih < nh; ih++) {
1289                 AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences.At(ih);
1290                 Int_t label = tr->Label();
1291                 // Skip primaries
1292                 if (label == ip) continue;
1293                 if (label > dau2 || label < dau1)
1294                     AliWarning(Form("Track Reference Label out of range !: %5d %5d %5d \n", label, dau1, dau2));
1295                 if (label == id) {
1296                     // secondary found
1297                     Int_t nref =  fTrackReferences.GetEntriesFast();
1298                     new(fTrackReferences[nref]) AliTrackReference(*tr);
1299                 }
1300             } // hits
1301             treeTR->Fill();
1302             fTrackReferences.Clear();
1303             ifills++;
1304         } // daughters
1305     } // tracks
1306     //
1307     // Now loop again and write the primaries
1308     it = nt - 1;
1309     for (Int_t ip = 0; ip < np; ip++) {
1310         TParticle* part = stack->Particle(ip);
1311 //      if ((part->GetFirstDaughter() == -1 && part->GetStatusCode() <= 1) || part->GetFirstDaughter() >= np)
1312         if (part->TestBit(kTransportBit))
1313         {
1314             // Skip particles that have not been transported
1315             fTmpTreeTR->GetEntry(it--);
1316             Int_t nh = fTmpTrackReferences.GetEntries();
1317             //
1318             // Loop over reference hits and find primary labels
1319             for (Int_t ih = 0; ih < nh; ih++) {
1320                 AliTrackReference* tr = (AliTrackReference*)  fTmpTrackReferences.At(ih);
1321                 Int_t label = tr->Label();
1322                 if (label == ip) {
1323                     Int_t nref = fTrackReferences.GetEntriesFast();
1324                     new(fTrackReferences[nref]) AliTrackReference(*tr);
1325                 }
1326             }
1327         }
1328         treeTR->Fill();
1329         fTrackReferences.Clear();
1330         ifills++;
1331     } // tracks
1332     // Check
1333     if (ifills != stack->GetNtrack())
1334         AliWarning(Form("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, stack->GetNtrack()));
1335 //
1336 //  Clean-up
1337     delete fTmpTreeTR;
1338     fTmpFileTR->Close();
1339     delete fTmpFileTR;
1340     fTmpTrackReferences.Clear();
1341     gSystem->Exec("rm -rf TrackRefsTmp.root");
1342 }
1343