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