Corrected protection.
[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 <TArrayI.h>
26 #include <TClonesArray.h>
27 #include <TFile.h>
28 #include <TGeoGlobalMagField.h>
29 #include <TGeoManager.h>
30 #include <TParticle.h>
31 #include <TROOT.h>
32 #include <TStopwatch.h>
33 #include <TSystem.h>
34 #include <TVirtualMC.h>
35  
36 #include "AliCDBEntry.h"
37 #include "AliCDBManager.h"
38 #include "AliCDBStorage.h"
39 #include "AliDetector.h"
40 #include "AliGenerator.h"
41 #include "AliGeomManager.h"
42 #include "AliHeader.h"
43 #include "AliHit.h"
44 #include "AliLego.h"
45 #include "AliLog.h"
46 #include "AliMC.h"
47 #include "AliMagF.h"
48 #include "AliRun.h"
49 #include "AliSimulation.h"
50 #include "AliStack.h"
51 #include "AliTrackReference.h"
52
53 ClassImp(AliMC)
54
55 //_______________________________________________________________________
56 AliMC::AliMC() :
57   fGenerator(0),
58   fEventEnergy(0),
59   fSummEnergy(0),
60   fSum2Energy(0),
61   fTrRmax(1.e10),
62   fTrZmax(1.e10),
63   fRDecayMax(1.e10),
64   fRDecayMin(-1.),
65   fDecayPdg(0),
66   fImedia(0),
67   fTransParName("\0"),
68   fHitLists(0),
69   fTmpTreeTR(0),
70   fTmpFileTR(0),
71   fTrackReferences(),
72   fTmpTrackReferences()
73
74 {
75   //default constructor
76   DecayLimits();
77 }
78
79 //_______________________________________________________________________
80 AliMC::AliMC(const char *name, const char *title) :
81   TVirtualMCApplication(name, title),
82   fGenerator(0),
83   fEventEnergy(0),
84   fSummEnergy(0),
85   fSum2Energy(0),
86   fTrRmax(1.e10),
87   fTrZmax(1.e10),
88   fRDecayMax(1.e10),
89   fRDecayMin(-1.),
90   fDecayPdg(0),
91   fImedia(new TArrayI(1000)),
92   fTransParName("\0"),
93   fHitLists(new TList()),
94   fTmpTreeTR(0),
95   fTmpFileTR(0),
96   fTrackReferences("AliTrackReference", 100),
97   fTmpTrackReferences("AliTrackReference", 100)
98 {
99   //constructor
100   // Set transport parameters
101   SetTransPar();
102   DecayLimits();
103   // Prepare the tracking medium lists
104   for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99;
105 }
106
107 //_______________________________________________________________________
108 AliMC::~AliMC()
109 {
110   //destructor
111   delete fGenerator;
112   delete fImedia;
113   delete fHitLists;
114   // Delete track references
115 }
116
117 //_______________________________________________________________________
118 void  AliMC::ConstructGeometry() 
119 {
120   //
121   // Either load geometry from file or create it through usual
122   // loop on detectors. In the first case the method
123   // AliModule::CreateMaterials() only builds fIdtmed and is postponed
124   // at InitGeometry().
125   //
126
127   if(AliSimulation::Instance()->IsGeometryFromFile()){ //load geometry either from CDB or from file
128     if(IsGeometryFromCDB()){
129       AliInfo("Loading geometry from CDB default storage");
130       AliCDBPath path("GRP","Geometry","Data");
131       AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath());
132       if(!entry) AliFatal("Unable to load geometry from CDB!");
133       entry->SetOwner(0);
134       gGeoManager = (TGeoManager*) entry->GetObject();
135       if (!gGeoManager) AliFatal("TGeoManager object not found in the specified CDB entry!");
136     }else{
137       // Load geometry
138       const char *geomfilename = AliSimulation::Instance()->GetGeometryFile();
139       if(gSystem->ExpandPathName(geomfilename)){
140         AliInfo(Form("Loading geometry from file:\n %40s",geomfilename));
141         TGeoManager::Import(geomfilename);
142       }else{
143         AliInfo(Form("Geometry file %40s not found!\n",geomfilename));
144         return;
145       }
146     }
147     gMC->SetRootGeometry();
148   }else{
149     // Create modules, materials, geometry
150     if (!gGeoManager) new TGeoManager("ALICE", "ALICE geometry");
151     TStopwatch stw;
152     TIter next(gAlice->Modules());
153     AliModule *detector;
154     AliDebug(1, "Geometry creation:");
155     while((detector = dynamic_cast<AliModule*>(next()))) {
156       stw.Start();
157       // Initialise detector materials and geometry
158       detector->CreateMaterials();
159       detector->CreateGeometry();
160       AliInfo(Form("%10s R:%.2fs C:%.2fs",
161                    detector->GetName(),stw.RealTime(),stw.CpuTime()));
162     }
163   }
164   
165 }
166
167 //_______________________________________________________________________
168 Bool_t  AliMC::MisalignGeometry() 
169 {
170   // Call misalignment code if AliSimulation object was defined.
171   
172   if(!AliSimulation::Instance()->IsGeometryFromFile()){
173     //Set alignable volumes for the whole geometry
174     SetAllAlignableVolumes();
175   }
176   // Misalign geometry via AliSimulation instance
177   if (!AliSimulation::Instance()) return kFALSE;
178   AliGeomManager::SetGeometry(gGeoManager);
179   if(!AliGeomManager::CheckSymNamesLUT("ALL"))
180     AliFatal("Current loaded geometry differs in the definition of symbolic names!");
181   
182   return AliSimulation::Instance()->MisalignGeometry(AliRunLoader::Instance());
183 }   
184
185 //_______________________________________________________________________
186 void  AliMC::ConstructOpGeometry() 
187 {
188   //
189   // Loop all detector modules and call DefineOpticalProperties() method 
190   //
191
192   TIter next(gAlice->Modules());
193   AliModule *detector;
194   AliInfo("Optical properties definition");
195   while((detector = dynamic_cast<AliModule*>(next()))) {
196     // Initialise detector geometry
197     if(AliSimulation::Instance()->IsGeometryFromFile()) detector->CreateMaterials();
198     // Initialise detector optical properties
199     detector->DefineOpticalProperties();
200   }  
201 }
202
203 //_______________________________________________________________________
204 void  AliMC::InitGeometry()
205
206   //
207   // Initialize detectors
208   //
209
210   AliInfo("Initialisation:");
211   TStopwatch stw;
212   TIter next(gAlice->Modules());
213   AliModule *detector;
214   while((detector = dynamic_cast<AliModule*>(next()))) {
215     stw.Start();
216     detector->Init();
217     AliInfo(Form("%10s R:%.2fs C:%.2fs",
218                  detector->GetName(),stw.RealTime(),stw.CpuTime()));
219   }
220 }
221
222 //_______________________________________________________________________
223 void AliMC::SetGeometryFromCDB()
224 {
225   // Set the loading of geometry from cdb instead of creating it
226   // A default CDB storage needs to be set before this method is called
227   if(AliCDBManager::Instance()->IsDefaultStorageSet() &&
228      AliCDBManager::Instance()->GetRun() >= 0)
229     AliSimulation::Instance()->SetGeometryFile("*OCDB*");
230   else
231     AliError("Loading of geometry from CDB ignored. First set a default CDB storage!");
232 }
233
234 //_______________________________________________________________________
235 Bool_t AliMC::IsGeometryFromCDB() const
236 {
237   return (strcmp(AliSimulation::Instance()->GetGeometryFile(),"*OCDB*")==0);
238 }
239
240 //_______________________________________________________________________
241 void  AliMC::SetAllAlignableVolumes()
242
243   //
244   // Add alignable volumes (TGeoPNEntries) looping on all
245   // active modules
246   //
247
248   AliInfo(Form("Setting entries for all alignable volumes of active detectors"));
249   AliModule *detector;
250   TIter next(gAlice->Modules());
251   while((detector = dynamic_cast<AliModule*>(next()))) {
252     detector->AddAlignableVolumes();
253   }
254 }
255
256 //_______________________________________________________________________
257 void  AliMC::GeneratePrimaries() 
258
259   //
260   // Generate primary particles and fill them in the stack.
261   //
262
263   Generator()->Generate();
264 }
265
266 //_______________________________________________________________________
267 void AliMC::SetGenerator(AliGenerator *generator)
268 {
269   //
270   // Load the event generator
271   //
272   if(!fGenerator) fGenerator = generator;
273 }
274
275 //_______________________________________________________________________
276 void AliMC::ResetGenerator(AliGenerator *generator)
277 {
278   //
279   // Load the event generator
280   //
281   if(fGenerator) {
282     if(generator) {
283       AliWarning(Form("Replacing generator %s with %s",
284                       fGenerator->GetName(),generator->GetName()));
285     }
286     else {
287       AliWarning(Form("Replacing generator %s with NULL",
288                       fGenerator->GetName()));
289     }
290   }
291   fGenerator = generator;
292 }
293
294 //_______________________________________________________________________
295 void AliMC::FinishRun()
296 {
297   // Clean generator information
298   AliDebug(1, "fGenerator->FinishRun()");
299   fGenerator->FinishRun();
300
301   //Output energy summary tables
302   AliDebug(1, "EnergySummary()");
303   ToAliDebug(1, EnergySummary());
304 }
305
306 //_______________________________________________________________________
307 void AliMC::BeginPrimary()
308 {
309   //
310   // Called  at the beginning of each primary track
311   //
312   
313   // Reset Hits info
314   ResetHits();
315   ResetTrackReferences();
316 }
317
318 //_______________________________________________________________________
319 void AliMC::PreTrack()
320 {
321   // Actions before the track's transport
322
323      TObjArray &dets = *gAlice->Modules();
324      AliModule *module;
325
326      for(Int_t i=0; i<=gAlice->GetNdets(); i++)
327        if((module = dynamic_cast<AliModule*>(dets[i])))
328          module->PreTrack();
329 }
330
331 //_______________________________________________________________________
332 void AliMC::Stepping() 
333 {
334   //
335   // Called at every step during transport
336   //
337   Int_t id = DetFromMate(gMC->CurrentMedium());
338   if (id < 0) return;
339
340
341   if ( gMC->IsNewTrack()            && 
342        gMC->TrackTime() == 0.       &&
343        fRDecayMin >= 0.             &&  
344        fRDecayMax > fRDecayMin      &&
345        gMC->TrackPid() == fDecayPdg ) 
346   {
347       FixParticleDecaytime();
348   } 
349     
350
351   
352   //
353   // --- If lego option, do it and leave 
354   if (AliSimulation::Instance()->Lego())
355     AliSimulation::Instance()->Lego()->StepManager();
356   else {
357     Int_t copy;
358     //Update energy deposition tables
359     AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
360     //
361     // write tracke reference for track which is dissapearing - MI
362
363     if (gMC->IsTrackDisappeared() && !(gMC->IsTrackAlive())) {      
364         if (gMC->Etot() > 0.05) AddTrackReference(GetCurrentTrackNumber(), 
365                                                 AliTrackReference::kDisappeared);
366         
367
368     }
369
370     //Call the appropriate stepping routine;
371     AliModule *det = dynamic_cast<AliModule*>(gAlice->Modules()->At(id));
372     if(det && det->StepManagerIsEnabled()) {
373       det->StepManager();
374     }
375   }
376 }
377
378 //_______________________________________________________________________
379 void AliMC::EnergySummary()
380 {
381   //e
382   // Print summary of deposited energy
383   //
384
385   Int_t ndep=0;
386   Float_t edtot=0;
387   Float_t ed, ed2;
388   Int_t kn, i, left, j, id;
389   const Float_t kzero=0;
390   Int_t ievent=AliRunLoader::Instance()->GetHeader()->GetEvent()+1;
391   //
392   // Energy loss information
393   if(ievent) {
394     printf("***************** Energy Loss Information per event (GEV) *****************\n");
395     for(kn=1;kn<fEventEnergy.GetSize();kn++) {
396       ed=fSummEnergy[kn];
397       if(ed>0) {
398         fEventEnergy[ndep]=kn;
399         if(ievent>1) {
400           ed=ed/ievent;
401           ed2=fSum2Energy[kn];
402           ed2=ed2/ievent;
403           ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
404         } else 
405           ed2=99;
406         fSummEnergy[ndep]=ed;
407         fSum2Energy[ndep]=TMath::Min(static_cast<Float_t>(99.),TMath::Max(ed2,kzero));
408         edtot+=ed;
409         ndep++;
410       }
411     }
412     for(kn=0;kn<(ndep-1)/3+1;kn++) {
413       left=ndep-kn*3;
414       for(i=0;i<(3<left?3:left);i++) {
415         j=kn*3+i;
416         id=Int_t (fEventEnergy[j]+0.1);
417         printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
418       }
419       printf("\n");
420     }
421     //
422     // Relative energy loss in different detectors
423     printf("******************** Relative Energy Loss per event ********************\n");
424     printf("Total energy loss per event %10.3f GeV\n",edtot);
425     for(kn=0;kn<(ndep-1)/5+1;kn++) {
426       left=ndep-kn*5;
427       for(i=0;i<(5<left?5:left);i++) {
428         j=kn*5+i;
429         id=Int_t (fEventEnergy[j]+0.1);
430         printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
431       }
432       printf("\n");
433     }
434     for(kn=0;kn<75;kn++) printf("*"); 
435     printf("\n");
436   }
437   //
438   // Reset the TArray's
439   //  fEventEnergy.Set(0);
440   //  fSummEnergy.Set(0);
441   //  fSum2Energy.Set(0);
442 }
443
444 //_____________________________________________________________________________
445 void AliMC::BeginEvent()
446 {
447   //
448   // Clean-up previous event
449   // Energy scores
450   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
451   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
452   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
453   AliDebug(1, "          BEGINNING EVENT               ");
454   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
455   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
456   AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>");
457     
458   AliRunLoader *runloader=AliRunLoader::Instance();
459
460   /*******************************/    
461   /*   Clean after eventual      */
462   /*   previous event            */
463   /*******************************/    
464
465   
466   //Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors),
467   gAlice->SetEventNrInRun(gAlice->GetEventNrInRun()+1);
468   runloader->SetEventNumber(gAlice->GetEventNrInRun());// sets new files, cleans the previous event stuff, if necessary, etc.,  
469   AliDebug(1, Form("EventNr is %d",gAlice->GetEventNrInRun()));
470      
471   fEventEnergy.Reset();  
472     // Clean detector information
473   
474   if (runloader->Stack())
475       runloader->Stack()->Reset();//clean stack -> tree is unloaded
476   else
477       runloader->MakeStack();//or make a new one
478   
479   
480   if(AliSimulation::Instance()->Lego() == 0x0)
481   { 
482       AliDebug(1, "fRunLoader->MakeTree(K)");
483       runloader->MakeTree("K");
484   }
485   
486   AliDebug(1, "gMC->SetStack(fRunLoader->Stack())");
487   gMC->SetStack(runloader->Stack());//Was in InitMC - but was moved here 
488                                      //because we don't have guarantee that 
489                                      //stack pointer is not going to change from event to event
490                          //since it bellobgs to header and is obtained via RunLoader
491   //
492   //  Reset all Detectors & kinematics & make/reset trees
493   //
494     
495   runloader->GetHeader()->Reset(AliCDBManager::Instance()->GetRun(),gAlice->GetEvNumber(),
496                                 gAlice->GetEventNrInRun());
497 //  fRunLoader->WriteKinematics("OVERWRITE");  is there any reason to rewrite here since MakeTree does so
498
499   if(AliSimulation::Instance()->Lego()) 
500   {
501       AliSimulation::Instance()->Lego()->BeginEvent();
502       return;
503   }
504   
505
506   AliDebug(1, "ResetHits()");
507   ResetHits();
508   
509   AliDebug(1, "fRunLoader->MakeTree(H)");
510   runloader->MakeTree("H");
511   
512
513
514   MakeTmpTrackRefsTree();
515   //create new branches and SetAdresses
516   TIter next(gAlice->Modules());
517   AliModule *detector;
518   while((detector = (AliModule*)next()))
519    {
520        AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName()));
521        detector->MakeBranch("H"); 
522    }
523 }
524
525 //_______________________________________________________________________
526 void AliMC::ResetHits()
527 {
528   //
529   //  Reset all Detectors hits
530   //
531   TIter next(gAlice->Modules());
532   AliModule *detector;
533   while((detector = dynamic_cast<AliModule*>(next()))) {
534      detector->ResetHits();
535   }
536 }
537
538 //_______________________________________________________________________
539 void AliMC::ResetDigits()
540 {
541   //
542   //  Reset all Detectors digits
543   //
544   TIter next(gAlice->Modules());
545   AliModule *detector;
546   while((detector = dynamic_cast<AliModule*>(next()))) {
547      detector->ResetDigits();
548   }
549 }
550
551 //_______________________________________________________________________
552 void AliMC::ResetSDigits()
553 {
554   //
555   //  Reset all Detectors digits
556   //
557   TIter next(gAlice->Modules());
558   AliModule *detector;
559   while((detector = dynamic_cast<AliModule*>(next()))) {
560      detector->ResetSDigits();
561   }
562 }
563
564 //_______________________________________________________________________
565 void AliMC::PostTrack()
566 {
567   // Posts tracks for each module
568
569   TObjArray &dets = *gAlice->Modules();
570   AliModule *module;
571   
572   for(Int_t i=0; i<=gAlice->GetNdets(); i++)
573     if((module = dynamic_cast<AliModule*>(dets[i])))
574       module->PostTrack();
575 }
576
577 //_______________________________________________________________________
578 void AliMC::FinishPrimary()
579 {
580   //
581   // Called  at the end of each primary track
582   //
583
584   AliRunLoader *runloader=AliRunLoader::Instance();
585   //  static Int_t count=0;
586   //  const Int_t times=10;
587   // This primary is finished, purify stack
588 #if ROOT_VERSION_CODE > 262152
589   if (!(gMC->SecondariesAreOrdered())) {
590       if (runloader->Stack()->ReorderKine()) RemapHits();
591   }
592 #endif
593   if (runloader->Stack()->PurifyKine()) RemapHits();
594   
595   TIter next(gAlice->Modules());
596   AliModule *detector;
597   while((detector = dynamic_cast<AliModule*>(next()))) {
598     detector->FinishPrimary();
599     AliLoader* loader = detector->GetLoader();
600     if(loader)
601      {
602        TTree* treeH = loader->TreeH();
603        if (treeH) treeH->Fill(); //can be Lego run and treeH can not exist
604      }
605   }
606
607   // Write out track references if any
608   if (fTmpTreeTR) fTmpTreeTR->Fill();
609 }
610
611 void AliMC::RemapHits()
612 {
613 //    
614 // Remaps the track labels of the hits
615     AliRunLoader *runloader=AliRunLoader::Instance();
616     AliStack* stack = runloader->Stack();
617     TList* hitLists = GetHitLists();
618     TIter next(hitLists);
619     TCollection *hitList;
620     
621     while((hitList = dynamic_cast<TCollection*>(next()))) {
622         TIter nexthit(hitList);
623         AliHit *hit;
624         while((hit = dynamic_cast<AliHit*>(nexthit()))) {
625             hit->SetTrack(stack->TrackLabel(hit->GetTrack()));
626         }
627     }
628     
629     // 
630     // This for detectors which have a special mapping mechanism
631     // for hits, such as TPC and TRD
632     //
633
634     
635     TObjArray* modules = gAlice->Modules();
636     TIter nextmod(modules);
637     AliModule *module;
638     while((module = (AliModule*) nextmod())) {
639         AliDetector* det = dynamic_cast<AliDetector*> (module);
640         if (det) det->RemapTrackHitIDs(stack->TrackLabelMap());
641     }
642     //
643     RemapTrackReferencesIDs(stack->TrackLabelMap());
644 }
645
646 //_______________________________________________________________________
647 void AliMC::FinishEvent()
648 {
649   //
650   // Called at the end of the event.
651   //
652     
653   if(AliSimulation::Instance()->Lego()) AliSimulation::Instance()->Lego()->FinishEvent();
654
655   TIter next(gAlice->Modules());
656   AliModule *detector;
657   while((detector = dynamic_cast<AliModule*>(next()))) {
658     detector->FinishEvent();
659   }
660
661   //Update the energy deposit tables
662   Int_t i;
663   for(i=0;i<fEventEnergy.GetSize();i++) 
664    {
665     fSummEnergy[i]+=fEventEnergy[i];
666     fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
667    }
668
669   AliRunLoader *runloader=AliRunLoader::Instance();
670
671   AliHeader* header = runloader->GetHeader();
672   AliStack* stack = runloader->Stack();
673   if ( (header == 0x0) || (stack == 0x0) )
674    {//check if we got header and stack. If not cry and exit aliroot
675     AliFatal("Can not get the stack or header from LOADER");
676     return;//never reached
677    }  
678   // Update Header information 
679   header->SetNprimary(stack->GetNprimary());
680   header->SetNtrack(stack->GetNtrack());  
681
682   // Write out the kinematics
683   if (!AliSimulation::Instance()->Lego()) stack->FinishEvent();
684
685   // Synchronize the TreeTR with TreeK
686   if (fTmpTreeTR) ReorderAndExpandTreeTR();
687    
688   // Write out the event Header information
689   TTree* treeE = runloader->TreeE();
690   if (treeE) 
691    {
692       header->SetStack(stack);
693       treeE->Fill();
694    }
695   else
696    {
697     AliError("Can not get TreeE from RL");
698    }
699   
700   if(AliSimulation::Instance()->Lego() == 0x0)
701    {
702      runloader->WriteKinematics("OVERWRITE");
703      runloader->WriteTrackRefs("OVERWRITE");
704      runloader->WriteHits("OVERWRITE");
705    }
706    
707   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
708   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
709   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
710   AliDebug(1, "          FINISHING EVENT               ");
711   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
712   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
713   AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<");
714 }
715
716 //_______________________________________________________________________
717 void AliMC::Init()
718 {
719   // MC initialization
720
721    //=================Create Materials and geometry
722    gMC->Init();
723   // Set alignable volumes for the whole geometry (with old root)
724 #if ROOT_VERSION_CODE < 331527
725   SetAllAlignableVolumes();
726 #endif
727    //Read the cuts for all materials
728    ReadTransPar();
729    //Build the special IMEDIA table
730    MediaTable();
731
732    //Compute cross-sections
733    gMC->BuildPhysics();
734    
735    //Initialise geometry deposition table
736    fEventEnergy.Set(gMC->NofVolumes()+1);
737    fSummEnergy.Set(gMC->NofVolumes()+1);
738    fSum2Energy.Set(gMC->NofVolumes()+1);
739
740    // Register MC in configuration 
741    AliConfig::Instance()->Add(gMC);
742
743 }
744
745 //_______________________________________________________________________
746 void AliMC::MediaTable()
747 {
748   //
749   // Built media table to get from the media number to
750   // the detector id
751   //
752
753   Int_t kz, nz, idt, lz, i, k, ind;
754   //  Int_t ibeg;
755   TObjArray &dets = *gAlice->Detectors();
756   AliModule *det;
757   Int_t ndets=gAlice->GetNdets();
758   //
759   // For all detectors
760   for (kz=0;kz<ndets;kz++) {
761     // If detector is defined
762     if((det=dynamic_cast<AliModule*>(dets[kz]))) {
763         TArrayI &idtmed = *(det->GetIdtmed()); 
764         for(nz=0;nz<100;nz++) {
765             
766         // Find max and min material number
767         if((idt=idtmed[nz])) {
768           det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
769           det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
770         }
771       }
772       if(det->LoMedium() > det->HiMedium()) {
773         det->LoMedium() = 0;
774         det->HiMedium() = 0;
775       } else {
776         if(det->HiMedium() > fImedia->GetSize()) {
777           AliError(Form("Increase fImedia from %d to %d",
778                         fImedia->GetSize(),det->HiMedium()));
779           return;
780         }
781         // Tag all materials in rage as belonging to detector kz
782         for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
783           (*fImedia)[lz]=kz;
784         }
785       }
786     }
787   }
788   //
789   // Print summary table
790   AliInfo("Tracking media ranges:");
791   ToAliInfo(
792   for(i=0;i<(ndets-1)/6+1;i++) {
793     for(k=0;k< (6<ndets-i*6?6:ndets-i*6);k++) {
794       ind=i*6+k;
795       det=dynamic_cast<AliModule*>(dets[ind]);
796       if(det)
797         printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
798                det->HiMedium());
799       else
800         printf(" %6s: %3d -> %3d;","NULL",0,0);
801     }
802     printf("\n");
803   }
804   )
805 }
806
807 //_______________________________________________________________________
808 void AliMC::ReadTransPar()
809 {
810   //
811   // Read filename to set the transport parameters
812   //
813
814
815   const Int_t kncuts=10;
816   const Int_t knflags=12;
817   const Int_t knpars=kncuts+knflags;
818   const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
819                                "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
820                                "BREM","COMP","DCAY","DRAY","HADR","LOSS",
821                                "MULS","PAIR","PHOT","RAYL","STRA"};
822   char line[256];
823   char detName[7];
824   char* filtmp;
825   Float_t cut[kncuts];
826   Int_t flag[knflags];
827   Int_t i, itmed, iret, ktmed, kz;
828   FILE *lun;
829   //
830   // See whether the file is there
831   filtmp=gSystem->ExpandPathName(fTransParName.Data());
832   lun=fopen(filtmp,"r");
833   delete [] filtmp;
834   if(!lun) {
835     AliWarning(Form("File %s does not exist!",fTransParName.Data()));
836     return;
837   }
838   //
839   while(1) {
840     // Initialise cuts and flags
841     for(i=0;i<kncuts;i++) cut[i]=-99;
842     for(i=0;i<knflags;i++) flag[i]=-99;
843     itmed=0;
844     for(i=0;i<256;i++) line[i]='\0';
845     // Read up to the end of line excluded
846     iret=fscanf(lun,"%[^\n]",line);
847     if(iret<0) {
848       //End of file
849       fclose(lun);
850       return;
851     }
852     // Read the end of line
853     fscanf(lun,"%*c");
854     if(!iret) continue;
855     if(line[0]=='*') continue;
856     // Read the numbers
857     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 %d",
858                 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
859                 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
860                 &flag[8],&flag[9],&flag[10],&flag[11]);
861     if(!iret) continue;
862     if(iret<0) {
863       //reading error
864       AliWarning(Form("Error reading file %s",fTransParName.Data()));
865       continue;
866     }
867     // Check that the module exist
868     AliModule *mod = gAlice->GetModule(detName);
869     if(mod) {
870       // Get the array of media numbers
871       TArrayI &idtmed = *mod->GetIdtmed();
872       // Check that the tracking medium code is valid
873       if(0<=itmed && itmed < 100) {
874         ktmed=idtmed[itmed];
875         if(!ktmed) {
876           AliWarning(Form("Invalid tracking medium code %d for %s",itmed,mod->GetName()));
877           continue;
878         }
879         // Set energy thresholds
880         for(kz=0;kz<kncuts;kz++) {
881           if(cut[kz]>=0) {
882             AliDebug(2, Form("%-6s set to %10.3E for tracking medium code %4d for %s",
883                              kpars[kz],cut[kz],itmed,mod->GetName()));
884             gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
885           }
886         }
887         // Set transport mechanisms
888         for(kz=0;kz<knflags;kz++) {
889           if(flag[kz]>=0) {
890             AliDebug(2, Form("%-6s set to %10d for tracking medium code %4d for %s",
891                              kpars[kncuts+kz],flag[kz],itmed,mod->GetName()));
892             gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
893           }
894         }
895       } else {
896         AliWarning(Form("Invalid medium code %d",itmed));
897         continue;
898       }
899     } else {
900       AliDebug(1, Form("%s not present",detName));
901       continue;
902     }
903   }
904 }
905
906 //_______________________________________________________________________
907 void AliMC::SetTransPar(const char *filename)
908 {
909   //
910   // Sets the file name for transport parameters
911   //
912   fTransParName = filename;
913 }
914
915 //_______________________________________________________________________
916 void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
917 {
918   //
919   //  Add a hit to detector id
920   //
921   TObjArray &dets = *gAlice->Modules();
922   if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddHit(track,vol,hits);
923 }
924
925 //_______________________________________________________________________
926 void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
927 {
928   //
929   // Add digit to detector id
930   //
931   TObjArray &dets = *gAlice->Modules();
932   if(dets[id]) dynamic_cast<AliModule*>(dets[id])->AddDigit(tracks,digits);
933 }
934
935 //_______________________________________________________________________
936 Int_t AliMC::GetCurrentTrackNumber() const {
937   //
938   // Returns current track
939   //
940   return AliRunLoader::Instance()->Stack()->GetCurrentTrackNumber();
941 }
942
943 //_______________________________________________________________________
944 void AliMC::DumpPart (Int_t i) const
945 {
946   //
947   // Dumps particle i in the stack
948   //
949   AliRunLoader * runloader = AliRunLoader::Instance();
950    if (runloader->Stack())
951     runloader->Stack()->DumpPart(i);
952 }
953
954 //_______________________________________________________________________
955 void AliMC::DumpPStack () const
956 {
957   //
958   // Dumps the particle stack
959   //
960   AliRunLoader * runloader = AliRunLoader::Instance();
961    if (runloader->Stack())
962     runloader->Stack()->DumpPStack();
963 }
964
965 //_______________________________________________________________________
966 Int_t AliMC::GetNtrack() const {
967   //
968   // Returns number of tracks in stack
969   //
970   Int_t ntracks = -1;
971   AliRunLoader * runloader = AliRunLoader::Instance();
972    if (runloader->Stack())
973      ntracks = runloader->Stack()->GetNtrack();
974    return ntracks;
975 }
976
977 //_______________________________________________________________________
978 Int_t AliMC::GetPrimary(Int_t track) const
979 {
980   //
981   // return number of primary that has generated track
982   //
983   Int_t nprimary = -999;
984   AliRunLoader * runloader = AliRunLoader::Instance();
985   if (runloader->Stack())
986     nprimary = runloader->Stack()->GetPrimary(track);
987   return nprimary;
988 }
989  
990 //_______________________________________________________________________
991 TParticle* AliMC::Particle(Int_t i) const
992 {
993   // Returns the i-th particle from the stack taking into account
994   // the remaping done by PurifyKine
995   AliRunLoader * runloader = AliRunLoader::Instance();
996   if (runloader)
997    if (runloader->Stack())
998     return runloader->Stack()->Particle(i);
999   return 0x0;   
1000 }
1001
1002 //_______________________________________________________________________
1003 const TObjArray* AliMC::Particles() const {
1004   //
1005   // Returns pointer to Particles array
1006   //
1007   AliRunLoader * runloader = AliRunLoader::Instance();
1008   if (runloader)
1009    if (runloader->Stack())
1010     return runloader->Stack()->Particles();
1011   return 0x0;
1012 }
1013
1014 //_______________________________________________________________________
1015 void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom,
1016                       const Float_t *vpos, const Float_t *polar, Float_t tof,
1017                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
1018
1019 // Delegate to stack
1020 //
1021   AliRunLoader * runloader = AliRunLoader::Instance();
1022   if (runloader)
1023     if (runloader->Stack())
1024       runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof,
1025                                     mech, ntr, weight, is);
1026 }
1027
1028 //_______________________________________________________________________
1029 void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg,
1030                       Double_t px, Double_t py, Double_t pz, Double_t e,
1031                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
1032                       Double_t polx, Double_t poly, Double_t polz,
1033                       TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const
1034
1035   // Delegate to stack
1036   //
1037   AliRunLoader * runloader = AliRunLoader::Instance();
1038   if (runloader)
1039     if (runloader->Stack())
1040       runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
1041                                     polx, poly, polz, mech, ntr, weight, is);
1042 }
1043
1044 //_______________________________________________________________________
1045 void AliMC::SetHighWaterMark(Int_t nt) const
1046 {
1047     //
1048     // Set high water mark for last track in event
1049   AliRunLoader * runloader = AliRunLoader::Instance();
1050   if (runloader)
1051     if (runloader->Stack())
1052       runloader->Stack()->SetHighWaterMark(nt);
1053 }
1054
1055 //_______________________________________________________________________
1056 void AliMC::KeepTrack(Int_t track) const
1057
1058   //
1059   // Delegate to stack
1060   //
1061   AliRunLoader * runloader = AliRunLoader::Instance();
1062   if (runloader)
1063     if (runloader->Stack())
1064       runloader->Stack()->KeepTrack(track);
1065 }
1066  
1067 //_______________________________________________________________________
1068 void AliMC::FlagTrack(Int_t track) const
1069 {
1070   // Delegate to stack
1071   //
1072   AliRunLoader * runloader = AliRunLoader::Instance();
1073   if (runloader)
1074     if (runloader->Stack())
1075       runloader->Stack()->FlagTrack(track);
1076 }
1077
1078 //_______________________________________________________________________
1079 void AliMC::SetCurrentTrack(Int_t track) const
1080
1081   //
1082   // Set current track number
1083   //
1084   AliRunLoader * runloader = AliRunLoader::Instance();
1085   if (runloader)
1086     if (runloader->Stack())
1087       runloader->Stack()->SetCurrentTrack(track); 
1088 }
1089
1090 //_______________________________________________________________________
1091 AliTrackReference*  AliMC::AddTrackReference(Int_t label, Int_t id) 
1092 {
1093   //
1094   // add a trackrefernce to the list
1095   Int_t primary = GetPrimary(label);
1096   Particle(primary)->SetBit(kKeepBit);
1097
1098   Int_t nref = fTmpTrackReferences.GetEntriesFast();
1099   return new(fTmpTrackReferences[nref]) AliTrackReference(label, id);
1100 }
1101
1102
1103
1104 //_______________________________________________________________________
1105 void AliMC::ResetTrackReferences()
1106 {
1107   //
1108   //  Reset all  references
1109   //
1110     fTmpTrackReferences.Clear();
1111 }
1112
1113 //_______________________________________________________________________
1114 void AliMC::RemapTrackReferencesIDs(Int_t *map)
1115 {
1116   // 
1117   // Remapping track reference
1118   // Called at finish primary
1119   //
1120     
1121   Int_t nEntries = fTmpTrackReferences.GetEntries();
1122   for (Int_t i=0; i < nEntries; i++){
1123       AliTrackReference * ref = dynamic_cast<AliTrackReference*>(fTmpTrackReferences.UncheckedAt(i));
1124       if (ref) {
1125           Int_t newID = map[ref->GetTrack()];
1126           if (newID>=0) ref->SetTrack(newID);
1127           else {
1128               ref->SetBit(kNotDeleted,kFALSE);
1129               fTmpTrackReferences.RemoveAt(i);  
1130           }      
1131       } // if ref
1132   }
1133   fTmpTrackReferences.Compress();
1134 }
1135
1136 //_______________________________________________________________________
1137 void AliMC::FixParticleDecaytime()
1138 {
1139     //
1140     // Fix the particle decay time according to rmin and rmax for decays
1141     //
1142
1143     TLorentzVector p;
1144     gMC->TrackMomentum(p);
1145     Double_t tmin, tmax;
1146     Double_t b;
1147
1148     // Transverse velocity 
1149     Double_t vt    = p.Pt() / p.E();
1150     
1151     if ((b = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField()) > 0.) {     // [kG]
1152
1153         // Radius of helix
1154         
1155         Double_t rho   = p.Pt() / 0.0003 / b; // [cm]
1156         
1157         // Revolution frequency
1158         
1159         Double_t omega = vt / rho;
1160         
1161         // Maximum and minimum decay time
1162         //
1163         // Check for curlers first
1164         const Double_t kOvRhoSqr2 = 1./(rho*TMath::Sqrt(2.));
1165         if (fRDecayMax * kOvRhoSqr2 > 1.) return;
1166         
1167         //
1168  
1169         tmax  = TMath::ACos((1.-fRDecayMax*kOvRhoSqr2)*(1.+fRDecayMax*kOvRhoSqr2)) / omega;   // [ct]
1170         tmin  = TMath::ACos((1.-fRDecayMin*kOvRhoSqr2)*(1.+fRDecayMin*kOvRhoSqr2)) / omega;   // [ct]
1171     } else {
1172         tmax =  fRDecayMax / vt;                                                      // [ct] 
1173         tmin =  fRDecayMin / vt;                                                      // [ct]
1174     }
1175     
1176     //
1177     // Dial t using the two limits
1178     Double_t t = tmin + (tmax - tmin) * gRandom->Rndm();                              // [ct]
1179     //
1180     //
1181     // Force decay time in transport code
1182     //
1183     gMC->ForceDecayTime(t / 2.99792458e10);
1184 }
1185
1186 void AliMC::MakeTmpTrackRefsTree()
1187 {
1188     // Make the temporary track reference tree
1189     fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
1190     fTmpTreeTR = new TTree("TreeTR", "Track References");
1191     TClonesArray* pRef = &fTmpTrackReferences;
1192     fTmpTreeTR->Branch("TrackReferences", &pRef, 4000);
1193 }
1194
1195 //_______________________________________________________________________
1196 void AliMC::ReorderAndExpandTreeTR()
1197 {
1198 //
1199 //  Reorder and expand the temporary track reference tree in order to match the kinematics tree
1200 //
1201
1202     AliRunLoader *rl = AliRunLoader::Instance();
1203 //
1204 //  TreeTR
1205     AliDebug(1, "fRunLoader->MakeTrackRefsContainer()");
1206     rl->MakeTrackRefsContainer(); 
1207     TTree * treeTR = rl->TreeTR();
1208     if (treeTR){
1209         // make branch for central track references
1210         TBranch *branch;
1211         TClonesArray* pRef = &fTrackReferences;
1212         branch = treeTR->Branch("TrackReferences", &pRef);
1213     }
1214
1215     AliStack* stack  = rl->Stack();
1216     Int_t np = stack->GetNprimary();
1217     Int_t nt = fTmpTreeTR->GetEntries();
1218     //
1219     // Loop over tracks and find the secondaries with the help of the kine tree
1220     Int_t ifills = 0;
1221     Int_t it = 0;
1222     for (Int_t ip = np - 1; ip > -1; ip--) {
1223         TParticle *part = stack->Particle(ip);
1224         //printf("Particle %5d %5d %5d %5d %5d \n", ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), part->GetLastDaughter());
1225         
1226         // Skip primaries that have not been transported
1227         Int_t dau1  = part->GetFirstDaughter();
1228         Int_t dau2  = -1;
1229         if (!part->TestBit(kTransportBit)) continue;
1230         //
1231         fTmpTreeTR->GetEntry(it++);
1232         Int_t nh = fTmpTrackReferences.GetEntries();
1233         // Determine range of secondaries produced by this primary
1234         if (dau1 > -1) {
1235             Int_t inext = ip - 1;
1236             while (dau2 < 0) {
1237                 if (inext >= 0) {
1238                     part = stack->Particle(inext);
1239                     dau2 =  part->GetFirstDaughter();
1240                     if (!(part->TestBit(kTransportBit)) || dau2 == -1 || dau2 < np) {
1241 //                  if (dau2 == -1 || dau2 < np) {
1242                         dau2 = -1;
1243                     } else {
1244                         dau2--;
1245                     }
1246                 } else {
1247                     dau2 = stack->GetNtrack() - 1;
1248                 }
1249                 inext--;
1250             } // find upper bound
1251         }  // dau2 < 0
1252 //      printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
1253         // 
1254         // Loop over reference hits and find secondary label
1255         for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
1256             for (Int_t ih = 0; ih < nh; ih++) {
1257                 AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences.At(ih);
1258                 Int_t label = tr->Label();
1259                 // Skip primaries
1260                 if (label == ip) continue;
1261                 if (label > dau2 || label < dau1) 
1262                     AliWarning(Form("Track Reference Label out of range !: %5d %5d %5d \n", label, dau1, dau2));
1263                 if (label == id) {
1264                     // secondary found
1265                     Int_t nref =  fTrackReferences.GetEntriesFast();
1266                     new(fTrackReferences[nref]) AliTrackReference(*tr);
1267                 }
1268             } // hits
1269             treeTR->Fill();
1270             fTrackReferences.Clear();
1271             ifills++;
1272         } // daughters
1273     } // tracks
1274     //
1275     // Now loop again and write the primaries
1276     it = nt - 1;
1277     for (Int_t ip = 0; ip < np; ip++) {
1278         TParticle* part = stack->Particle(ip);
1279 //      if ((part->GetFirstDaughter() == -1 && part->GetStatusCode() <= 1) || part->GetFirstDaughter() >= np) 
1280         if (part->TestBit(kTransportBit))
1281         {
1282             // Skip particles that have not been transported
1283             fTmpTreeTR->GetEntry(it--);
1284             Int_t nh = fTmpTrackReferences.GetEntries();
1285             // 
1286             // Loop over reference hits and find primary labels
1287             for (Int_t ih = 0; ih < nh; ih++) {
1288                 AliTrackReference* tr = (AliTrackReference*)  fTmpTrackReferences.At(ih);
1289                 Int_t label = tr->Label();
1290                 if (label == ip) {
1291                     Int_t nref = fTrackReferences.GetEntriesFast();
1292                     new(fTrackReferences[nref]) AliTrackReference(*tr);
1293                 }
1294             } 
1295         }
1296         treeTR->Fill();
1297         fTrackReferences.Clear();
1298         ifills++;
1299     } // tracks
1300     // Check
1301     if (ifills != stack->GetNtrack()) 
1302         AliWarning(Form("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, stack->GetNtrack()));
1303 //
1304 //  Clean-up
1305     delete fTmpTreeTR;
1306     fTmpFileTR->Close();
1307     delete fTmpFileTR;
1308     fTmpTrackReferences.Clear();
1309     gSystem->Exec("rm -rf TrackRefsTmp.root");
1310 }
1311