Improved merging.
[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   //
120   // Merging
121   //
122   fTotalTime = (fTotalTime + volM->GetTotalTime());
123   if (fTimeRZ && volM->GetHistogram()) {
124     fTimeRZ->Add(volM->GetHistogram()); 
125   } else if (volM->GetHistogram()) {
126     fTimeRZ = (TH2F*)(volM->GetHistogram()->Clone());
127   }
128
129   Int_t ntypes = volM->GetNtypes();
130   for (Int_t i = 0; i < ntypes; i++) {
131     Int_t pdg = volM->GetPDG(i);
132      PMonData &data  = GetPMonData(pdg);
133      data.fEdt  += (volM->GetEmed(i) * volM->GetTotalTime());
134      data.fTime += (volM->GetTime(i));
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 %g [s]\n", volMon->GetName(), 100.*timepervol/fTotalTime, fTotalTime,
207            volMon->GetTotalTime());
208     TMath::Sort(ntypes, timeperpart, isort, kTRUE);
209     TString particle;
210     TDatabasePDG *pdgDB = TDatabasePDG::Instance();    
211     if (!pdgDB->ParticleList()) AliPDG::AddParticlesToPdgDataBase();
212     for (i=0; i<ntypes; i++)  {
213        timeperpart[i] /=  timepervol;
214        particle = pdgDB->GetParticle(volMon->GetPDG(isort[i]))->GetName();
215        printf("   %s: %g%%  mean energy: %g\n", particle.Data(), 100.*timeperpart[i], volMon->GetEmed(i));
216     }
217     if (volMon->GetHistogram()) {
218       TCanvas *c1 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("crz");
219       if (!c1) c1 = new TCanvas("crz");
220       c1->cd();
221       volMon->GetHistogram()->GetXaxis()->SetTitle("z [cm]");
222       volMon->GetHistogram()->GetYaxis()->SetTitle("r [cm]");
223       volMon->GetHistogram()->SetTitle(Form("RZ plot weighted by time spent in %s",volMon->GetName()));
224       volMon->GetHistogram()->Draw();
225     }  
226     return;
227   }
228   // General view
229   TIter next(fVolumeMon);
230   AliTransportMonitorVol *volMon;
231   Int_t ncrossed = 0;
232   TH1F *hnames = new TH1F("volume timing", "relative volume timing", 3,0,3);
233   hnames->SetStats(0);
234   hnames->SetFillColor(38);
235   hnames->SetBit(TH1::kCanRebin);
236   while ((volMon=(AliTransportMonitorVol*)next())) {
237     if (volMon->GetNtypes()) {
238       hnames->Fill(volMon->GetName(), volMon->GetTotalTime());
239       ncrossed++;
240     }
241   }
242   
243   hnames->LabelsDeflate();
244   hnames->GetXaxis()->LabelsOption(">");
245   
246   TCanvas *c = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("cvol_timing");
247   if (!c) c = new TCanvas("cvol_timing");
248   c->cd();
249   c->SetLogy();
250   hnames->Draw();
251   
252   printf("=============================================================================\n");
253   printf("Effective transport time:  %6.2f minutes\n", fTotalTime/60.);
254   printf("Number of crossed volumes: %d from %d\n", ncrossed, fVolumeMon->GetEntriesFast());
255   printf("=============================================================================\n");  
256 }     
257
258 //______________________________________________________________________________
259 void AliTransportMonitor::DummyStep()
260 {
261 // Reset timer for zero-length steps
262    fTimer.Stop();
263    fTimer.Start(kTRUE);
264 }   
265
266 //______________________________________________________________________________
267 void AliTransportMonitor::StepInfo( Int_t volId,
268                                     Int_t pdg,
269                                     Double_t energy, 
270                                     Double_t x, Double_t y, Double_t z)
271 {
272 // This method is called at each N steps to store timing info.
273   fTimer.Stop();
274   Double_t dt = fTimer.RealTime();
275   fTotalTime += dt;
276   AliTransportMonitorVol *volMon = (AliTransportMonitorVol*)fVolumeMon->At(volId);
277   volMon->StepInfo(pdg,energy,dt,x,y,z);
278   fTimer.Start(kTRUE);
279 }
280
281 //______________________________________________________________________________
282 void AliTransportMonitor::Start()
283 {
284 // Start collecting timing information
285   if (fTotalTime > 0) {
286     Info("Start", "Cannot start twice");
287     return;
288   }
289   fTimer.Start(kTRUE);
290 }
291
292 //______________________________________________________________________________
293 void AliTransportMonitor::Stop()
294 {
295 // Stop the global timer
296   fTimer.Stop();
297 }
298
299 //______________________________________________________________________________
300 void AliTransportMonitor::Export(const char *fname)
301 {
302 // Export information to file
303   TFile *file = TFile::Open(fname, "RECREATE");
304   Write();
305   file->Write();
306   file->Close();
307 }  
308 //______________________________________________________________________________
309 void AliTransportMonitor::Merge(AliTransportMonitor* mergeMon)
310 {
311
312   // merge with monitor 
313   if (!fVolumeMon) 
314     {
315       TObjArray* arr = mergeMon->GetVolumes();
316       Int_t nvol = arr->GetEntriesFast();
317       fVolumeMon = new TObjArray(nvol);
318       fVolumeMon->SetOwner();
319       for (Int_t i = 0; i < nvol; i++) {
320         AliTransportMonitorVol *volMon = new AliTransportMonitorVol();
321         volMon->SetName(arr->At(i)->GetName());
322         fVolumeMon->Add(volMon);
323       }
324     } // first time
325
326
327   Int_t n = fVolumeMon->GetEntriesFast();
328   TObjArray* mergeVols = mergeMon->GetVolumes();
329   fTotalTime = 0;
330   for (Int_t i = 0; i < n; i++)
331     {
332       AliTransportMonitorVol *volMon1 = (AliTransportMonitorVol*)fVolumeMon->At(i);      
333       AliTransportMonitorVol *volMon2 = (AliTransportMonitorVol*)mergeVols->At(i);      
334       volMon1->Merge(volMon2);
335       fTotalTime += (volMon1->GetTotalTime());
336     }
337 }
338 //______________________________________________________________________________
339 AliTransportMonitor *AliTransportMonitor::Import(const char *fname)
340 {
341 // Import information from a file
342   TFile *file = TFile::Open(fname);
343   if (!file) {
344     ::Error("Import", "File %s could not be opened", fname);
345     return 0;
346   }
347   AliTransportMonitor *mon = (AliTransportMonitor *)file->Get("AliTransportMonitor");
348   if (!mon) {
349     ::Error("Import", "No AliTransportMonitor object found n file %s", fname);
350     return 0;
351   }
352   AliPDG::AddParticlesToPdgDataBase();
353   return mon;
354 }
355
356