tutorials/PHG4DSTReader/Example.C

0001 // $Id: $
0002
0003 /*!
0004  * \file Example.C
0005  * \brief See the main function Example()
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include <TFile.h>
0012 #include <TPad.h>
0013 #include <TString.h>
0014 #include <TSystem.h>
0015 #include <TTree.h>
0016 #include <cassert>
0017 #include <cmath>
0018
0019 // it is good to load the libg4dst lib and the header files,
0020 // which allow us to use the compiled function in the output class
0021 #include <calobase/RawTower.h>
0022 #include <calobase/RawTowerv1.h>
0023 #include <g4main/PHG4HitEval.h>
0024 #include <g4main/PHG4Particlev1.h>
0025 #include <g4main/PHG4Particlev2.h>
0026 #include <g4main/PHG4VtxPointv1.h>
0027 using namespace std;
0028
0029 // allow you to access them in command line after running this macro
0030 TFile *_file0 = NULL;
0031 TTree *T = NULL;
0032
0033 //! First, checkout what is in there
0034 void CheckItOut()
0035 {
0036   // This print out the content of the first event to screen.
0037   // It contain a bunch of nodes here:
0038   //  G4Hit_*.*         - PHG4Hit for various subsystems, https://www.phenix.bnl.gov/WWW/sPHENIX/doxygen/html/d3/d9e/classPHG4Hit.html
0039   //  TOWER_*.*         - RawTower for calorimeters, https://www.phenix.bnl.gov/WWW/sPHENIX/doxygen/html/d9/dd8/classRawTower.html
0040   //  PHG4Particle.*    - Geant4 particles, https://www.phenix.bnl.gov/WWW/sPHENIX/doxygen/html/de/dc9/classPHG4Particle.html
0041   //  PHG4VtxPoint.*    - Vertex point Geant4 particles were produced, https://www.phenix.bnl.gov/WWW/sPHENIX/doxygen/html/d6/d81/classPHG4VtxPoint.html
0042   T->Show(0);
0043
0044   // Anything you see from the last print can be directly plotted using T->Draw();
0045   new TCanvas("CheckItOut", "CheckItOut");
0046   T->Draw(
0047       "Sum$(G4HIT_HCALOUT.edep) + Sum$(G4HIT_ABSORBER_HCALOUT.edep)>>hSum_CEMC(30,0,30)");
0048   TH1 *hSum_CEMC = static_cast<TH1 *>(gDirectory->GetObjectChecked("hSum_CEMC", "TH1"));
0049   hSum_CEMC->SetTitle(
0050       "Total energy deposition in Outer HCal;Energy deposition (GeV)");
0051 }
0052
0053 //! Access information for the truth Geant4 particles
0054 void AcessGeant4Particles()
0055 {
0056   // I am interested in first particle in the array, which is the generated pion.
0057   // Here our life would be easier with some new definition of variables
0058   T->SetAlias("PHG4Particle_p",
0059               "1*sqrt(PHG4Particle[].fpx**2+PHG4Particle[].fpy**2+PHG4Particle[].fpz**2)");
0060   T->SetAlias("PHG4Particle_eta",
0061               "0.5*log((PHG4Particle_p+PHG4Particle[].fpz)/(PHG4Particle_p-PHG4Particle[].fpz))");
0062   T->SetAlias("PHG4Particle_pT",
0063               "1*sqrt(PHG4Particle[].fpx**2+PHG4Particle[].fpy**2)");
0064
0065   new TCanvas("AcessGeant4Particles", "AcessGeant4Particles");
0066   T->Draw("PHG4Particle_pT>>hPHG4Particle_pT(30,0,50)",
0067           "PHG4Particle[].trkid>0"  // select only the primary particle to draw (positive track ID)
0068   );
0069   TH1 *hPHG4Particle_pT = static_cast<TH1 *>(gDirectory->GetObjectChecked("hPHG4Particle_pT", "TH1"));
0070   hPHG4Particle_pT->SetTitle("pT for the first Geant4 particle;pT (GeV/c)");
0071 }
0072
0073 //! X-Y distribution for the hits
0074 //! note, as long as we load libg4eval.so, the compiled functions can be used in T->Draw. Example, I will use
0075 //! PHG4Hit::get_avg_x(), https://www.phenix.bnl.gov/WWW/sPHENIX/doxygen/html/d3/d9e/classPHG4Hit.html#a99663034e71d3f324eb878eb0e3b36ba
0076 void WhereIsTheHits()
0077 {
0078   new TCanvas("WhereIsTheHits", "WhereIsTheHits");
0079   gPad->DrawFrame(-250, -250, 250, 250,
0080                   "Geant4 Hit distribution;X (cm);Y (cm)");
0081
0082   T->Draw("G4HIT_SVTX.get_avg_y():G4HIT_SVTX.get_avg_x()", "", "* same");
0083   T->Draw(
0084       "G4HIT_ABSORBER_CEMC.get_avg_y():G4HIT_ABSORBER_CEMC.get_avg_x()>>hG4HIT_ABSORBER_CEMC(500,-250,250,500,-250,250)",
0085       "", "scat same");
0086   T->Draw(
0087       "G4HIT_ABSORBER_HCALIN.get_avg_y():G4HIT_ABSORBER_HCALIN.get_avg_x()>>hG4HIT_ABSORBER_HCALIN(500,-250,250,500,-250,250)",
0088       "", "scat same");
0089   T->Draw(
0090       "G4HIT_ABSORBER_HCALOUT.get_avg_y():G4HIT_ABSORBER_HCALOUT.get_avg_x()>>hG4HIT_ABSORBER_HCALOUT(500,-250,250,500,-250,250)",
0091       "", "scat same");
0092 }
0093
0094 void PlotCalorimeterSamplingFraction()
0095 {
0096   // Here our life would be easier with some new definition of variables
0097   T->SetAlias("CEMC_Sample",
0098               "Sum$(G4HIT_CEMC.edep)/(Sum$(G4HIT_CEMC.edep) + Sum$(G4HIT_ABSORBER_CEMC.edep))");
0099   T->SetAlias("HCALIN_Sample",
0100               "Sum$(G4HIT_HCALIN.edep)/(Sum$(G4HIT_HCALIN.edep) + Sum$(G4HIT_ABSORBER_HCALIN.edep))");
0101   T->SetAlias("HCALOUT_Sample",
0102               "Sum$(G4HIT_HCALOUT.edep)/(Sum$(G4HIT_HCALOUT.edep) + Sum$(G4HIT_ABSORBER_HCALOUT.edep))");
0103
0104   new TCanvas("PlotCalorimeterSamplingFraction", "PlotCalorimeterSamplingFraction");
0105
0106   T->Draw("CEMC_Sample>>CEMC_Sample(30,0,.15)");
0107   TH1 *CEMC_Sample = static_cast<TH1 *>(gDirectory->GetObjectChecked("CEMC_Sample", "TH1"));
0108   CEMC_Sample->SetLineColor(kRed);
0109   T->Draw("HCALIN_Sample>>HCALIN_Sample(30,0,.15)", "", "same");
0110   TH1 *HCALIN_Sample = static_cast<TH1 *>(gDirectory->GetObjectChecked("HCALIN_Sample", "TH1"));
0111   HCALIN_Sample->SetLineColor(kMagenta);
0112   T->Draw("HCALOUT_Sample>>HCALOUT_Sample(30,0,.15)", "", "same");
0113   TH1 *HCALOUT_Sample = static_cast<TH1 *>(gDirectory->GetObjectChecked("HCALOUT_Sample", "TH1"));
0114   HCALOUT_Sample->SetLineColor(kBlue);
0115   CEMC_Sample->SetTitle("Sampling fraction for three calorimeters; E_{Scint}/(E_{Scint} + E_{Absorber})");
0116 }
0117
0118 void AccessCalorimeterTowers()
0119 {
0120   new TCanvas("AccessCalorimeterTowers", "AccessCalorimeterTowers");
0121
0122   T->Draw("TOWER_CALIB_HCALOUT.get_bineta()>>hTowerBin(22,-0.5,21.5)", "TOWER_CALIB_HCALOUT.get_energy()");
0123   TH1 *hTowerBin = static_cast<TH1 *>(gDirectory->GetObjectChecked("hTowerBin", "TH1"));
0124   hTowerBin->SetTitle("HCal tower eta bin distribution; Eta Bin # ; Sum energy in scintillator (GeV)");
0125 }
0126
0127 /*
0128  * \brief Example macro to use output for PHG4DSTReader
0129  *
0130  * (s)PHENIX DST files are designed for high speed IO for batch computing,
0131  * the content is not necessarily easily readable through ROOT commandline.
0132  * This maybe some headache during testing/development stage.
0133  *
0134  * PHG4DSTReader is designed to solve this problem, by exporting a copy of
0135  * DST tree into a ROOT readable TTree file. User can directly inspect the content
0136  * using ROOT commandlines or simple macros.
0137  * https://www.phenix.bnl.gov/WWW/sPHENIX/doxygen/html/d4/dc9/classPHG4DSTReader.html
0138  *
0139  * PHG4DSTReader is by-default disabled in sPHENIX simulation, but can be enabled
0140  * for small-scale test simulation runs by following switch in  Fun4All_G4_sPHENIX.C
0141  * bool do_DSTReader = true;
0142  *
0143  * This example macro work on its output, usually named *_DSTReader.root.
0144  * A test file is prepared at https://www.phenix.bnl.gov/phenix/WWW/sPHENIX/tutorial/G4sPHENIX.root_DSTReader.root
0145  * which is output of the default Fun4All_G4_sPHENIX.C macro with 100 events and do_DSTReader = true;
0146  *
0147  * */
0148 void Example(  //
0149     const TString infile =
0150         "https://www.phenix.bnl.gov/phenix/WWW/sPHENIX/tutorial/G4sPHENIX.root_DSTReader.root" // this is output of the default Fun4All_G4_sPHENIX.C macro with 100 events and do_DSTReader = true;
0151         )
0152 {
0153   gSystem->Load("libg4dst.so");
0154
0155   // open file, get the tree
0156   _file0 = TFile::Open(infile);
0157   assert(_file0);
0158   T = (TTree *) _file0->GetObjectChecked("T", "TTree");
0159   assert(T);
0160
0161   // Now I have bunch of example to show how to use this file.
0162   // They are self-explanatory
0163   //  and all the lines can be copy-paste in ROOT command lines too
0164   CheckItOut();
0165   AcessGeant4Particles();
0166   WhereIsTheHits();
0167   PlotCalorimeterSamplingFraction();
0168   AccessCalorimeterTowers();
0169 }