dna/phasespace/plot.C

0001 // *********************************************************************
0002 // To execute this macro under ROOT after your simulation ended,
0003 //   1 - launch ROOT (usually type 'root' at your machine's prompt)
0004 //   2 - type '.X plot.C' at the ROOT session prompt
0005 // *********************************************************************
0006
0007 void SetLeafAddress(TNtuple* ntuple, const char* name, void* address);
0008
0009 void plot()
0010 {
0011   gROOT->Reset();
0012   gStyle->SetPalette(1);
0013   gROOT->SetStyle("Plain");
0014
0015   TCanvas* c1 = new TCanvas ("c1","",20,20,1500,500);
0016   c1->Divide(3,1);
0017
0018   // Uncomment if merging should be done
0019   //system ("rm -rf psp.root");
0020   //system ("hadd psp.root psp_*.root");
0021
0022   TFile* f = new TFile("psp.root");
0023
0024   TNtuple* ntuple;
0025   ntuple = (TNtuple*)f->Get("scorer");
0026   bool rowWise = true;
0027   TBranch* eventBranch = ntuple->FindBranch("row_wise_branch");
0028   if ( ! eventBranch ) rowWise = false;
0029   // std::cout <<  "rowWise: " << rowWise << std::endl;
0030
0031   //*********************************************************************
0032   // canvas tab 1
0033   //*********************************************************************
0034
0035   c1->cd(1);
0036   gStyle->SetOptStat(10);
0037   ntuple->SetFillStyle(1001);
0038   ntuple->SetFillColor(2);
0039   ntuple->Draw("PDGCode","","B");
0040   gPad->SetLogy();
0041
0042   //*********************************************************************
0043   // canvas tab 2
0044   //*********************************************************************
0045
0046   c1->cd(2);
0047   ntuple->SetMarkerColor(4);
0048   ntuple->SetMarkerSize(5);
0049   ntuple->Draw("x:y:z");
0050
0051   //*********************************************************************
0052   // canvas tab 3
0053   //*********************************************************************
0054
0055   c1->cd(3);
0056   ntuple->SetMarkerColor(4);
0057   ntuple->Draw("kineticEnergy");
0058   gPad->SetLogy();
0059
0060   //*********************************************************************
0061   // CREATION OF CSV FILE IN GRAS STYLE
0062   //*********************************************************************
0063
0064   Double_t PDGCode;
0065   Double_t x;
0066   Double_t y;
0067   Double_t z;
0068   Double_t xMom;
0069   Double_t yMom;
0070   Double_t zMom;
0071   Double_t kineticEnergy;
0072
0073   if ( ! rowWise ) {
0074     ntuple->SetBranchAddress("PDGCode",&PDGCode);
0075     ntuple->SetBranchAddress("x",&x);
0076     ntuple->SetBranchAddress("y",&y);
0077     ntuple->SetBranchAddress("z",&z);
0078     ntuple->SetBranchAddress("xMom",&xMom);
0079     ntuple->SetBranchAddress("yMom",&yMom);
0080     ntuple->SetBranchAddress("zMom",&zMom);
0081     ntuple->SetBranchAddress("kineticEnergy",&kineticEnergy);
0082   }
0083   else {
0084     SetLeafAddress(ntuple, "PDGCode",&PDGCode);
0085     SetLeafAddress(ntuple, "x",&x);
0086     SetLeafAddress(ntuple, "y",&y);
0087     SetLeafAddress(ntuple, "z",&z);
0088     SetLeafAddress(ntuple, "xMom",&xMom);
0089     SetLeafAddress(ntuple, "yMom",&yMom);
0090     SetLeafAddress(ntuple, "zMom",&zMom);
0091     SetLeafAddress(ntuple, "kineticEnergy",&kineticEnergy);
0092   }
0093
0094   system ("rm -rf GRAS_tmp.csv");
0095   std::ofstream csv("GRAS_tmp.csv");
0096
0097   // REMINDER: GRAS CSV FORMAT
0098   /*
0099   'eventID','',    1,'Event ID'
0100   'trackID','',    1,'Track ID'
0101   'PDGEncoding','',    1,'PDG Encoding of particle type'
0102   'Z','',    1,'Atomic number of nucleus'
0103   'A','',    1,'Nucleon number of nucleus'
0104   'Q','',    1,'Charge'
0105   'excitation','MeV',    1,'Nuclear excitation energy'
0106   'weight','',    1,'Particle weight'
0107   'KE','MeV',    1,'Kinetic energy'
0108   'T','s',    1,'Time coordinate'
0109   'XG','mm',    3,'Position coordinate (global)'
0110   'XL','mm',    3,'Position coordinate (local)'
0111   'VG','',    3,'Direction cosine (global)'
0112   'VL','',    3,'Direction cosine (local)'
0113   'D','mm',    1,'Distance from last volume/scatter'
0114   'instanceID','',    1,'Instance ID of PV'
0115   */
0116
0117   cout << "-> Number of events = " << ntuple->GetEntries() << endl;
0118
0119   for (Int_t j=0;j<ntuple->GetEntries(); j++)
0120   {
0121     ntuple->GetEntry(j);
0122
0123     csv
0124     << "       "
0125     << j+1 << "," << "             "
0126     << 1 << "," << "             "
0127     << PDGCode << "," << "             "
0128     << 1 << "," << "             "
0129     << 1 << "," << "             "
0130     << -999 << "," << "             "
0131     << -999 << "," << "             "
0132     << -999 << "," << "             "
0133     << kineticEnergy*1e-3 << ","  << "             "// keV to MeV
0134     << 0 << "," << "             "
0135     << x*1e-3 << ","  << "             "// um to mm
0136     << y*1e-3 << ","  << "             "// um to mm
0137     << z*1e-3 << ","  << "             "// um to mm
0138     << -1 << "," << "             "
0139     << -1 << "," << "             "
0140     << -1 << "," << "             "
0141     << xMom << "," << "             "
0142     << yMom << "," << "             "
0143     << zMom << "," << "             "
0144     << -1 << "," << "             "
0145     << -1 << "," << "             "
0146     << -1 << "," << "             "
0147     << 0 << "," << "             "
0148     << 1 << endl;
0149   }
0150
0151   system ("rm -rf GRAS.csv");
0152   system ("cat GRAS_header.txt GRAS_tmp.csv GRAS_eof.txt >> GRAS.csv");
0153 }
0154
0155 void SetLeafAddress(TNtuple* ntuple, const char* name, void* address) {
0156   TLeaf* leaf = ntuple->FindLeaf(name);
0157   if ( ! leaf ) {
0158     std::cerr << "Error in <SetLeafAddress>: unknown leaf --> " << name << std::endl;
0159     return;
0160   }
0161   leaf->SetAddress(address);
0162 }