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