Scripts/MaterialMapping/Mat_map_detector_plot.C

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <TROOT.h>
0010
0011 #include "Acts/Utilities/Helpers.hpp"
0012
0013 #include "materialPlotHelper.cpp"
0014
0015 #include <fstream>
0016 #include <iostream>
0017 #include <sstream>
0018
0019 /// Draw and save the histograms.
0020
0021 void plot(std::vector<TH2F*> Map, std::vector<int> detectors, const std::string& name){
0022
0023   std::string sVol = "Detector volumes :";
0024   for(auto const& det: detectors) {
0025     sVol += " ";
0026     sVol += std::to_string(det);
0027   }
0028   TText *vol = new TText(.1, .95, sVol.c_str());
0029   vol->SetNDC();
0030
0031   TCanvas *c1 = new TCanvas("c1","mat_X0",1200,1200);
0032   c1->SetRightMargin(0.14);
0033   c1->SetTopMargin(0.14);
0034   c1->SetLeftMargin(0.14);
0035   c1->SetBottomMargin(0.14);
0036   Map[0]->Draw("COLZ");
0037   vol->Draw();
0038   c1->Print( (name+"X0.pdf").c_str());
0039
0040   TH2F *Unit_Map = (TH2F*) Map[2]->Clone();
0041   Unit_Map->Divide(Map[2]);
0042   TCanvas *c2 = new TCanvas("c2", "mat_X0/eta", 1200, 1200);
0043   c2->SetRightMargin(0.14);
0044   c2->SetTopMargin(0.14);
0045   c2->SetLeftMargin(0.14);
0046   c2->SetBottomMargin(0.14);
0047   TH1D *Proj_eta = Map[0]->ProjectionX();
0048   Proj_eta->Divide(Unit_Map->ProjectionX());
0049   Proj_eta->GetYaxis()->SetTitle("X0");
0050   Proj_eta->SetMarkerStyle(7);
0051   Proj_eta->Draw("HIST PC");
0052   c2->Print( (name+"X0_eta.pdf").c_str());
0053   TCanvas *c3 = new TCanvas("c3", "mat_X0/phi", 1200, 1200);
0054   c3->SetRightMargin(0.14);
0055   c3->SetTopMargin(0.14);
0056   c3->SetLeftMargin(0.14);
0057   c3->SetBottomMargin(0.14);
0058   TH1D *Proj_phi = Map[0]->ProjectionY();
0059   Proj_phi->Divide(Unit_Map->ProjectionY());
0060   Proj_phi->GetYaxis()->SetTitle("X0");
0061   Proj_phi->SetMarkerStyle(7);
0062   Proj_phi->Draw("HIST PC");
0063   c3->Print( (name+"X0_phi.pdf").c_str());
0064
0065   delete c1;
0066   delete c2;
0067   delete c3;
0068   delete vol;
0069   delete Unit_Map;
0070   return;
0071 }
0072
0073 /// Initialise the histograms for the detector.
0074
0075 void Initialise_hist(std::vector<TH2F*>& detector_hist){
0076
0077   TH2F * Map_X0;
0078   TH2F * Map_L0;
0079   TH2F * Map_scale;
0080
0081   Map_X0    = new TH2F("Map_X0_detector","Map_X0_detector",
0082                        100,-4,4,50,-3.2,3.2);
0083   Map_L0    = new TH2F("Map_L0_detector","Map_L0_detector",
0084                        100,-4,4,50,-3.2,3.2);
0085   Map_scale = new TH2F("Map_Scale_detector","Map_Scale_detector",
0086                        100,-4,4,50,-3.2,3.2);
0087   Map_X0->GetXaxis()->SetTitle("Eta");
0088   Map_X0->GetYaxis()->SetTitle("Phi");
0089   Map_X0->GetZaxis()->SetTitle("X0");
0090   Map_L0->GetXaxis()->SetTitle("Eta");
0091   Map_L0->GetYaxis()->SetTitle("Phi");
0092   Map_L0->GetZaxis()->SetTitle("L0");
0093   std::vector<TH2F*> v_hist;
0094   v_hist.push_back(Map_X0);
0095   v_hist.push_back(Map_L0);
0096   v_hist.push_back(Map_scale);
0097   detector_hist = v_hist;
0098 }
0099
0100 /// Fill the histograms for the detector.
0101
0102 void Fill(std::vector<TH2F*>& detector_hist, const std::string& input_file, std::vector<int> detectors,  const int& nbprocess){
0103
0104
0105   Initialise_hist(detector_hist);
0106
0107   //Get file, tree and set top branch address
0108   TFile *tfile = new TFile(input_file.c_str());
0109   TTree *tree = (TTree*)tfile->Get("material_tracks");
0110
0111   float v_phi  = 0;
0112   float v_eta  = 0;
0113   float t_X0   = 0;
0114   std::vector<float> *mat_X0   = 0;
0115   std::vector<float> *mat_L0   = 0;
0116   std::vector<float> *mat_step_length = 0;
0117
0118   std::vector<std::uint64_t> *sur_id = 0;
0119
0120   std::vector<std::uint64_t> *vol_id = 0;
0121
0122   tree->SetBranchAddress("v_phi",&v_phi);
0123   tree->SetBranchAddress("v_eta",&v_eta);
0124   tree->SetBranchAddress("t_X0",&t_X0);
0125   tree->SetBranchAddress("mat_X0",&mat_X0);
0126   tree->SetBranchAddress("mat_L0",&mat_L0);
0127   tree->SetBranchAddress("mat_step_length",&mat_step_length);
0128
0129   tree->SetBranchAddress("sur_id",&sur_id);
0130
0131   tree->SetBranchAddress("vol_id",&vol_id);
0132
0133   int nentries = tree->GetEntries();
0134   if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
0135   // Loop over all the material tracks.
0136   for (Long64_t i=0;i<nentries; i++) {
0137     if(i%10000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
0138     tree->GetEntry(i);
0139
0140     double matX0 = 0;
0141     double matL0 = 0;
0142
0143     // loop over all the material hits
0144     for(int j=0; j<mat_X0->size(); j++ ){
0145
0146       Acts::GeometryIdentifier ID;
0147
0148       if(sur_id->at(j) != 0){
0149         ID = Acts::GeometryIdentifier(sur_id->at(j));
0150       }
0151       else if(vol_id->at(j) != 0){
0152         ID = Acts::GeometryIdentifier(vol_id->at(j));
0153       }
0154
0155       // Check if the volume/surface is part of the selected ones
0156       if(rangeContainsValue(detectors, ID.volume())) {
0157         matX0 += mat_step_length->at(j) / mat_X0->at(j);
0158         matL0 += mat_step_length->at(j) / mat_L0->at(j);
0159       }
0160     }
0161
0162     if (matX0 != 0 && matL0 != 0){
0163       detector_hist[0]->Fill(v_eta, v_phi, matX0);
0164       detector_hist[1]->Fill(v_eta, v_phi, matL0);
0165       detector_hist[2]->Fill(v_eta, v_phi);
0166     }
0167   }
0168   detector_hist[0]->Divide(detector_hist[2]);
0169   detector_hist[1]->Divide(detector_hist[2]);
0170 }
0171
0172 /// Plot the material as function of eta and phi for a given Geometry/sub-detector
0173 /// detectors : list of the ID of the volume constitutive of the Geometry/sub-detector
0174 /// nbprocess : number of parameter to be processed.
0175 /// name : name of the output directory.
0176
0177 void Mat_map_detector_plot(std::string input_file = "",
0178                          std::vector<int> detectors = vector<int>(), int nbprocess = -1,
0179                          std::string name = "") {
0180   gStyle->SetOptStat(0);
0181   gStyle->SetOptTitle(0);
0182
0183   std::vector<TH2F*> detector_hist;
0184
0185   Fill(detector_hist, input_file, detectors, nbprocess);
0186   plot(detector_hist, detectors, name);
0187 }