Scripts/MaterialMapping/Mat_map_surface_plot_1D.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 two 1D histograms for each surface.
0018
0019 void plot(std::vector<TH1F*> Map, const sinfo& surface_info, const std::string& name){
0020
0021   std::string out_name = name+"/"+surface_info.name+"/"+surface_info.name+"_"+surface_info.idname;
0022   gSystem->Exec( Form("mkdir %s", (name+"/"+surface_info.name).c_str()) );
0023
0024   // Disk
0025   if(surface_info.type == 2 || surface_info.type == 4){
0026
0027     TText *vol = new TText(.1,.95,surface_info.name.c_str());
0028     vol->SetNDC();
0029     TText *surface = new TText(.1,.9,surface_info.id.c_str());
0030     surface->SetNDC();
0031     TText *surface_z = new TText(.1,.85,("Z = " + to_string(surface_info.pos)).c_str() );
0032     surface_z->SetNDC();
0033
0034     TCanvas *c1 = new TCanvas("c1","mat_R",1200,1200);
0035     c1->SetRightMargin(0.14);
0036     c1->SetTopMargin(0.14);
0037     c1->SetLeftMargin(0.14);
0038     c1->SetBottomMargin(0.14);
0039     Map[0]->Draw("P HIST");
0040     vol->Draw();
0041     surface->Draw();
0042     surface_z->Draw();
0043     c1->Print( (out_name+"_R.pdf").c_str());
0044     //c1->Print( (out_name+"_X0.root").c_str());
0045
0046     TCanvas *c2 = new TCanvas("c2","mat_Phi",1200,1200);
0047     c2->SetRightMargin(0.14);
0048     c2->SetTopMargin(0.14);
0049     c2->SetLeftMargin(0.14);
0050     c2->SetBottomMargin(0.14);
0051     Map[1]->Draw("P HIST");
0052     vol->Draw();
0053     surface->Draw();
0054     surface_z->Draw();
0055     c2->Print( (out_name+"_Phi.pdf").c_str());
0056     //c2->Print( (out_name+"_Phi.root").c_str());
0057
0058     delete c1;
0059     delete c2;
0060
0061     delete vol;
0062     delete surface;
0063     delete surface_z;
0064   }
0065
0066   // Cylinder
0067   if(surface_info.type == 1){
0068
0069     TText *vol = new TText(.1,.95,surface_info.name.c_str());
0070     vol->SetNDC();
0071     TText *surface = new TText(.1,.9,surface_info.id.c_str());
0072     surface->SetNDC();
0073     TText *surface_r = new TText(.1,.85,("R = " + to_string(surface_info.pos)).c_str() );
0074     surface_r->SetNDC();
0075     TCanvas *c1 = new TCanvas("c1","mat_Z",1200,1200);
0076     c1->SetRightMargin(0.14);
0077     c1->SetTopMargin(0.14);
0078     c1->SetLeftMargin(0.14);
0079     c1->SetBottomMargin(0.14);
0080     Map[0]->Draw("P HIST");
0081     vol->Draw();
0082     surface->Draw();
0083     surface_r->Draw();
0084     c1->Print( (out_name+"_Z.pdf").c_str());
0085     //c1->Print( (out_name+"_X0.root").c_str());
0086
0087     TCanvas *c2 = new TCanvas("c2","mat_Phi",1200,1200);
0088     c2->SetRightMargin(0.14);
0089     c2->SetTopMargin(0.14);
0090     c2->SetLeftMargin(0.14);
0091     c2->SetBottomMargin(0.14);
0092     Map[1]->Draw("P HIST");
0093     vol->Draw();
0094     surface->Draw();
0095     surface_r->Draw();
0096     c2->Print( (out_name+"_Phi.pdf").c_str());
0097     //c2->Print( (out_name+"_Phi.root").c_str());
0098
0099     delete c1;
0100     delete c2;
0101
0102     delete vol;
0103     delete surface;
0104     delete surface_r;
0105   }
0106   return;
0107 }
0108
0109 /// Initialise the two 1D histograms for each surface.
0110
0111 void Initialise_hist(std::vector<TH1F*>& surface_hist,
0112   const sinfo& surface_info){
0113
0114   TH1F * Map;
0115   TH1F * Map_Phi;
0116   TH1F * Map_scale;
0117   TH1F * Map_scale_Phi;
0118
0119   if(surface_info.type == 1){
0120     Map           = new TH1F(("Map_"+surface_info.idname).c_str(),("Map_"+surface_info.idname).c_str(),
0121                               50,-1*surface_info.range_max, surface_info.range_max);
0122     Map_Phi       = new TH1F(("Map_Phi_"+surface_info.idname).c_str(),("Map_Phi_"+surface_info.idname).c_str(),
0123                               50,-3.2,3.2);
0124     Map_scale     = new TH1F(("Map_scale_"+surface_info.idname).c_str(),("Map_scale_"+surface_info.idname).c_str(),
0125                               50,-1*surface_info.range_max, surface_info.range_max);
0126     Map_scale_Phi = new TH1F(("Map_scale_Phi_"+surface_info.idname).c_str(),("Map_scale_Phi_"+surface_info.idname).c_str(),
0127                                 50,-3.2,3.2);
0128     Map->GetXaxis()->SetTitle("Z [mm]");
0129     Map->GetYaxis()->SetTitle("X0");
0130     Map_Phi->GetXaxis()->SetTitle("Phi");
0131     Map_Phi->GetYaxis()->SetTitle("X0");
0132   }
0133
0134   if(surface_info.type == 2 || surface_info.type == 4){
0135     Map           = new TH1F(("Map_"+surface_info.idname).c_str(),("Map_"+surface_info.idname).c_str(),
0136                               50,surface_info.range_min, surface_info.range_max);
0137     Map_Phi       = new TH1F(("Map_Phi_"+surface_info.idname).c_str(),("Map_Phi_"+surface_info.idname).c_str(),
0138                               50,-3.2,3.2);
0139     Map_scale     = new TH1F(("Map_scale_"+surface_info.idname).c_str(),("Map_scale_"+surface_info.idname).c_str(),
0140                               50,surface_info.range_min, surface_info.range_max);
0141     Map_scale_Phi = new TH1F(("Map_scale_Phi_"+surface_info.idname).c_str(),("Map_scale_Phi_"+surface_info.idname).c_str(),
0142                                 50,-3.2,3.2);
0143     Map->GetXaxis()->SetTitle("R [mm]");
0144     Map->GetYaxis()->SetTitle("X0");
0145     Map_Phi->GetXaxis()->SetTitle("Phi");
0146     Map_Phi->GetYaxis()->SetTitle("X0");
0147   }
0148
0149   Map->SetMarkerStyle(20);
0150   Map_Phi->SetMarkerStyle(20);
0151
0152   std::vector<TH1F*> v_hist;
0153   v_hist.push_back(Map);
0154   v_hist.push_back(Map_Phi);
0155   v_hist.push_back(Map_scale);
0156   v_hist.push_back(Map_scale_Phi);
0157   surface_hist = v_hist;
0158 }
0159
0160   /// Fill the two 1D histograms for each surfaces.
0161
0162 void Fill(std::map<std::uint64_t,std::vector<TH1F*>>& surface_hist,  std::map<std::uint64_t,sinfo>& surface_info,
0163   const std::string& input_file, const std::string& geometry_file, const int& nbprocess){
0164   std::map<std::string,std::string> surface_name;
0165   std::map<std::uint64_t,float> surface_weight;
0166
0167   //Get file, tree and set top branch address
0168   TFile *tfile = new TFile(input_file.c_str());
0169   TTree *tree = (TTree*)tfile->Get("material_tracks");
0170
0171   float v_phi   = 0;
0172   float v_eta   = 0;
0173   std::vector<float> *mat_X0   = 0;
0174   std::vector<float> *mat_L0   = 0;
0175   std::vector<float> *mat_step_length = 0;
0176
0177   std::vector<std::uint64_t> *sur_id = 0;
0178   std::vector<float> *sur_x = 0;
0179   std::vector<float> *sur_y = 0;
0180   std::vector<float> *sur_z = 0;
0181
0182   tree->SetBranchAddress("v_phi",&v_phi);
0183   tree->SetBranchAddress("v_eta",&v_eta);
0184   tree->SetBranchAddress("mat_X0",&mat_X0);
0185   tree->SetBranchAddress("mat_L0",&mat_L0);
0186   tree->SetBranchAddress("mat_step_length",&mat_step_length);
0187
0188   tree->SetBranchAddress("sur_id",&sur_id);
0189   tree->SetBranchAddress("sur_x",&sur_x);
0190   tree->SetBranchAddress("sur_y",&sur_y);
0191   tree->SetBranchAddress("sur_z",&sur_z);
0192
0193   int nentries = tree->GetEntries();
0194   if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
0195   std::unordered_map<std::uint64_t, json> surface_bounds = load_geometry_file(geometry_file);
0196
0197   // Loop over all the material tracks.
0198   for (Long64_t i=0;i<nentries; i++) {
0199     if(i%10000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
0200     tree->GetEntry(i);
0201
0202     // Reset the weight
0203     for (auto weight_it = surface_weight.begin(); weight_it != surface_weight.end(); weight_it++){
0204       weight_it->second = 0;
0205     }
0206     // loop over all the material hits to do initialisation and compute weight
0207     for(int j=0; j<mat_X0->size(); j++ ){
0208
0209       // If a valid surface id was never encountered initialise the hist, info and weight
0210       if(surface_hist.find(sur_id->at(j))==surface_hist.end()){
0211         int type = -1;
0212         float range_min = 0.;
0213         float range_max = 0.;
0214         if (surface_bounds.count(sur_id->at(j))) {
0215           json &bounds = surface_bounds[sur_id->at(j)];
0216           std::string btype = bounds["type"].get<std::string>();
0217           const auto &values = bounds["values"];
0218           if (btype == "CylinderBounds") {
0219             type = 1;
0220             range_min = -values[1].get<float>();
0221             range_max = values[1].get<float>();
0222           } else if (btype == "RadialBounds") {
0223             type = 2;
0224             range_min = values[0].get<float>();
0225             range_max = values[1].get<float>();
0226           } else {
0227             type = -1;
0228           }
0229         }
0230
0231         float pos;
0232         if(type == 1){
0233           pos = sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j));
0234         }
0235         if(type == 2 || type == 4){
0236           pos = sur_z->at(j);
0237         }
0238
0239         // Ignore surface of incorrect type
0240         if(type == -1) continue;
0241
0242         // Weight for each surface = number of hit associated to it.
0243         surface_weight[sur_id->at(j)] = 0;
0244         Initialise_info(surface_info[sur_id->at(j)], surface_name, sur_id->at(j), type, pos, range_min, range_max);
0245         Initialise_hist(surface_hist[sur_id->at(j)], surface_info[sur_id->at(j)]);
0246       }
0247       // Weight for each surface = number of hit associated to it.
0248       surface_weight[sur_id->at(j)]++;
0249     }
0250
0251     // loop over all the material hit to fill the histogram
0252     for(int j=0; j<mat_X0->size(); j++ ){
0253
0254       // Ignore surface of incorrect type
0255       if(!surface_info.count(sur_id->at(j))) continue;
0256
0257       if(surface_info[sur_id->at(j)].type == 1){
0258         surface_hist[sur_id->at(j)][0]->Fill(sur_z->at(j), (mat_step_length->at(j)/mat_X0->at(j)));
0259         surface_hist[sur_id->at(j)][1]->Fill(v_phi, (mat_step_length->at(j)/mat_L0->at(j)));
0260         surface_hist[sur_id->at(j)][2]->Fill(sur_z->at(j), (1/surface_weight[sur_id->at(j)]));
0261         surface_hist[sur_id->at(j)][3]->Fill(v_phi, (1/surface_weight[sur_id->at(j)]));
0262       }
0263       if(surface_info[sur_id->at(j)].type == 2 || surface_info[sur_id->at(j)].type == 4){
0264         surface_hist[sur_id->at(j)][0]->Fill(sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j)), (mat_step_length->at(j)/mat_X0->at(j)));
0265         surface_hist[sur_id->at(j)][1]->Fill(v_phi, (mat_step_length->at(j)/mat_L0->at(j)));
0266         surface_hist[sur_id->at(j)][2]->Fill(sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j)), (1/surface_weight[sur_id->at(j)]));
0267         surface_hist[sur_id->at(j)][3]->Fill(v_phi, (1/surface_weight[sur_id->at(j)]));
0268       }
0269     }
0270   }
0271   // Normalise the histograms
0272   for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
0273     hist_it->second[0]->Divide(hist_it->second[2]);
0274     hist_it->second[1]->Divide(hist_it->second[3]);
0275   }
0276 }
0277
0278 /// Plot the material on each surface as function of Phi and Eta (two 1D plots).
0279 /// nbprocess : number of parameter to be processed.
0280 /// name : name of the output directory.
0281
0282 void Mat_map_surface_plot_1D(std::string input_file = "", int nbprocess = -1, std::string name = "", std::string geometry_file = "geometry-map.json"){
0283
0284   gStyle->SetOptStat(0);
0285   gStyle->SetOptTitle(0);
0286
0287   std::map<std::uint64_t,std::vector<TH1F*>> surface_hist;
0288   std::map<std::uint64_t,sinfo> surface_info;
0289
0290   Fill(surface_hist, surface_info, input_file, geometry_file, nbprocess);
0291   for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
0292     plot(hist_it->second, surface_info[hist_it->first], name);
0293     for (auto hist : hist_it->second){
0294       delete hist;
0295     }
0296     hist_it->second.clear();
0297   }
0298 }