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