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 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<std::int32_t> *sur_type = 0;
0179   std::vector<float> *sur_x = 0;
0180   std::vector<float> *sur_y = 0;
0181   std::vector<float> *sur_z = 0;
0182   std::vector<float> *sur_range_min = 0;
0183   std::vector<float> *sur_range_max = 0;
0184
0185   tree->SetBranchAddress("v_phi",&v_phi);
0186   tree->SetBranchAddress("v_eta",&v_eta);
0187   tree->SetBranchAddress("mat_X0",&mat_X0);
0188   tree->SetBranchAddress("mat_L0",&mat_L0);
0189   tree->SetBranchAddress("mat_step_length",&mat_step_length);
0190
0191   tree->SetBranchAddress("sur_id",&sur_id);
0192   tree->SetBranchAddress("sur_type",&sur_type);
0193   tree->SetBranchAddress("sur_x",&sur_x);
0194   tree->SetBranchAddress("sur_y",&sur_y);
0195   tree->SetBranchAddress("sur_z",&sur_z);
0196   tree->SetBranchAddress("sur_range_min",&sur_range_min);
0197   tree->SetBranchAddress("sur_range_max",&sur_range_max);
0198
0199   int nentries = tree->GetEntries();
0200   if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
0201   // Loop over all the material tracks.
0202   for (Long64_t i=0;i<nentries; i++) {
0203     if(i%10000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
0204     tree->GetEntry(i);
0205
0206     // Reset the weight
0207     for (auto weight_it = surface_weight.begin(); weight_it != surface_weight.end(); weight_it++){
0208       weight_it->second = 0;
0209     }
0210     // loop over all the material hits to do initialisation and compute weight
0211     for(int j=0; j<mat_X0->size(); j++ ){
0212
0213       // Ignore surface of incorrect type
0214       if(sur_type->at(j) == -1) continue;
0215       // If a surface was never encountered initialise the hist, info and weight
0216       if(surface_hist.find(sur_id->at(j))==surface_hist.end()){
0217
0218         float pos;
0219         if(sur_type->at(j) == 1){
0220           pos = sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j));
0221         }
0222         if(sur_type->at(j) == 2 || sur_type->at(j) == 4){
0223           pos = sur_z->at(j);
0224         }
0225         // Weight for each surface = number of hit associated to it.
0226         surface_weight[sur_id->at(j)] = 0;
0227         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));
0228         Initialise_hist(surface_hist[sur_id->at(j)], surface_info[sur_id->at(j)]);
0229       }
0230       // Weight for each surface = number of hit associated to it.
0231       surface_weight[sur_id->at(j)]++;
0232     }
0233
0234     // loop over all the material hit to fill the histogram
0235     for(int j=0; j<mat_X0->size(); j++ ){
0236
0237       // Ignore surface of incorrect type
0238       if(sur_type->at(j) == -1) continue;
0239
0240       if(sur_type->at(j) == 1){
0241         surface_hist[sur_id->at(j)][0]->Fill(sur_z->at(j), (mat_step_length->at(j)/mat_X0->at(j)));
0242         surface_hist[sur_id->at(j)][1]->Fill(v_phi, (mat_step_length->at(j)/mat_L0->at(j)));
0243         surface_hist[sur_id->at(j)][2]->Fill(sur_z->at(j), (1/surface_weight[sur_id->at(j)]));
0244         surface_hist[sur_id->at(j)][3]->Fill(v_phi, (1/surface_weight[sur_id->at(j)]));
0245       }
0246       if(sur_type->at(j) == 2 || sur_type->at(j) == 4){
0247         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)));
0248         surface_hist[sur_id->at(j)][1]->Fill(v_phi, (mat_step_length->at(j)/mat_L0->at(j)));
0249         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)]));
0250         surface_hist[sur_id->at(j)][3]->Fill(v_phi, (1/surface_weight[sur_id->at(j)]));
0251       }
0252     }
0253   }
0254   // Normalise the histograms
0255   for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
0256     hist_it->second[0]->Divide(hist_it->second[2]);
0257     hist_it->second[1]->Divide(hist_it->second[3]);
0258   }
0259 }
0260
0261 /// Plot the material on each surface as function of Phi and Eta (two 1D plots).
0262 /// nbprocess : number of parameter to be processed.
0263 /// name : name of the output directory.
0264
0265 void Mat_map_surface_plot_1D(std::string input_file = "", int nbprocess = -1, std::string name = ""){
0266
0267   gStyle->SetOptStat(0);
0268   gStyle->SetOptTitle(0);
0269
0270   std::map<std::uint64_t,std::vector<TH1F*>> surface_hist;
0271   std::map<std::uint64_t,sinfo> surface_info;
0272
0273   Fill(surface_hist, surface_info, input_file, nbprocess);
0274   for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
0275     plot(hist_it->second, surface_info[hist_it->first], name);
0276     for (auto hist : hist_it->second){
0277       delete hist;
0278     }
0279     hist_it->second.clear();
0280   }
0281 }