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