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