Scripts/MaterialMapping/Mat_map_surface_plot_dist.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 "materialPlotHelper.cpp"
0012
0013 #include <fstream>
0014 #include <iostream>
0015 #include <sstream>
0016
0017 // Draw the plot for each surface.
0018 void plot(TGraph* Dist, const sinfo& surface_info, const std::string& name){
0019
0020   std::string out_name = name+"/"+surface_info.name+"/"+surface_info.name+"_"+surface_info.idname;
0021   gSystem->Exec( Form("mkdir %s", (name+"/"+surface_info.name).c_str()) );
0022
0023   TCanvas *c = new TCanvas("c","dist",1200,1200);
0024   c->SetRightMargin(0.14);
0025   c->SetTopMargin(0.14);
0026   c->SetLeftMargin(0.14);
0027   c->SetBottomMargin(0.14);
0028
0029   Dist->Draw("AP");
0030   TLine *line_pos;
0031
0032   TText *vol = new TText(.1,.95,surface_info.name.c_str());
0033   TText *surface = new TText(.1,.9,surface_info.id.c_str());
0034   TText *surface_z = new TText(.1,.85,("Z = " + to_string(surface_info.pos)).c_str() );
0035   TText *surface_r = new TText(.1,.85,("R = " + to_string(surface_info.pos)).c_str() );
0036   vol->SetNDC();
0037   surface->SetNDC();
0038   surface_z->SetNDC();
0039   surface_r->SetNDC();
0040
0041   // Disk
0042   if(surface_info.type == 2 || surface_info.type == 4){
0043     vol->Draw();
0044     surface->Draw();
0045     surface_z->Draw();
0046
0047     Dist->GetYaxis()->SetRangeUser(surface_info.range_min-(surface_info.range_max-surface_info.range_min)/10,
0048                                    surface_info.range_max+(surface_info.range_max-surface_info.range_min)/10);
0049
0050     // Position of the disk surface
0051     line_pos = new TLine(surface_info.pos,surface_info.range_min,surface_info.pos,surface_info.range_max);
0052     line_pos->SetLineColor(kRed);
0053     line_pos->Draw("Same");
0054   }
0055
0056   // Cylinder
0057   if(surface_info.type == 1){
0058     vol->Draw();
0059     surface->Draw();
0060     surface_r->Draw();
0061
0062     Dist->GetYaxis()->SetRangeUser(surface_info.range_min-(surface_info.range_max-surface_info.range_min)/20,
0063                                    surface_info.range_max+(surface_info.range_max-surface_info.range_min)/20);
0064
0065     // Position of the cylinder surface
0066     line_pos = new TLine(-1*surface_info.range_max, surface_info.pos, surface_info.range_max, surface_info.pos);
0067     line_pos->SetLineColor(kRed);
0068     line_pos->Draw("Same");
0069   }
0070
0071   c->Print( (out_name+"_Dist.pdf").c_str());
0072   //c->Print( (out_name+"_Dist.root").c_str());
0073
0074   delete c;
0075
0076   delete vol;
0077   delete surface;
0078   delete surface_z;
0079   delete line_pos;
0080 }
0081
0082
0083 /// Create the Tgraphs from each surface based on a vector of positions.
0084
0085 void Initialise_hist(TGraph*& surface_hist,
0086   const std::pair<std::vector<float>,std::vector<float>>& surface_pos, const sinfo& surface_info){
0087
0088   if(surface_info.type != -1){
0089     TGraph * Dist = new TGraph(surface_pos.first.size(), &surface_pos.second[0], &surface_pos.first[0]);
0090     Dist->Draw();
0091     Dist->GetXaxis()->SetTitle("Z [mm]");
0092     Dist->GetYaxis()->SetTitle("R [mm]");
0093     surface_hist = Dist;
0094   }
0095 }
0096
0097 /// Fill the histograms for each surfaces.
0098
0099 void Fill(std::map<std::uint64_t,TGraph*>& surface_hist,  std::map<std::uint64_t,sinfo>& surface_info,
0100   const std::string& input_file, const std::string& geometry_file, const int& nbprocess){
0101   std::map<std::string,std::string> surface_name;
0102   std::map<std::uint64_t,std::pair<std::vector<float>,std::vector<float>>> surface_pos;
0103
0104   //Get file, tree and set top branch address
0105   TFile *tfile = new TFile(input_file.c_str());
0106   TTree *tree = (TTree*)tfile->Get("material_tracks");
0107
0108   std::vector<float> *mat_x = 0;
0109   std::vector<float> *mat_y = 0;
0110   std::vector<float> *mat_z = 0;
0111
0112   std::vector<std::uint64_t> *sur_id = 0;
0113   std::vector<float> *sur_x = 0;
0114   std::vector<float> *sur_y = 0;
0115   std::vector<float> *sur_z = 0;
0116
0117   tree->SetBranchAddress("mat_x",&mat_x);
0118   tree->SetBranchAddress("mat_y",&mat_y);
0119   tree->SetBranchAddress("mat_z",&mat_z);
0120
0121   tree->SetBranchAddress("sur_id",&sur_id);
0122   tree->SetBranchAddress("sur_x",&sur_x);
0123   tree->SetBranchAddress("sur_y",&sur_y);
0124   tree->SetBranchAddress("sur_z",&sur_z);
0125
0126   int nentries = tree->GetEntries();
0127   if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
0128   std::unordered_map<std::uint64_t, json> surface_bounds = load_geometry_file(geometry_file);
0129
0130   // Limit the number of event processed event to 10000
0131   // more could lead to errors with the TGraphs
0132   if(nentries > 10000){
0133     nentries = 10000;
0134     std::cout << "Number of event reduced to 10000" << std::endl;
0135   }
0136   // Loop over all the material tracks.
0137   for (Long64_t i=0;i<nentries; i++) {
0138     if(i%1000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
0139     tree->GetEntry(i);
0140
0141     // loop over all the material hits.
0142     for(int j=0; j<mat_x->size(); j++ ){
0143
0144       // If a valid surface id was never encountered initialise the hist, info and weight
0145       if(surface_hist.find(sur_id->at(j))==surface_hist.end()){
0146         int type = -1;
0147         float range_min = 0.;
0148         float range_max = 0.;
0149         if (surface_bounds.count(sur_id->at(j))) {
0150           json &bounds = surface_bounds[sur_id->at(j)];
0151           std::string btype = bounds["type"].get<std::string>();
0152           const auto &values = bounds["values"];
0153           if (btype == "CylinderBounds") {
0154             type = 1;
0155             range_min = -values[1].get<float>();
0156             range_max = values[1].get<float>();
0157           } else if (btype == "RadialBounds") {
0158             type = 2;
0159             range_min = values[0].get<float>();
0160             range_max = values[1].get<float>();
0161           } else {
0162             type = -1;
0163           }
0164         }
0165
0166         float pos;
0167         float range;
0168         if(type == 1){
0169           pos = sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j));
0170         }
0171         if(type == 2 || type == 4){
0172           pos = sur_z->at(j);
0173         }
0174
0175         // Ignore surface of incorrect type
0176         if(type == -1) continue;
0177
0178         Initialise_info(surface_info[sur_id->at(j)], surface_name, sur_id->at(j), type, pos, range_min, range_max);
0179       }
0180       // Fill the vector of positions for each layer.
0181       surface_pos[sur_id->at(j)].first.push_back(sqrt(mat_y->at(j)*mat_y->at(j)+mat_x->at(j)*mat_x->at(j)));
0182       surface_pos[sur_id->at(j)].second.push_back(mat_z->at(j));
0183
0184     }
0185   }
0186   // Use the vector of positions to create the TGraphs
0187   for (auto pos_it = surface_pos.begin(); pos_it != surface_pos.end(); pos_it++){
0188     Initialise_hist(surface_hist[pos_it->first], pos_it->second, surface_info[pos_it->first]);
0189   }
0190 }
0191
0192
0193 /// Plot the position of material interaction with respect with the associated surface.
0194 /// nbprocess : number of parameter to be processed.
0195 /// name : name of the output directory.
0196
0197
0198 void Mat_map_surface_plot_dist(std::string input_file = "", int nbprocess = -1, std::string name = "", std::string geometry_file = "geometry-map.json"){
0199
0200   gStyle->SetOptStat(0);
0201   gStyle->SetOptTitle(0);
0202
0203   std::map<std::uint64_t,TGraph*> surface_hist;
0204   std::map<std::uint64_t,sinfo> surface_info;
0205   Fill(surface_hist, surface_info, input_file, geometry_file, nbprocess);
0206   for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
0207     if(hist_it->second)
0208     plot(hist_it->second, surface_info[hist_it->first], name);
0209   }
0210 }