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