Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Mat_map_detector_plot.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Mat_map_detector_plot.C
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <TROOT.h>
10 
11 #include "materialPlotHelper.cpp"
12 
13 #include <fstream>
14 #include <iostream>
15 #include <sstream>
16 
18 
19 void plot(std::vector<TH2F*> Map, std::vector<int> detectors, const std::string& name){
20 
21  std::string sVol = "Detector volumes :";
22  for(auto const& det: detectors) {
23  sVol += " ";
24  sVol += std::to_string(det);
25  }
26  TText *vol = new TText(.1, .95, sVol.c_str());
27  vol->SetNDC();
28 
29  TCanvas *c1 = new TCanvas("c1","mat_X0",1200,1200);
30  c1->SetRightMargin(0.14);
31  c1->SetTopMargin(0.14);
32  c1->SetLeftMargin(0.14);
33  c1->SetBottomMargin(0.14);
34  Map[0]->Draw("COLZ");
35  vol->Draw();
36  c1->Print( (name+"X0.pdf").c_str());
37 
38  TH2F *Unit_Map = (TH2F*) Map[2]->Clone();
39  Unit_Map->Divide(Map[2]);
40  TCanvas *c2 = new TCanvas("c2", "mat_X0/eta", 1200, 1200);
41  c2->SetRightMargin(0.14);
42  c2->SetTopMargin(0.14);
43  c2->SetLeftMargin(0.14);
44  c2->SetBottomMargin(0.14);
45  TH1D *Proj_eta = Map[0]->ProjectionX();
46  Proj_eta->Divide(Unit_Map->ProjectionX());
47  Proj_eta->GetYaxis()->SetTitle("X0");
48  Proj_eta->SetMarkerStyle(7);
49  Proj_eta->Draw("HIST PC");
50  c2->Print( (name+"X0_eta.pdf").c_str());
51  TCanvas *c3 = new TCanvas("c3", "mat_X0/phi", 1200, 1200);
52  c3->SetRightMargin(0.14);
53  c3->SetTopMargin(0.14);
54  c3->SetLeftMargin(0.14);
55  c3->SetBottomMargin(0.14);
56  TH1D *Proj_phi = Map[0]->ProjectionY();
57  Proj_phi->Divide(Unit_Map->ProjectionY());
58  Proj_phi->GetYaxis()->SetTitle("X0");
59  Proj_phi->SetMarkerStyle(7);
60  Proj_phi->Draw("HIST PC");
61  c3->Print( (name+"X0_phi.pdf").c_str());
62 
63  delete c1;
64  delete c2;
65  delete c3;
66  delete vol;
67  delete Unit_Map;
68  return;
69 }
70 
72 
73 void Initialise_hist(std::vector<TH2F*>& detector_hist){
74 
75  TH2F * Map_X0;
76  TH2F * Map_L0;
77  TH2F * Map_scale;
78 
79  Map_X0 = new TH2F("Map_X0_detector","Map_X0_detector",
80  100,-4,4,50,-3.2,3.2);
81  Map_L0 = new TH2F("Map_L0_detector","Map_L0_detector",
82  100,-4,4,50,-3.2,3.2);
83  Map_scale = new TH2F("Map_Scale_detector","Map_Scale_detector",
84  100,-4,4,50,-3.2,3.2);
85  Map_X0->GetXaxis()->SetTitle("Eta");
86  Map_X0->GetYaxis()->SetTitle("Phi");
87  Map_X0->GetZaxis()->SetTitle("X0");
88  Map_L0->GetXaxis()->SetTitle("Eta");
89  Map_L0->GetYaxis()->SetTitle("Phi");
90  Map_L0->GetZaxis()->SetTitle("L0");
91  std::vector<TH2F*> v_hist;
92  v_hist.push_back(Map_X0);
93  v_hist.push_back(Map_L0);
94  v_hist.push_back(Map_scale);
95  detector_hist = v_hist;
96 }
97 
99 
100 void Fill(std::vector<TH2F*>& detector_hist, const std::string& input_file, std::vector<int> detectors, const int& nbprocess){
101 
102 
103  Initialise_hist(detector_hist);
104 
105  //Get file, tree and set top branch address
106  TFile *tfile = new TFile(input_file.c_str());
107  TTree *tree = (TTree*)tfile->Get("material-tracks");
108 
109  float v_phi = 0;
110  float v_eta = 0;
111  float t_X0 = 0;
112  std::vector<float> *mat_X0 = 0;
113  std::vector<float> *mat_L0 = 0;
114  std::vector<float> *mat_step_length = 0;
115 
116  std::vector<uint64_t> *sur_id = 0;
117  std::vector<int32_t> *sur_type = 0;
118 
119  std::vector<uint64_t> *vol_id = 0;
120 
121  tree->SetBranchAddress("v_phi",&v_phi);
122  tree->SetBranchAddress("v_eta",&v_eta);
123  tree->SetBranchAddress("t_X0",&t_X0);
124  tree->SetBranchAddress("mat_X0",&mat_X0);
125  tree->SetBranchAddress("mat_L0",&mat_L0);
126  tree->SetBranchAddress("mat_step_length",&mat_step_length);
127 
128  tree->SetBranchAddress("sur_id",&sur_id);
129  tree->SetBranchAddress("sur_type",&sur_type);
130 
131  tree->SetBranchAddress("vol_id",&vol_id);
132 
133  int nentries = tree->GetEntries();
134  if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
135  // Loop over all the material tracks.
136  for (Long64_t i=0;i<nentries; i++) {
137  if(i%10000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
138  tree->GetEntry(i);
139 
140  double matX0 = 0;
141  double matL0 = 0;
142 
143  // loop over all the material hits
144  for(int j=0; j<mat_X0->size(); j++ ){
145 
146  Acts::GeometryIdentifier ID;
147 
148  if(sur_id->at(j) != 0){
149  ID = Acts::GeometryIdentifier(sur_id->at(j));
150  }
151  else if(vol_id->at(j) != 0){
152  ID = Acts::GeometryIdentifier(vol_id->at(j));
153  }
154 
155  // Check if the volume/surface is part of the selected ones
156  if(std::find(detectors.begin(), detectors.end(), ID.volume()) != detectors.end()) {
157  matX0 += mat_step_length->at(j) / mat_X0->at(j);
158  matL0 += mat_step_length->at(j) / mat_L0->at(j);
159  }
160  }
161 
162  if (matX0 != 0 && matL0 != 0){
163  detector_hist[0]->Fill(v_eta, v_phi, matX0);
164  detector_hist[1]->Fill(v_eta, v_phi, matL0);
165  detector_hist[2]->Fill(v_eta, v_phi);
166  }
167  }
168  detector_hist[0]->Divide(detector_hist[2]);
169  detector_hist[1]->Divide(detector_hist[2]);
170 }
171 
176 
177 void Mat_map_detector_plot(std::string input_file = "",
178  std::vector<int> detectors = vector<int>(), int nbprocess = -1,
179  std::string name = "") {
180  gStyle->SetOptStat(0);
181  gStyle->SetOptTitle(0);
182 
183  std::vector<TH2F*> detector_hist;
184 
185  Fill(detector_hist, input_file, detectors, nbprocess);
186  plot(detector_hist, detectors, name);
187 }