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