b86eff3816cc61e73b497583554c568ec5257663
[u/mrichter/AliRoot.git] / STEER / STEER / AliTransportMonitor.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 // Class that can be plugged in the simulation to monitor transport timing per 
17 // particle for each geometry volume.
18 //
19 //  andrei.gheata@cern.ch 
20
21 #include "AliTransportMonitor.h"
22
23 #include "TDatabasePDG.h"
24 #include "TCanvas.h"
25 #include "TMath.h"
26 #include "TFile.h"
27 #include "TROOT.h"
28 #include "THashList.h"
29 #include "TH2F.h"
30 #include "TGeoManager.h"
31 #include "AliPDG.h"
32
33 ClassImp(AliTransportMonitor)
34 ClassImp(AliTransportMonitor::AliTransportMonitorVol)
35 ClassImp(AliTransportMonitor::AliTransportMonitorVol::AliPMonData)
36
37 typedef AliTransportMonitor::AliTransportMonitorVol::AliPMonData PMonData;
38
39 //______________________________________________________________________________
40 AliTransportMonitor::AliTransportMonitorVol::AliTransportMonitorVol()
41                     :TNamed(),
42                      fNtypes(0),
43                      fTotalTime(0),
44                      fPData(0),
45                      fTimeRZ(0),
46                      fParticles()
47 {
48 // Default constructor
49 }
50
51 //______________________________________________________________________________
52 AliTransportMonitor::AliTransportMonitorVol::~AliTransportMonitorVol()
53 {
54 // Destructor
55   delete [] fPData;
56   delete fTimeRZ;
57 }
58
59 //______________________________________________________________________________
60 void AliTransportMonitor::AliTransportMonitorVol::StepInfo(
61                                     Int_t pdg,
62                                     Double_t energy, 
63                                     Double_t dt,
64                                     Double_t x, Double_t y, Double_t z)
65 {
66 // This method is called at each N steps to store timing info.
67   PMonData &data = GetPMonData(pdg);
68   data.fEdt += energy*dt;
69   data.fTime += dt;
70   fTotalTime += dt;
71   if (!fTimeRZ) {
72     Bool_t status = TH1::AddDirectoryStatus();
73     TH1::AddDirectory(kFALSE);
74     fTimeRZ = new TH2F("h2rz", "", 100, -5000., 5000, 100, 0., 1000.);
75     TH1::AddDirectory(status);
76   }
77   Double_t r = TMath::Sqrt(x*x+y*y);
78   fTimeRZ->Fill(z,r,dt);
79 }
80
81 //______________________________________________________________________________
82 PMonData &AliTransportMonitor::AliTransportMonitorVol::GetPMonData(Int_t pdg)
83 {
84 // Retrieve stored monitoring object for a given pdg type. If not existing 
85 // create one.
86
87   // The object could have been retrieved from file, in which case we have to 
88   // build the map.
89   PMonData *data;
90   if (fNtypes) {
91     if (fParticles.empty()) {
92       for (Int_t i=0; i<fNtypes; i++) {
93         data = &fPData[i];
94         fParticles[i] = data->fPDG;
95       }
96     }
97     ParticleMapIt_t it = fParticles.find(pdg);
98     if (it == fParticles.end()) {
99       data = &fPData[fNtypes]; 
100       data->fPDG = pdg;
101       fParticles[pdg] = fNtypes++;
102     } else {
103       data = &fPData[it->second];
104     }  
105   } else {
106      TDatabasePDG *pdgDB = TDatabasePDG::Instance();
107      if (!pdgDB->ParticleList()) AliPDG::AddParticlesToPdgDataBase();
108      Int_t size = pdgDB->ParticleList()->GetSize();
109      fPData = new PMonData[size];
110      data = &fPData[fNtypes];
111      data->fPDG = pdg;
112      fParticles[pdg] = fNtypes++;
113   }
114   return *data;
115 }
116
117 void AliTransportMonitor::AliTransportMonitorVol::Merge(AliTransportMonitorVol* volM) 
118 {
119   // Merging
120   fTotalTime = (fTotalTime + volM->GetTotalTime());
121   if (fTimeRZ && volM->GetHistogram()) {
122     fTimeRZ->Add(volM->GetHistogram()); 
123   } else if (volM->GetHistogram()) {
124     fTimeRZ = (TH2F*)(volM->GetHistogram()->Clone());
125   }
126
127   Int_t ntypes = volM->GetNtypes();
128   Int_t ntypes2 = GetNtypes();
129   for (Int_t i = 0; i < ntypes; i++) {
130     Int_t pdg = volM->GetPDG(i);
131     //PMonData &data  = GetPMonData(pdg);
132     //PMonData &dataM = volM->GetPMonData(pdg);
133     //data.fEdt  += dataM.fEdt;
134     //data.fTime += dataM.fTime;
135   }
136 }
137
138 ClassImp(AliTransportMonitor)
139
140 //______________________________________________________________________________
141 AliTransportMonitor::AliTransportMonitor()
142                     :TObject(),
143                      fTotalTime(0),
144                      fTimer(),
145                      fVolumeMon(0)
146 {
147 // Default constructor
148 }
149
150 //______________________________________________________________________________
151 AliTransportMonitor::AliTransportMonitor(Int_t nvolumes)
152                     :TObject(),
153                      fTotalTime(0),
154                      fTimer(),
155                      fVolumeMon(0)
156 {
157 // Default constructor
158   fVolumeMon = new TObjArray(nvolumes);
159   fVolumeMon->SetOwner();
160   for (Int_t i=0; i<nvolumes; i++) {
161     AliTransportMonitorVol *volMon = new AliTransportMonitorVol();
162     if (gGeoManager) volMon->SetName(gGeoManager->GetListOfUVolumes()->At(i)->GetName());
163     fVolumeMon->Add(volMon);
164   }   
165 }
166
167 //______________________________________________________________________________
168 AliTransportMonitor::~AliTransportMonitor()
169 {
170 // Destructor
171   delete fVolumeMon;
172 }
173
174 //______________________________________________________________________________
175 void AliTransportMonitor::Print(Option_t *volName) const
176 {
177 // Inspect the timing statistics for a single volume or for all the setup
178   Int_t uid = -1;
179   Int_t ntotal = 0;
180   if (!fVolumeMon || !(ntotal=fVolumeMon->GetEntriesFast())) {
181      Info("Inspect", "Transport monitor is empty !");
182      return;
183   }   
184   if (strlen(volName)) {
185     TString svname = volName;
186     Int_t i = 0;
187     for (i=1; i<ntotal; i++) if (svname == fVolumeMon->At(i)->GetName()) break;
188     if (i==ntotal) {
189        Error("Inspect", "No monitoring info stored for volume %s", volName);
190        return;
191     }
192     uid = i;
193     AliTransportMonitorVol *volMon = (AliTransportMonitorVol*)fVolumeMon->At(uid);
194     Int_t ntypes = volMon->GetNtypes();
195     if (!ntypes) {
196       Info("Inspect", "No particles crossed volume %s", volName);
197       return;
198     }  
199     Double_t *timeperpart = new Double_t[ntypes];
200     Int_t *isort = new Int_t[ntypes];
201     Double_t timepervol = 0.;
202     for (i=0; i<ntypes; i++) {
203       timeperpart[i] = volMon->GetTime(i); 
204       timepervol += timeperpart[i];
205     }
206     printf("Volume %s: Transport time: %g%% of %g [s]\n", volMon->GetName(), 100.*timepervol/fTotalTime, fTotalTime);
207     TMath::Sort(ntypes, timeperpart, isort, kTRUE);
208     TString particle;
209     TDatabasePDG *pdgDB = TDatabasePDG::Instance();    
210     if (!pdgDB->ParticleList()) AliPDG::AddParticlesToPdgDataBase();
211     for (i=0; i<ntypes; i++)  {
212        timeperpart[i] /=  timepervol;
213        particle = pdgDB->GetParticle(volMon->GetPDG(isort[i]))->GetName();
214        printf("   %s: %g%%  mean energy: %g\n", particle.Data(), 100.*timeperpart[i], volMon->GetEmed(i));
215     }
216     if (volMon->GetHistogram()) {
217       TCanvas *c1 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("crz");
218       if (!c1) c1 = new TCanvas("crz");
219       c1->cd();
220       volMon->GetHistogram()->GetXaxis()->SetTitle("z [cm]");
221       volMon->GetHistogram()->GetYaxis()->SetTitle("r [cm]");
222       volMon->GetHistogram()->SetTitle(Form("RZ plot weighted by time spent in %s",volMon->GetName()));
223       volMon->GetHistogram()->Draw();
224     }  
225     return;
226   }
227   // General view
228   TIter next(fVolumeMon);
229   AliTransportMonitorVol *volMon;
230   Int_t ncrossed = 0;
231   TH1F *hnames = new TH1F("volume timing", "relative volume timing", 3,0,3);
232   hnames->SetStats(0);
233   hnames->SetFillColor(38);
234   hnames->SetBit(TH1::kCanRebin);
235   while ((volMon=(AliTransportMonitorVol*)next())) {
236     if (volMon->GetNtypes()) {
237       hnames->Fill(volMon->GetName(), volMon->GetTotalTime());
238       ncrossed++;
239     }
240   }
241   
242   hnames->LabelsDeflate();
243   hnames->GetXaxis()->LabelsOption(">");
244   
245   TCanvas *c = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("cvol_timing");
246   if (!c) c = new TCanvas("cvol_timing");
247   c->cd();
248   c->SetLogy();
249   hnames->Draw();
250   
251   printf("=============================================================================\n");
252   printf("Effective transport time:  %6.2f minutes\n", fTotalTime/60.);
253   printf("Number of crossed volumes: %d from %d\n", ncrossed, fVolumeMon->GetEntriesFast());
254   printf("=============================================================================\n");  
255 }     
256
257 //______________________________________________________________________________
258 void AliTransportMonitor::DummyStep()
259 {
260 // Reset timer for zero-length steps
261    fTimer.Stop();
262    fTimer.Start(kTRUE);
263 }   
264
265 //______________________________________________________________________________
266 void AliTransportMonitor::StepInfo( Int_t volId,
267                                     Int_t pdg,
268                                     Double_t energy, 
269                                     Double_t x, Double_t y, Double_t z)
270 {
271 // This method is called at each N steps to store timing info.
272   fTimer.Stop();
273   Double_t dt = fTimer.RealTime();
274   fTotalTime += dt;
275   AliTransportMonitorVol *volMon = (AliTransportMonitorVol*)fVolumeMon->At(volId);
276   volMon->StepInfo(pdg,energy,dt,x,y,z);
277   fTimer.Start(kTRUE);
278 }
279
280 //______________________________________________________________________________
281 void AliTransportMonitor::Start()
282 {
283 // Start collecting timing information
284   if (fTotalTime > 0) {
285     Info("Start", "Cannot start twice");
286     return;
287   }
288   fTimer.Start(kTRUE);
289 }
290
291 //______________________________________________________________________________
292 void AliTransportMonitor::Stop()
293 {
294 // Stop the global timer
295   fTimer.Stop();
296 }
297
298 //______________________________________________________________________________
299 void AliTransportMonitor::Export(const char *fname)
300 {
301 // Export information to file
302   TFile *file = TFile::Open(fname, "RECREATE");
303   Write();
304   file->Write();
305   file->Close();
306 }  
307 //______________________________________________________________________________
308 void AliTransportMonitor::Merge(AliTransportMonitor* mergeMon)
309 {
310   // merge with monitor 
311   Int_t n = fVolumeMon->GetEntriesFast();
312   TObjArray* mergeVols = mergeMon->GetVolumes();
313   fTotalTime = 0;
314   for (Int_t i = 0; i < n; i++)
315     {
316       AliTransportMonitorVol *volMon1 = (AliTransportMonitorVol*)fVolumeMon->At(i);      
317       AliTransportMonitorVol *volMon2 = (AliTransportMonitorVol*)mergeVols->At(i);      
318       volMon1->Merge(volMon2);
319       fTotalTime += (volMon1->GetTotalTime());
320     }
321 }
322 //______________________________________________________________________________
323 AliTransportMonitor *AliTransportMonitor::Import(const char *fname)
324 {
325 // Import information from a file
326   TFile *file = TFile::Open(fname);
327   if (!file) {
328     ::Error("Import", "File %s could not be opened", fname);
329     return 0;
330   }
331   AliTransportMonitor *mon = (AliTransportMonitor *)file->Get("AliTransportMonitor");
332   if (!mon) {
333     ::Error("Import", "No AliTransportMonitor object found n file %s", fname);
334     return 0;
335   }
336   AliPDG::AddParticlesToPdgDataBase();
337   return mon;
338 }
339
340