MC-dependent part of AliRun extracted in AliMC (F.Carminati)
[u/mrichter/AliRoot.git] / 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 }