Scripts/MaterialMapping/Mat_map_detector_plot.C

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include <TROOT.h>
0010
0011 #include "materialPlotHelper.cpp"
0012
0013 #include <fstream>
0014 #include <iostream>
0015 #include <sstream>
0016
0017 /// Draw and save the histograms.
0018
0019 void plot(std::vector<TH2F*> Map, std::vector<int> detectors, const std::string& name){
0020
0021   std::string sVol = "Detector volumes :";
0022   for(auto const& det: detectors) {
0023     sVol += " ";
0024     sVol += std::to_string(det);
0025   }
0026   TText *vol = new TText(.1, .95, sVol.c_str());
0027   vol->SetNDC();
0028
0029   TCanvas *c1 = new TCanvas("c1","mat_X0",1200,1200);
0030   c1->SetRightMargin(0.14);
0031   c1->SetTopMargin(0.14);
0032   c1->SetLeftMargin(0.14);
0033   c1->SetBottomMargin(0.14);
0034   Map[0]->Draw("COLZ");
0035   vol->Draw();
0036   c1->Print( (name+"X0.pdf").c_str());
0037
0038   TH2F *Unit_Map = (TH2F*) Map[2]->Clone();
0039   Unit_Map->Divide(Map[2]);
0040   TCanvas *c2 = new TCanvas("c2", "mat_X0/eta", 1200, 1200);
0041   c2->SetRightMargin(0.14);
0042   c2->SetTopMargin(0.14);
0043   c2->SetLeftMargin(0.14);
0044   c2->SetBottomMargin(0.14);
0045   TH1D *Proj_eta = Map[0]->ProjectionX();
0046   Proj_eta->Divide(Unit_Map->ProjectionX());
0047   Proj_eta->GetYaxis()->SetTitle("X0");
0048   Proj_eta->SetMarkerStyle(7);
0049   Proj_eta->Draw("HIST PC");
0050   c2->Print( (name+"X0_eta.pdf").c_str());
0051   TCanvas *c3 = new TCanvas("c3", "mat_X0/phi", 1200, 1200);
0052   c3->SetRightMargin(0.14);
0053   c3->SetTopMargin(0.14);
0054   c3->SetLeftMargin(0.14);
0055   c3->SetBottomMargin(0.14);
0056   TH1D *Proj_phi = Map[0]->ProjectionY();
0057   Proj_phi->Divide(Unit_Map->ProjectionY());
0058   Proj_phi->GetYaxis()->SetTitle("X0");
0059   Proj_phi->SetMarkerStyle(7);
0060   Proj_phi->Draw("HIST PC");
0061   c3->Print( (name+"X0_phi.pdf").c_str());
0062
0063   delete c1;
0064   delete c2;
0065   delete c3;
0066   delete vol;
0067   delete Unit_Map;
0068   return;
0069 }
0070
0071 /// Initialise the histograms for the detector.
0072
0073 void Initialise_hist(std::vector<TH2F*>& detector_hist){
0074
0075   TH2F * Map_X0;
0076   TH2F * Map_L0;
0077   TH2F * Map_scale;
0078
0079   Map_X0    = new TH2F("Map_X0_detector","Map_X0_detector",
0080                        100,-4,4,50,-3.2,3.2);
0081   Map_L0    = new TH2F("Map_L0_detector","Map_L0_detector",
0082                        100,-4,4,50,-3.2,3.2);
0083   Map_scale = new TH2F("Map_Scale_detector","Map_Scale_detector",
0084                        100,-4,4,50,-3.2,3.2);
0085   Map_X0->GetXaxis()->SetTitle("Eta");
0086   Map_X0->GetYaxis()->SetTitle("Phi");
0087   Map_X0->GetZaxis()->SetTitle("X0");
0088   Map_L0->GetXaxis()->SetTitle("Eta");
0089   Map_L0->GetYaxis()->SetTitle("Phi");
0090   Map_L0->GetZaxis()->SetTitle("L0");
0091   std::vector<TH2F*> v_hist;
0092   v_hist.push_back(Map_X0);
0093   v_hist.push_back(Map_L0);
0094   v_hist.push_back(Map_scale);
0095   detector_hist = v_hist;
0096 }
0097
0098 /// Fill the histograms for the detector.
0099
0100 void Fill(std::vector<TH2F*>& detector_hist, const std::string& input_file, std::vector<int> detectors,  const int& nbprocess){
0101
0102
0103   Initialise_hist(detector_hist);
0104
0105   //Get file, tree and set top branch address
0106   TFile *tfile = new TFile(input_file.c_str());
0107   TTree *tree = (TTree*)tfile->Get("material-tracks");
0108
0109   float v_phi  = 0;
0110   float v_eta  = 0;
0111   float t_X0   = 0;
0112   std::vector<float> *mat_X0   = 0;
0113   std::vector<float> *mat_L0   = 0;
0114   std::vector<float> *mat_step_length = 0;
0115
0116   std::vector<uint64_t> *sur_id = 0;
0117   std::vector<int32_t> *sur_type = 0;
0118
0119   std::vector<uint64_t> *vol_id = 0;
0120
0121   tree->SetBranchAddress("v_phi",&v_phi);
0122   tree->SetBranchAddress("v_eta",&v_eta);
0123   tree->SetBranchAddress("t_X0",&t_X0);
0124   tree->SetBranchAddress("mat_X0",&mat_X0);
0125   tree->SetBranchAddress("mat_L0",&mat_L0);
0126   tree->SetBranchAddress("mat_step_length",&mat_step_length);
0127
0128   tree->SetBranchAddress("sur_id",&sur_id);
0129   tree->SetBranchAddress("sur_type",&sur_type);
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(std::find(detectors.begin(), detectors.end(), ID.volume()) != detectors.end()) {
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 detector/sub-detector
0173 /// detectors : list of the ID of the volume constitutive of the detector/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 }