Scripts/MaterialMapping/Mat_map_surface_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 "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, 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_X0",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("COLZ");
0040     vol->Draw();
0041     surface->Draw();
0042     surface_z->Draw();
0043     c1->Print( (out_name+"_X0.pdf").c_str());
0044     //c1->Print( (out_name+"_X0.root").c_str());
0045
0046     delete c1;
0047
0048     delete vol;
0049     delete surface;
0050     delete surface_z;
0051   }
0052
0053   // Cylinder
0054   if(surface_info.type == 1){
0055
0056     TText *vol = new TText(.1,.95,surface_info.name.c_str());
0057     vol->SetNDC();
0058     TText *surface = new TText(.1,.9,surface_info.id.c_str());
0059     surface->SetNDC();
0060     TText *surface_r = new TText(.1,.85,("R = " + to_string(surface_info.pos)).c_str() );
0061     surface_r->SetNDC();
0062     TCanvas *c1 = new TCanvas("c1","mat_X0",1200,1200);
0063     c1->SetRightMargin(0.14);
0064     c1->SetTopMargin(0.14);
0065     c1->SetLeftMargin(0.14);
0066     c1->SetBottomMargin(0.14);
0067     Map[0]->Draw("COLZ");
0068     vol->Draw();
0069     surface->Draw();
0070     surface_r->Draw();
0071     c1->Print( (out_name+"_X0.pdf").c_str());
0072     //c1->Print( (out_name+"_X0.root").c_str());
0073
0074     delete c1;
0075
0076     delete vol;
0077     delete surface;
0078     delete surface_r;
0079   }
0080   return;
0081 }
0082
0083 /// Initialise the histograms for each surface.
0084
0085 void Initialise_hist(std::vector<TH2F*>& surface_hist,
0086   const sinfo& surface_info){
0087
0088   TH2F * Map_X0;
0089   TH2F * Map_L0;
0090
0091   TH2F * Map_scale;
0092
0093   if(surface_info.type == 1){
0094     Map_X0    = new TH2F(("Map_X0_"+surface_info.idname).c_str(),("Map_X0_"+surface_info.idname).c_str(),
0095                          50,-6,6,50,-3.2,3.2);
0096     Map_L0    = new TH2F(("Map_L0_"+surface_info.idname).c_str(),("Map_L0_"+surface_info.idname).c_str(),
0097                          50,-6,6,50,-3.2,3.2);
0098     Map_scale = new TH2F(("Map_scale_"+surface_info.idname).c_str(),("Map_scale_"+surface_info.idname).c_str(),
0099                           50,-6,6,50,-3.2,3.2);
0100     Map_X0->GetXaxis()->SetTitle("Eta");
0101     Map_X0->GetYaxis()->SetTitle("Phi");
0102     Map_X0->GetZaxis()->SetTitle("X0");
0103     Map_L0->GetXaxis()->SetTitle("Eta");
0104     Map_L0->GetYaxis()->SetTitle("Phi");
0105     Map_L0->GetZaxis()->SetTitle("L0");
0106   }
0107
0108   if(surface_info.type == 2 || surface_info.type == 4){
0109     Map_X0    = new TH2F(("Map_X0_"+surface_info.idname).c_str(),("Map_X0_"+surface_info.idname).c_str(),
0110                           50,-1*surface_info.range_max,surface_info.range_max,50,-1*surface_info.range_max, surface_info.range_max);
0111     Map_L0    = new TH2F(("Map_L0_"+surface_info.idname).c_str(),("Map_L0_"+surface_info.idname).c_str(),
0112                           50,-1*surface_info.range_max,surface_info.range_max,50,-1*surface_info.range_max, surface_info.range_max);
0113     Map_scale = new TH2F(("Map_scale_"+surface_info.idname).c_str(),("Map_scale_"+surface_info.idname).c_str(),
0114                           50,-1*surface_info.range_max,surface_info.range_max,50,-1*surface_info.range_max, surface_info.range_max);
0115     Map_X0->GetXaxis()->SetTitle("X [mm]");
0116     Map_X0->GetYaxis()->SetTitle("Y [mm]");
0117     Map_X0->GetZaxis()->SetTitle("X0");
0118     Map_L0->GetXaxis()->SetTitle("X [mm]");
0119     Map_L0->GetYaxis()->SetTitle("Y [mm]");
0120     Map_L0->GetZaxis()->SetTitle("L0");
0121   }
0122   std::vector<TH2F*> v_hist;
0123   v_hist.push_back(Map_X0);
0124   v_hist.push_back(Map_L0);
0125   v_hist.push_back(Map_scale);
0126   surface_hist = v_hist;
0127 }
0128
0129 /// Fill the histograms for each surfaces.
0130
0131 void Fill(std::map<std::uint64_t,std::vector<TH2F*>>& surface_hist,  std::map<std::uint64_t,sinfo>& surface_info,
0132   const std::string& input_file, const int& nbprocess){
0133
0134   std::map<std::string,std::string> surface_name;
0135
0136   std::map<std::uint64_t,float> surface_weight;
0137
0138   //Get file, tree and set top branch address
0139   TFile *tfile = new TFile(input_file.c_str());
0140   TTree *tree = (TTree*)tfile->Get("material-tracks");
0141
0142   float v_phi   = 0;
0143   float v_eta   = 0;
0144   std::vector<float> *mat_X0   = 0;
0145   std::vector<float> *mat_L0   = 0;
0146   std::vector<float> *mat_step_length = 0;
0147
0148   std::vector<std::uint64_t> *sur_id = 0;
0149   std::vector<std::int32_t> *sur_type = 0;
0150   std::vector<float> *sur_x = 0;
0151   std::vector<float> *sur_y = 0;
0152   std::vector<float> *sur_z = 0;
0153   std::vector<float> *sur_range_min = 0;
0154   std::vector<float> *sur_range_max = 0;
0155
0156   tree->SetBranchAddress("v_phi",&v_phi);
0157   tree->SetBranchAddress("v_eta",&v_eta);
0158   tree->SetBranchAddress("mat_X0",&mat_X0);
0159   tree->SetBranchAddress("mat_L0",&mat_L0);
0160   tree->SetBranchAddress("mat_step_length",&mat_step_length);
0161
0162   tree->SetBranchAddress("sur_id",&sur_id);
0163   tree->SetBranchAddress("sur_type",&sur_type);
0164   tree->SetBranchAddress("sur_x",&sur_x);
0165   tree->SetBranchAddress("sur_y",&sur_y);
0166   tree->SetBranchAddress("sur_z",&sur_z);
0167   tree->SetBranchAddress("sur_range_min",&sur_range_min);
0168   tree->SetBranchAddress("sur_range_max",&sur_range_max);
0169
0170   int nentries = tree->GetEntries();
0171   if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
0172   // Loop over all the material tracks.
0173   for (Long64_t i=0;i<nentries; i++) {
0174     if(i%10000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
0175     tree->GetEntry(i);
0176
0177     // Reset the weight
0178     for (auto weight_it = surface_weight.begin(); weight_it != surface_weight.end(); weight_it++){
0179       weight_it->second = 0;
0180     }
0181     // loop over all the material hits to do initialisation and compute weight
0182     for(int j=0; j<mat_X0->size(); j++ ){
0183
0184       // Ignore surface of incorrect type
0185       if(sur_type->at(j) == -1) continue;
0186       // If a surface was never encountered initialise the hist, info and weight
0187       if(surface_hist.find(sur_id->at(j))==surface_hist.end()){
0188         float pos;
0189         if(sur_type->at(j) == 1){
0190           pos = sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j));
0191         }
0192         if(sur_type->at(j) == 2 || sur_type->at(j) == 4){
0193           pos = sur_z->at(j);
0194         }
0195         surface_weight[sur_id->at(j)] = 0;
0196         Initialise_info(surface_info[sur_id->at(j)], surface_name, sur_id->at(j), sur_type->at(j), pos, sur_range_min->at(j), sur_range_max->at(j));
0197         Initialise_hist(surface_hist[sur_id->at(j)], surface_info[sur_id->at(j)]);
0198       }
0199       // Weight for each surface = number of hit associated to it.
0200       surface_weight[sur_id->at(j)]++;
0201     }
0202
0203     // loop over all the material hit to fill the histogram
0204     for(int j=0; j<mat_X0->size(); j++ ){
0205
0206       // Ignore surface of incorrect type
0207       if(sur_type->at(j) == -1) continue;
0208
0209       if(sur_type->at(j) == 1){
0210         surface_hist[sur_id->at(j)][0]->Fill(v_eta, v_phi, (mat_step_length->at(j)/mat_X0->at(j)));
0211         surface_hist[sur_id->at(j)][1]->Fill(v_eta, v_phi, (mat_step_length->at(j)/mat_L0->at(j)));
0212         surface_hist[sur_id->at(j)][2]->Fill(v_eta, v_phi, (1/surface_weight[sur_id->at(j)]));
0213       }
0214       if(sur_type->at(j) == 2 || sur_type->at(j) == 4){
0215         surface_hist[sur_id->at(j)][0]->Fill(sur_x->at(j), sur_y->at(j), (mat_step_length->at(j)/mat_X0->at(j)));
0216         surface_hist[sur_id->at(j)][1]->Fill(sur_x->at(j), sur_y->at(j), (mat_step_length->at(j)/mat_L0->at(j)));
0217         surface_hist[sur_id->at(j)][2]->Fill(sur_x->at(j), sur_y->at(j), (1/surface_weight[sur_id->at(j)]));
0218       }
0219     }
0220   }
0221   // Normalise the histograms
0222   for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
0223     hist_it->second[0]->Divide(hist_it->second[2]);
0224     hist_it->second[1]->Divide(hist_it->second[2]);
0225   }
0226 }
0227
0228 /// Plot the material on each surface.
0229 /// nbprocess : number of parameter to be processed.
0230 /// name : name of the output directory.
0231
0232 void Mat_map_surface_plot(std::string input_file = "", int nbprocess = -1, std::string name = ""){
0233
0234   gStyle->SetOptStat(0);
0235   gStyle->SetOptTitle(0);
0236
0237   std::map<std::uint64_t,std::vector<TH2F*>> surface_hist;
0238   std::map<std::uint64_t,sinfo> surface_info;
0239
0240   Fill(surface_hist, surface_info, input_file, nbprocess);
0241   for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
0242     plot(hist_it->second, surface_info[hist_it->first], name);
0243     for (auto hist : hist_it->second){
0244       delete hist;
0245     }
0246     hist_it->second.clear();
0247   }
0248 }