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 <TBrowser.h>
  25 #include <TStopwatch.h>
  26 #include <TSystem.h>
  27 #include <TVirtualMC.h>
  28
  29 #include "AliDetector.h"
  30 #include "AliGenerator.h"
  31 #include "AliHeader.h"
  32 #include "AliLego.h"
  33 #include "AliMC.h"
  34 #include "AliMCQA.h"
  35 #include "AliRun.h"
  36 #include "AliStack.h"
  37 #include "AliTrackReference.h"
  38
  39
  40 ClassImp(AliMC)
  41
  42 //_______________________________________________________________________
  43 AliMC::AliMC() :
  44   fGenerator(0),
  45   fEventEnergy(0),
  46   fSummEnergy(0),
  47   fSum2Energy(0),
  48   fTrRmax(1.e10),
  49   fTrZmax(1.e10),
  50   fDebug(0),
  51   fImedia(0),
  52   fTransParName("\0"),
  53   fMCQA(0),
  54   fHitLists(0),
  55   fTrackReferences(0)
  56
  57 {
  58   //default constructor
  59 }
  60
  61 //_______________________________________________________________________
  62 AliMC::AliMC(const char *name, const char *title) :
  63   TVirtualMCApplication(name, title),
  64   fGenerator(0),
  65   fEventEnergy(0),
  66   fSummEnergy(0),
  67   fSum2Energy(0),
  68   fTrRmax(1.e10),
  69   fTrZmax(1.e10),
  70   fDebug(0),
  71   fImedia(new TArrayI(1000)),
  72   fTransParName("\0"),
  73   fMCQA(0),
  74   fHitLists(new TList()),
  75   fTrackReferences(new TClonesArray("AliTrackReference", 100))
  76 {
  77   //constructor
  78   // Set transport parameters
  79   SetTransPar();
  80
  81   // Prepare the tracking medium lists
  82   for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
  83
  84 }
  85
  86 //_______________________________________________________________________
  87 AliMC::AliMC(const AliMC &mc) :
  88   TVirtualMCApplication(mc),
  89   fGenerator(0),
  90   fEventEnergy(0),
  91   fSummEnergy(0),
  92   fSum2Energy(0),
  93   fTrRmax(1.e10),
  94   fTrZmax(1.e10),
  95   fDebug(0),
  96   fImedia(0),
  97   fTransParName("\0"),
  98   fMCQA(0),
  99   fHitLists(0),
 100   fTrackReferences(0)
 101 {
 102   //
 103   // Copy constructor for AliRun
 104   //
 105   mc.Copy(*this);
 106 }
 107
 108 //_______________________________________________________________________
 109 AliMC::~AliMC()
 110 {
 111   //destructor
 112   delete fGenerator;
 113   delete fImedia;
 114   delete fMCQA;
 115   delete fHitLists;
 116   // Delete track references
 117   if (fTrackReferences) {
 118     fTrackReferences->Delete();
 119     delete fTrackReferences;
 120     fTrackReferences     = 0;
 121   }
 122
 123 }
 124
 125 //_______________________________________________________________________
 126 void AliMC::Copy(AliMC &) const
 127 {
 128   //dummy Copy function
 129   Fatal("Copy","Not implemented!\n");
 130 }
 131
 132 //_______________________________________________________________________
 133 void  AliMC::ConstructGeometry()
 134 {
 135   //
 136   // Create modules, materials, geometry
 137   //
 138
 139     TStopwatch stw;
 140     TIter next(gAlice->Modules());
 141     AliModule *detector;
 142     if (GetDebug()) Info("ConstructGeometry","Geometry creation:");
 143     while((detector = dynamic_cast<AliModule*>(next()))) {
 144       stw.Start();
 145       // Initialise detector materials and geometry
 146       detector->CreateMaterials();
 147       detector->CreateGeometry();
 148       printf("%10s R:%.2fs C:%.2fs\n",
 149              detector->GetName(),stw.RealTime(),stw.CpuTime());
 150     }
 151 }
 152
 153 //_______________________________________________________________________
 154 void  AliMC::InitGeometry()
 155 {
 156   //
 157   // Initialize detectors and display geometry
 158   //
 159
 160    printf("Initialisation:\n");
 161     TStopwatch stw;
 162     TIter next(gAlice->Modules());
 163     AliModule *detector;
 164     while((detector = dynamic_cast<AliModule*>(next()))) {
 165       stw.Start();
 166       // Initialise detector and display geometry
 167       detector->Init();
 168       detector->BuildGeometry();
 169       printf("%10s R:%.2fs C:%.2fs\n",
 170              detector->GetName(),stw.RealTime(),stw.CpuTime());
 171     }
 172
 173 }
 174
 175 //_______________________________________________________________________
 176 void  AliMC::GeneratePrimaries()
 177 {
 178   //
 179   // Generate primary particles and fill them in the stack.
 180   //
 181
 182   Generator()->Generate();
 183 }
 184
 185 //_______________________________________________________________________
 186 void AliMC::SetGenerator(AliGenerator *generator)
 187 {
 188   //
 189   // Load the event generator
 190   //
 191   if(!fGenerator) fGenerator = generator;
 192 }
 193
 194 //_______________________________________________________________________
 195 void AliMC::ResetGenerator(AliGenerator *generator)
 196 {
 197   //
 198   // Load the event generator
 199   //
 200   if(fGenerator)
 201     if(generator)
 202       Warning("ResetGenerator","Replacing generator %s with %s\n",
 203               fGenerator->GetName(),generator->GetName());
 204     else
 205       Warning("ResetGenerator","Replacing generator %s with NULL\n",
 206               fGenerator->GetName());
 207   fGenerator = generator;
 208 }
 209
 210 //_______________________________________________________________________
 211 void AliMC::FinishRun()
 212 {
 213   // Clean generator information
 214   if (GetDebug()) Info("FinishRun"," fGenerator->FinishRun()");
 215   fGenerator->FinishRun();
 216
 217   //Output energy summary tables
 218   if (GetDebug()) Info("FinishRun"," EnergySummary()");
 219   EnergySummary();
 220
 221 }
 222
 223 //_______________________________________________________________________
 224 void AliMC::BeginPrimary()
 225 {
 226   //
 227   // Called  at the beginning of each primary track
 228   //
 229
 230   // Reset Hits info
 231   ResetHits();
 232   ResetTrackReferences();
 233
 234 }
 235
 236 //_______________________________________________________________________
 237 void AliMC::PreTrack()
 238 {
 239   // Actions before the track's transport
 240      TObjArray &dets = *gAlice->Modules();
 241      AliModule *module;
 242
 243      for(Int_t i=0; i<=gAlice->GetNdets(); i++)
 244        if((module = dynamic_cast<AliModule*>(dets[i])))
 245          module->PreTrack();
 246
 247      fMCQA->PreTrack();
 248 }
 249
 250 //_______________________________________________________________________
 251 void AliMC::Stepping()
 252 {
 253   //
 254   // Called at every step during transport
 255   //
 256
 257   Int_t id = DetFromMate(gMC->GetMedium());
 258   if (id < 0) return;
 259
 260   //
 261   // --- If lego option, do it and leave
 262   if (gAlice->Lego())
 263     gAlice->Lego()->StepManager();
 264   else {
 265     Int_t copy;
 266     //Update energy deposition tables
 267     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
 268
 269     //Call the appropriate stepping routine;
 270     AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
 271     if(det && det->StepManagerIsEnabled()) {
 272       fMCQA->StepManager(id);
 273       det->StepManager();
 274     }
 275   }
 276 }
 277
 278 //_______________________________________________________________________
 279 void AliMC::EnergySummary()
 280 {
 281   //
 282   // Print summary of deposited energy
 283   //
 284
 285   Int_t ndep=0;
 286   Float_t edtot=0;
 287   Float_t ed, ed2;
 288   Int_t kn, i, left, j, id;
 289   const Float_t kzero=0;
 290   Int_t ievent=gAlice->GetRunLoader()->GetHeader()->GetEvent()+1;
 291   //
 292   // Energy loss information
 293   if(ievent) {
 294     printf("***************** Energy Loss Information per event (GEV) *****************\n");
 295     for(kn=1;kn<fEventEnergy.GetSize();kn++) {
 296       ed=fSummEnergy[kn];
 297       if(ed>0) {
 298         fEventEnergy[ndep]=kn;
 299         if(ievent>1) {
 300           ed=ed/ievent;
 301           ed2=fSum2Energy[kn];
 302           ed2=ed2/ievent;
 303           ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
 304         } else
 305           ed2=99;
 306         fSummEnergy[ndep]=ed;
 307         fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
 308         edtot+=ed;
 309         ndep++;
 310       }
 311     }
 312     for(kn=0;kn<(ndep-1)/3+1;kn++) {
 313       left=ndep-kn*3;
 314       for(i=0;i<(3<left?3:left);i++) {
 315         j=kn*3+i;
 316         id=Int_t (fEventEnergy[j]+0.1);
 317         printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
 318       }
 319       printf("\n");
 320     }
 321     //
 322     // Relative energy loss in different detectors
 323     printf("******************** Relative Energy Loss per event ********************\n");
 324     printf("Total energy loss per event %10.3f GeV\n",edtot);
 325     for(kn=0;kn<(ndep-1)/5+1;kn++) {
 326       left=ndep-kn*5;
 327       for(i=0;i<(5<left?5:left);i++) {
 328         j=kn*5+i;
 329         id=Int_t (fEventEnergy[j]+0.1);
 330         printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
 331       }
 332       printf("\n");
 333     }
 334     for(kn=0;kn<75;kn++) printf("*");
 335     printf("\n");
 336   }
 337   //
 338   // Reset the TArray's
 339   //  fEventEnergy.Set(0);
 340   //  fSummEnergy.Set(0);
 341   //  fSum2Energy.Set(0);
 342 }
 343
 344 //_____________________________________________________________________________
 345 void AliMC::BeginEvent()
 346 {
 347   //
 348   // Clean-up previous event
 349   // Energy scores
 350   if (GetDebug())
 351    {
 352      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 353      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 354      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 355      Info("BeginEvent","          BEGINNING EVENT               ");
 356      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 357      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 358      Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
 359    }
 360
 361     AliRunLoader *runloader=gAlice->GetRunLoader();
 362
 363   /*******************************/
 364   /*   Clean after eventual      */
 365   /*   previous event            */
 366   /*******************************/
 367
 368
 369   //Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors),
 370   gAlice->SetEventNrInRun(gAlice->GetEventNrInRun()+1);
 371   runloader->SetEventNumber(gAlice->GetEventNrInRun());// sets new files, cleans the previous event stuff, if necessary, etc.,
 372   if (GetDebug()) Info("BeginEvent","EventNr is %d",gAlice->GetEventNrInRun());
 373
 374   fEventEnergy.Reset();
 375     // Clean detector information
 376
 377   if (runloader->Stack())
 378     runloader->Stack()->Reset();//clean stack -> tree is unloaded
 379   else
 380     runloader->MakeStack();//or make a new one
 381
 382   if(gAlice->Lego() == 0x0)
 383    {
 384      if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTree(K)");
 385      runloader->MakeTree("K");
 386    }
 387
 388   if (GetDebug()) Info("BeginEvent","  gMC->SetStack(fRunLoader->Stack())");
 389   gMC->SetStack(gAlice->GetRunLoader()->Stack());//Was in InitMC - but was moved here
 390                                      //because we don't have guarantee that
 391                                      //stack pointer is not going to change from event to event
 392                          //since it bellobgs to header and is obtained via RunLoader
 393   //
 394   //  Reset all Detectors & kinematics & make/reset trees
 395   //
 396
 397   runloader->GetHeader()->Reset(gAlice->GetRunNumber(),gAlice->GetEvNumber(),
 398     gAlice->GetEventNrInRun());
 399 //  fRunLoader->WriteKinematics("OVERWRITE");  is there any reason to rewrite here since MakeTree does so
 400
 401   if(gAlice->Lego())
 402    {
 403     gAlice->Lego()->BeginEvent();
 404     return;
 405    }
 406
 407
 408   if (GetDebug()) Info("BeginEvent","  ResetHits()");
 409   ResetHits();
 410
 411   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTree(H)");
 412   runloader->MakeTree("H");
 413
 414   if (GetDebug()) Info("BeginEvent","  fRunLoader->MakeTrackRefsContainer()");
 415   runloader->MakeTrackRefsContainer();//for insurance
 416
 417
 418   //create new branches and SetAdresses
 419   TIter next(gAlice->Modules());
 420   AliModule *detector;
 421   while((detector = (AliModule*)next()))
 422    {
 423     if (GetDebug()) Info("BeginEvent","  %s->MakeBranch(H)",detector->GetName());
 424     detector->MakeBranch("H");
 425     if (GetDebug()) Info("BeginEvent","  %s->MakeBranchTR()",detector->GetName());
 426     detector->MakeBranchTR();
 427     if (GetDebug()) Info("BeginEvent","  %s->SetTreeAddress()",detector->GetName());
 428     detector->SetTreeAddress();
 429    }
 430 }
 431
 432 //_______________________________________________________________________
 433 void AliMC::ResetHits()
 434 {
 435   //
 436   //  Reset all Detectors hits
 437   //
 438   TIter next(gAlice->Modules());
 439   AliModule *detector;
 440   while((detector = dynamic_cast<AliModule*>(next()))) {
 441      detector->ResetHits();
 442   }
 443 }
 444
 445 //_______________________________________________________________________
 446 void AliMC::PostTrack()
 447 {
 448   // Posts tracks for each module
 449      TObjArray &dets = *gAlice->Modules();
 450      AliModule *module;
 451
 452      for(Int_t i=0; i<=gAlice->GetNdets(); i++)
 453        if((module = dynamic_cast<AliModule*>(dets[i])))
 454          module->PostTrack();
 455 }
 456
 457 //_______________________________________________________________________
 458 void AliMC::FinishPrimary()
 459 {
 460   //
 461   // Called  at the end of each primary track
 462   //
 463   AliRunLoader *runloader=gAlice->GetRunLoader();
 464   //  static Int_t count=0;
 465   //  const Int_t times=10;
 466   // This primary is finished, purify stack
 467   runloader->Stack()->PurifyKine();
 468
 469   TIter next(gAlice->Modules());
 470   AliModule *detector;
 471   while((detector = dynamic_cast<AliModule*>(next()))) {
 472     detector->FinishPrimary();
 473     AliLoader* loader = detector->GetLoader();
 474     if(loader)
 475      {
 476        TTree* treeH = loader->TreeH();
 477        if (treeH) treeH->Fill(); //can be Lego run and treeH can not exist
 478      }
 479   }
 480
 481   // Write out track references if any
 482   if (runloader->TreeTR()) runloader->TreeTR()->Fill();
 483 }
 484
 485 //_______________________________________________________________________
 486 void AliMC::FinishEvent()
 487 {
 488   //
 489   // Called at the end of the event.
 490   //
 491
 492   //
 493   if(gAlice->Lego()) gAlice->Lego()->FinishEvent();
 494
 495   TIter next(gAlice->Modules());
 496   AliModule *detector;
 497   while((detector = dynamic_cast<AliModule*>(next()))) {
 498     detector->FinishEvent();
 499   }
 500
 501   //Update the energy deposit tables
 502   Int_t i;
 503   for(i=0;i<fEventEnergy.GetSize();i++)
 504    {
 505     fSummEnergy[i]+=fEventEnergy[i];
 506     fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
 507    }
 508
 509   AliRunLoader *runloader=gAlice->GetRunLoader();
 510
 511   AliHeader* header = runloader->GetHeader();
 512   AliStack* stack = runloader->Stack();
 513   if ( (header == 0x0) || (stack == 0x0) )
 514    {//check if we got header and stack. If not cry and exit aliroot
 515     Fatal("AliRun","Can not get the stack or header from LOADER");
 516     return;//never reached
 517    }
 518   // Update Header information
 519   header->SetNprimary(stack->GetNprimary());
 520   header->SetNtrack(stack->GetNtrack());
 521
 522
 523   // Write out the kinematics
 524   stack->FinishEvent();
 525
 526   // Write out the event Header information
 527   TTree* treeE = runloader->TreeE();
 528   if (treeE)
 529    {
 530       header->SetStack(stack);
 531       treeE->Fill();
 532    }
 533   else
 534    {
 535     Error("FinishEvent","Can not get TreeE from RL");
 536    }
 537
 538   if(gAlice->Lego() == 0x0)
 539    {
 540      runloader->WriteKinematics("OVERWRITE");
 541      runloader->WriteTrackRefs("OVERWRITE");
 542      runloader->WriteHits("OVERWRITE");
 543    }
 544
 545   if (GetDebug())
 546    {
 547      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 548      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 549      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 550      Info("FinishEvent","          FINISHING EVENT               ");
 551      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 552      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 553      Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
 554    }
 555 }
 556
 557 //_______________________________________________________________________
 558 void AliMC::Field(const Double_t* x, Double_t* b) const
 559 {
 560   // Calculates field "b" at point "x"
 561     gAlice->Field(x,b);
 562 }
 563
 564 //_______________________________________________________________________
 565 void AliMC::Init()
 566 {
 567   // MC initialization
 568
 569    //=================Create Materials and geometry
 570    gMC->Init();
 571    gMC->DefineParticles();  //Create standard MC particles
 572
 573    //Read the cuts for all materials
 574    ReadTransPar();
 575    //Build the special IMEDIA table
 576    MediaTable();
 577
 578    //Compute cross-sections
 579    gMC->BuildPhysics();
 580
 581    //Write Geometry object to current file.
 582    gAlice->GetRunLoader()->WriteGeometry();
 583
 584    //Initialise geometry deposition table
 585    fEventEnergy.Set(gMC->NofVolumes()+1);
 586    fSummEnergy.Set(gMC->NofVolumes()+1);
 587    fSum2Energy.Set(gMC->NofVolumes()+1);
 588
 589    //
 590    fMCQA = new AliMCQA(gAlice->GetNdets());
 591
 592    // Register MC in configuration
 593    AliConfig::Instance()->Add(gMC);
 594
 595 }
 596
 597 //_______________________________________________________________________
 598 void AliMC::MediaTable()
 599 {
 600   //
 601   // Built media table to get from the media number to
 602   // the detector id
 603   //
 604
 605   Int_t kz, nz, idt, lz, i, k, ind;
 606   //  Int_t ibeg;
 607   TObjArray &dets = *gAlice->Detectors();
 608   AliModule *det;
 609   Int_t ndets=gAlice->GetNdets();
 610   //
 611   // For all detectors
 612   for (kz=0;kz<ndets;kz++) {
 613     // If detector is defined
 614     if((det=dynamic_cast<AliModule*>(dets[kz]))) {
 615         TArrayI &idtmed = *(det->GetIdtmed());
 616         for(nz=0;nz<100;nz++) {
 617         // Find max and min material number
 618         if((idt=idtmed[nz])) {
 619           det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
 620           det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
 621         }
 622       }
 623       if(det->LoMedium() > det->HiMedium()) {
 624         det->LoMedium() = 0;
 625         det->HiMedium() = 0;
 626       } else {
 627         if(det->HiMedium() > fImedia->GetSize()) {
 628           Error("MediaTable","Increase fImedia from %d to %d",
 629                 fImedia->GetSize(),det->HiMedium());
 630           return;
 631         }
 632         // Tag all materials in rage as belonging to detector kz
 633         for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
 634           (*fImedia)[lz]=kz;
 635         }
 636       }
 637     }
 638   }
 639   //
 640   // Print summary table
 641   printf(" Traking media ranges:\n");
 642   for(i=0;i<(ndets-1)/6+1;i++) {
 643     for(k=0;k< (6<ndets-i*6?6:ndets-i*6);k++) {
 644       ind=i*6+k;
 645       det=dynamic_cast<AliModule*>(dets[ind]);
 646       if(det)
 647         printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
 648                det->HiMedium());
 649       else
 650         printf(" %6s: %3d -> %3d;","NULL",0,0);
 651     }
 652     printf("\n");
 653   }
 654 }
 655
 656 //_______________________________________________________________________
 657 void AliMC::ReadTransPar()
 658 {
 659   //
 660   // Read filename to set the transport parameters
 661   //
 662
 663
 664   const Int_t kncuts=10;
 665   const Int_t knflags=11;
 666   const Int_t knpars=kncuts+knflags;
 667   const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
 668                                "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
 669                                "BREM","COMP","DCAY","DRAY","HADR","LOSS",
 670                                "MULS","PAIR","PHOT","RAYL"};
 671   char line[256];
 672   char detName[7];
 673   char* filtmp;
 674   Float_t cut[kncuts];
 675   Int_t flag[knflags];
 676   Int_t i, itmed, iret, ktmed, kz;
 677   FILE *lun;
 678   //
 679   // See whether the file is there
 680   filtmp=gSystem->ExpandPathName(fTransParName.Data());
 681   lun=fopen(filtmp,"r");
 682   delete [] filtmp;
 683   if(!lun) {
 684     Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
 685     return;
 686   }
 687   //
 688   if(fDebug) {
 689     printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
 690     printf(" *%59s\n","*");
 691     printf(" *       Please check carefully what you are doing!%10s\n","*");
 692     printf(" *%59s\n","*");
 693   }
 694   //
 695   while(1) {
 696     // Initialise cuts and flags
 697     for(i=0;i<kncuts;i++) cut[i]=-99;
 698     for(i=0;i<knflags;i++) flag[i]=-99;
 699     itmed=0;
 700     for(i=0;i<256;i++) line[i]='\0';
 701     // Read up to the end of line excluded
 702     iret=fscanf(lun,"%[^\n]",line);
 703     if(iret<0) {
 704       //End of file
 705       fclose(lun);
 706       if(fDebug){
 707         printf(" *%59s\n","*");
 708         printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
 709       }
 710       return;
 711     }
 712     // Read the end of line
 713     fscanf(lun,"%*c");
 714     if(!iret) continue;
 715     if(line[0]=='*') continue;
 716     // Read the numbers
 717     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",
 718                 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
 719                 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
 720                 &flag[8],&flag[9],&flag[10]);
 721     if(!iret) continue;
 722     if(iret<0) {
 723       //reading error
 724       Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
 725       continue;
 726     }
 727     // Check that the module exist
 728     AliModule *mod = gAlice->GetModule(detName);
 729     if(mod) {
 730       // Get the array of media numbers
 731       TArrayI &idtmed = *mod->GetIdtmed();
 732       // Check that the tracking medium code is valid
 733       if(0<=itmed && itmed < 100) {
 734         ktmed=idtmed[itmed];
 735         if(!ktmed) {
 736           Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
 737           continue;
 738         }
 739         // Set energy thresholds
 740         for(kz=0;kz<kncuts;kz++) {
 741           if(cut[kz]>=0) {
 742             if(fDebug) printf(" *  %-6s set to %10.3E for tracking medium code %4d for %s\n",
 743                    kpars[kz],cut[kz],itmed,mod->GetName());
 744             gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
 745           }
 746         }
 747         // Set transport mechanisms
 748         for(kz=0;kz<knflags;kz++) {
 749           if(flag[kz]>=0) {
 750             if(fDebug) printf(" *  %-6s set to %10d for tracking medium code %4d for %s\n",
 751                    kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
 752             gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
 753           }
 754         }
 755       } else {
 756         Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
 757         continue;
 758       }
 759     } else {
 760       if(fDebug) printf("%s::ReadTransParModule: %s not present\n",ClassName(),detName);
 761       continue;
 762     }
 763   }
 764 }
 765
 766 //_______________________________________________________________________
 767 void AliMC::SetTransPar(const char *filename)
 768 {
 769   //
 770   // Sets the file name for transport parameters
 771   //
 772   fTransParName = filename;
 773 }
 774
 775 //_______________________________________________________________________
 776 void AliMC::Browse(TBrowser *b)
 777 {
 778   //
 779   // Called when the item "Run" is clicked on the left pane
 780   // of the Root browser.
 781   // It displays the Root Trees and all detectors.
 782   //
 783   //detectors are in folders anyway
 784   b->Add(fMCQA,"AliMCQA");
 785 }
 786
 787 //PH
 788 //_______________________________________________________________________
 789 void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
 790 {
 791   //
 792   //  Add a hit to detector id
 793   //
 794   TObjArray &dets = *gAlice->Modules();
 795   if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
 796 }
 797
 798 //_______________________________________________________________________
 799 void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
 800 {
 801   //
 802   // Add digit to detector id
 803   //
 804   TObjArray &dets = *gAlice->Modules();
 805   if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
 806 }
 807
 808 //_______________________________________________________________________
 809 Int_t AliMC::GetCurrentTrackNumber() const {
 810   //
 811   // Returns current track
 812   //
 813   return gAlice->GetRunLoader()->Stack()->GetCurrentTrackNumber();
 814 }
 815
 816 //_______________________________________________________________________
 817 void AliMC::DumpPart (Int_t i) const
 818 {
 819   //
 820   // Dumps particle i in the stack
 821   //
 822   AliRunLoader * runloader = gAlice->GetRunLoader();
 823    if (runloader->Stack())
 824     runloader->Stack()->DumpPart(i);
 825 }
 826
 827 //_______________________________________________________________________
 828 void AliMC::DumpPStack () const
 829 {
 830   //
 831   // Dumps the particle stack
 832   //
 833   AliRunLoader * runloader = gAlice->GetRunLoader();
 834    if (runloader->Stack())
 835     runloader->Stack()->DumpPStack();
 836 }
 837
 838 //_______________________________________________________________________
 839 Int_t AliMC::GetNtrack() const {
 840   //
 841   // Returns number of tracks in stack
 842   //
 843   Int_t ntracks = -1;
 844   AliRunLoader * runloader = gAlice->GetRunLoader();
 845    if (runloader->Stack())
 846      ntracks = runloader->Stack()->GetNtrack();
 847    return ntracks;
 848 }
 849
 850 //_______________________________________________________________________
 851 Int_t AliMC::GetPrimary(Int_t track) const
 852 {
 853   //
 854   // return number of primary that has generated track
 855   //
 856   Int_t nprimary = -999;
 857   AliRunLoader * runloader = gAlice->GetRunLoader();
 858   if (runloader->Stack())
 859     nprimary = runloader->Stack()->GetPrimary(track);
 860   return nprimary;
 861 }
 862
 863 //_______________________________________________________________________
 864 TParticle* AliMC::Particle(Int_t i) const
 865 {
 866   // Returns the i-th particle from the stack taking into account
 867   // the remaping done by PurifyKine
 868   AliRunLoader * runloader = gAlice->GetRunLoader();
 869   if (runloader)
 870    if (runloader->Stack())
 871     return runloader->Stack()->Particle(i);
 872   return 0x0;
 873 }
 874
 875 //_______________________________________________________________________
 876 TObjArray* AliMC::Particles() const {
 877   //
 878   // Returns pointer to Particles array
 879   //
 880   AliRunLoader * runloader = gAlice->GetRunLoader();
 881   if (runloader)
 882    if (runloader->Stack())
 883     return runloader->Stack()->Particles();
 884   return 0x0;
 885 }
 886
 887 //_______________________________________________________________________
 888 void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
 889                       Float_t *vpos, Float_t *polar, Float_t tof,
 890                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
 891 {
 892 // Delegate to stack
 893 //
 894   AliRunLoader * runloader = gAlice->GetRunLoader();
 895   if (runloader)
 896     if (runloader->Stack())
 897       runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
 898                                     mech, ntr, weight, is);
 899 }
 900
 901 //_______________________________________________________________________
 902 void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg,
 903                       Double_t px, Double_t py, Double_t pz, Double_t e,
 904                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
 905                       Double_t polx, Double_t poly, Double_t polz,
 906                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
 907 {
 908   // Delegate to stack
 909   //
 910   AliRunLoader * runloader = gAlice->GetRunLoader();
 911   if (runloader)
 912     if (runloader->Stack())
 913       runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
 914                                     polx, poly, polz, mech, ntr, weight, is);
 915 }
 916
 917 //_______________________________________________________________________
 918 void AliMC::SetHighWaterMark(const Int_t nt)
 919 {
 920     //
 921     // Set high water mark for last track in event
 922   AliRunLoader * runloader = gAlice->GetRunLoader();
 923   if (runloader)
 924     if (runloader->Stack())
 925       runloader->Stack()->SetHighWaterMark(nt);
 926 }
 927
 928 //_______________________________________________________________________
 929 void AliMC::KeepTrack(const Int_t track)
 930 {
 931   //
 932   // Delegate to stack
 933   //
 934   AliRunLoader * runloader = gAlice->GetRunLoader();
 935   if (runloader)
 936     if (runloader->Stack())
 937       runloader->Stack()->KeepTrack(track);
 938 }
 939
 940 //_______________________________________________________________________
 941 void AliMC::FlagTrack(Int_t track)
 942 {
 943   // Delegate to stack
 944   //
 945   AliRunLoader * runloader = gAlice->GetRunLoader();
 946   if (runloader)
 947     if (runloader->Stack())
 948       runloader->Stack()->FlagTrack(track);
 949 }
 950
 951 //_______________________________________________________________________
 952 void AliMC::SetCurrentTrack(Int_t track)
 953 {
 954   //
 955   // Set current track number
 956   //
 957   AliRunLoader * runloader = gAlice->GetRunLoader();
 958   if (runloader)
 959     if (runloader->Stack())
 960       runloader->Stack()->SetCurrentTrack(track);
 961 }
 962
 963 //_______________________________________________________________________
 964 void  AliMC::AddTrackReference(Int_t label){
 965   //
 966   // add a trackrefernce to the list
 967   if (!fTrackReferences) {
 968     cerr<<"Container trackrefernce not active\n";
 969     return;
 970   }
 971   Int_t nref = fTrackReferences->GetEntriesFast();
 972   TClonesArray &lref = *fTrackReferences;
 973   new(lref[nref]) AliTrackReference(label);
 974 }
 975
 976
 977
 978 //_______________________________________________________________________
 979 void AliMC::ResetTrackReferences()
 980 {
 981   //
 982   //  Reset all  references
 983   //
 984   if (fTrackReferences)   fTrackReferences->Clear();
 985
 986   TIter next(gAlice->Modules());
 987   AliModule *detector;
 988   while((detector = dynamic_cast<AliModule*>(next()))) {
 989      detector->ResetTrackReferences();
 990   }
 991 }
 992
 993 void AliMC::RemapTrackReferencesIDs(Int_t *map)
 994 {
 995   //
 996   // Remapping track reference
 997   // Called at finish primary
 998   //
 999   if (!fTrackReferences) return;
1000   for (Int_t i=0;i<fTrackReferences->GetEntries();i++){
1001     AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTrackReferences->UncheckedAt(i));
1002     if (ref) {
1003       Int_t newID = map[ref->GetTrack()];
1004       if (newID>=0) ref->SetTrack(newID);
1005       else {
1006         //ref->SetTrack(-1);
1007         ref->SetBit(kNotDeleted,kFALSE);
1008         fTrackReferences->RemoveAt(i);
1009       }
1010     }
1011   }
1012   fTrackReferences->Compress();
1013 }