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   std::vector<std::int32_t> *sur_type = 0;
0120
0121   std::vector<std::uint64_t> *vol_id = 0;
0122
0123   tree->SetBranchAddress("v_phi",&v_phi);
0124   tree->SetBranchAddress("v_eta",&v_eta);
0125   tree->SetBranchAddress("t_X0",&t_X0);
0126   tree->SetBranchAddress("mat_X0",&mat_X0);
0127   tree->SetBranchAddress("mat_L0",&mat_L0);
0128   tree->SetBranchAddress("mat_step_length",&mat_step_length);
0129
0130   tree->SetBranchAddress("sur_id",&sur_id);
0131   tree->SetBranchAddress("sur_type",&sur_type);
0132
0133   tree->SetBranchAddress("vol_id",&vol_id);
0134
0135   int nentries = tree->GetEntries();
0136   if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
0137   // Loop over all the material tracks.
0138   for (Long64_t i=0;i<nentries; i++) {
0139     if(i%10000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
0140     tree->GetEntry(i);
0141
0142     double matX0 = 0;
0143     double matL0 = 0;
0144
0145     // loop over all the material hits
0146     for(int j=0; j<mat_X0->size(); j++ ){
0147
0148       Acts::GeometryIdentifier ID;
0149
0150       if(sur_id->at(j) != 0){
0151         ID = Acts::GeometryIdentifier(sur_id->at(j));
0152       }
0153       else if(vol_id->at(j) != 0){
0154         ID = Acts::GeometryIdentifier(vol_id->at(j));
0155       }
0156
0157       // Check if the volume/surface is part of the selected ones
0158       if(rangeContainsValue(detectors, ID.volume())) {
0159         matX0 += mat_step_length->at(j) / mat_X0->at(j);
0160         matL0 += mat_step_length->at(j) / mat_L0->at(j);
0161       }
0162     }
0163
0164     if (matX0 != 0 && matL0 != 0){
0165       detector_hist[0]->Fill(v_eta, v_phi, matX0);
0166       detector_hist[1]->Fill(v_eta, v_phi, matL0);
0167       detector_hist[2]->Fill(v_eta, v_phi);
0168     }
0169   }
0170   detector_hist[0]->Divide(detector_hist[2]);
0171   detector_hist[1]->Divide(detector_hist[2]);
0172 }
0173
0174 /// Plot the material as function of eta and phi for a given detector/sub-detector
0175 /// detectors : list of the ID of the volume constitutive of the detector/sub-detector
0176 /// nbprocess : number of parameter to be processed.
0177 /// name : name of the output directory.
0178
0179 void Mat_map_detector_plot(std::string input_file = "",
0180                          std::vector<int> detectors = vector<int>(), int nbprocess = -1,
0181                          std::string name = "") {
0182   gStyle->SetOptStat(0);
0183   gStyle->SetOptTitle(0);
0184
0185   std::vector<TH2F*> detector_hist;
0186
0187   Fill(detector_hist, input_file, detectors, nbprocess);
0188   plot(detector_hist, detectors, name);
0189 }