Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-18 08:23:26 

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 #include <unordered_map>
0017 
0018 #include <nlohmann/json.hpp>
0019 using json = nlohmann::json;
0020 
0021 /// Draw and save the histograms.
0022 
0023 void plot(std::vector<TH2F*> Map, const sinfo& surface_info, const std::string& name){
0024 
0025   std::string out_name = name+"/"+surface_info.name+"/"+surface_info.name+"_"+surface_info.idname;
0026   gSystem->Exec( Form("mkdir %s", (name+"/"+surface_info.name).c_str()) );
0027 
0028   // Disk
0029   if(surface_info.type == 2  || surface_info.type == 4){
0030 
0031     TText *vol = new TText(.1,.95,surface_info.name.c_str());
0032     vol->SetNDC();
0033     TText *surface = new TText(.1,.9,surface_info.id.c_str());
0034     surface->SetNDC();
0035     TText *surface_z = new TText(.1,.85,("Z = " + to_string(surface_info.pos)).c_str() );
0036     surface_z->SetNDC();
0037 
0038     TCanvas *c1 = new TCanvas("c1","mat_X0",1200,1200);
0039     c1->SetRightMargin(0.14);
0040     c1->SetTopMargin(0.14);
0041     c1->SetLeftMargin(0.14);
0042     c1->SetBottomMargin(0.14);
0043     Map[0]->Draw("COLZ");
0044     vol->Draw();
0045     surface->Draw();
0046     surface_z->Draw();
0047     c1->Print( (out_name+"_X0.pdf").c_str());
0048     //c1->Print( (out_name+"_X0.root").c_str());
0049 
0050     delete c1;
0051 
0052     delete vol;
0053     delete surface;
0054     delete surface_z;
0055   }
0056 
0057   // Cylinder
0058   if(surface_info.type == 1){
0059 
0060     TText *vol = new TText(.1,.95,surface_info.name.c_str());
0061     vol->SetNDC();
0062     TText *surface = new TText(.1,.9,surface_info.id.c_str());
0063     surface->SetNDC();
0064     TText *surface_r = new TText(.1,.85,("R = " + to_string(surface_info.pos)).c_str() );
0065     surface_r->SetNDC();
0066     TCanvas *c1 = new TCanvas("c1","mat_X0",1200,1200);
0067     c1->SetRightMargin(0.14);
0068     c1->SetTopMargin(0.14);
0069     c1->SetLeftMargin(0.14);
0070     c1->SetBottomMargin(0.14);
0071     Map[0]->Draw("COLZ");
0072     vol->Draw();
0073     surface->Draw();
0074     surface_r->Draw();
0075     c1->Print( (out_name+"_X0.pdf").c_str());
0076     //c1->Print( (out_name+"_X0.root").c_str());
0077 
0078     delete c1;
0079 
0080     delete vol;
0081     delete surface;
0082     delete surface_r;
0083   }
0084   return;
0085 }
0086 
0087 /// Initialise the histograms for each surface.
0088 
0089 void Initialise_hist(std::vector<TH2F*>& surface_hist,
0090   const sinfo& surface_info){
0091 
0092   TH2F * Map_X0;
0093   TH2F * Map_L0;
0094 
0095   TH2F * Map_scale;
0096 
0097   if(surface_info.type == 1){
0098     Map_X0    = new TH2F(("Map_X0_"+surface_info.idname).c_str(),("Map_X0_"+surface_info.idname).c_str(),
0099                          50,-6,6,50,-3.2,3.2);
0100     Map_L0    = new TH2F(("Map_L0_"+surface_info.idname).c_str(),("Map_L0_"+surface_info.idname).c_str(),
0101                          50,-6,6,50,-3.2,3.2);
0102     Map_scale = new TH2F(("Map_scale_"+surface_info.idname).c_str(),("Map_scale_"+surface_info.idname).c_str(),
0103                           50,-6,6,50,-3.2,3.2);
0104     Map_X0->GetXaxis()->SetTitle("Eta");
0105     Map_X0->GetYaxis()->SetTitle("Phi");
0106     Map_X0->GetZaxis()->SetTitle("X0");
0107     Map_L0->GetXaxis()->SetTitle("Eta");
0108     Map_L0->GetYaxis()->SetTitle("Phi");
0109     Map_L0->GetZaxis()->SetTitle("L0");
0110   }
0111 
0112   if(surface_info.type == 2 || surface_info.type == 4){
0113     Map_X0    = new TH2F(("Map_X0_"+surface_info.idname).c_str(),("Map_X0_"+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_L0    = new TH2F(("Map_L0_"+surface_info.idname).c_str(),("Map_L0_"+surface_info.idname).c_str(),
0116                           50,-1*surface_info.range_max,surface_info.range_max,50,-1*surface_info.range_max, surface_info.range_max);
0117     Map_scale = new TH2F(("Map_scale_"+surface_info.idname).c_str(),("Map_scale_"+surface_info.idname).c_str(),
0118                           50,-1*surface_info.range_max,surface_info.range_max,50,-1*surface_info.range_max, surface_info.range_max);
0119     Map_X0->GetXaxis()->SetTitle("X [mm]");
0120     Map_X0->GetYaxis()->SetTitle("Y [mm]");
0121     Map_X0->GetZaxis()->SetTitle("X0");
0122     Map_L0->GetXaxis()->SetTitle("X [mm]");
0123     Map_L0->GetYaxis()->SetTitle("Y [mm]");
0124     Map_L0->GetZaxis()->SetTitle("L0");
0125   }
0126   std::vector<TH2F*> v_hist;
0127   v_hist.push_back(Map_X0);
0128   v_hist.push_back(Map_L0);
0129   v_hist.push_back(Map_scale);
0130   surface_hist = v_hist;
0131 }
0132 
0133 /// Fill the histograms for each surfaces.
0134 
0135 void Fill(std::map<std::uint64_t,std::vector<TH2F*>>& surface_hist,  std::map<std::uint64_t,sinfo>& surface_info,
0136   const std::string& input_file, const std::string& geometry_file, const int& nbprocess){
0137 
0138   std::map<std::string,std::string> surface_name;
0139 
0140   std::map<std::uint64_t,float> surface_weight;
0141 
0142   //Get file, tree and set top branch address
0143   TFile *tfile = new TFile(input_file.c_str());
0144   TTree *tree = (TTree*)tfile->Get("material_tracks");
0145 
0146   float v_phi   = 0;
0147   float v_eta   = 0;
0148   std::vector<float> *mat_X0   = 0;
0149   std::vector<float> *mat_L0   = 0;
0150   std::vector<float> *mat_step_length = 0;
0151 
0152   std::vector<std::uint64_t> *sur_id = 0;
0153   std::vector<float> *sur_x = 0;
0154   std::vector<float> *sur_y = 0;
0155   std::vector<float> *sur_z = 0;
0156 
0157   tree->SetBranchAddress("v_phi",&v_phi);
0158   tree->SetBranchAddress("v_eta",&v_eta);
0159   tree->SetBranchAddress("mat_X0",&mat_X0);
0160   tree->SetBranchAddress("mat_L0",&mat_L0);
0161   tree->SetBranchAddress("mat_step_length",&mat_step_length);
0162 
0163   tree->SetBranchAddress("sur_id",&sur_id);
0164   tree->SetBranchAddress("sur_x",&sur_x);
0165   tree->SetBranchAddress("sur_y",&sur_y);
0166   tree->SetBranchAddress("sur_z",&sur_z);
0167 
0168   int nentries = tree->GetEntries();
0169   if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
0170   std::unordered_map<std::uint64_t, json> surface_bounds = load_geometry_file(geometry_file);
0171 
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       // If a valid surface id was never encountered initialise the hist, info and weight
0185       if(surface_hist.find(sur_id->at(j))==surface_hist.end()){
0186         int type = -1;
0187         float range_min = 0.;
0188         float range_max = 0.;
0189         if (surface_bounds.count(sur_id->at(j))) {
0190           json &bounds = surface_bounds[sur_id->at(j)];
0191           std::string btype = bounds["type"].get<std::string>();
0192           const auto &values = bounds["values"];
0193           if (btype == "CylinderBounds") {
0194             type = 1;
0195             range_min = -values[1].get<float>();
0196             range_max = values[1].get<float>();
0197           } else if (btype == "RadialBounds") {
0198             type = 2;
0199             range_min = values[0].get<float>();
0200             range_max = values[1].get<float>();
0201           } else {
0202             type = -1;
0203           }
0204         }
0205 
0206         float pos;
0207         if(type == 1){
0208           pos = sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j));
0209         }
0210         if(type == 2 || type == 4){
0211           pos = sur_z->at(j);
0212         }
0213 
0214         // Ignore surface of incorrect type
0215         if(type == -1) continue;
0216 
0217         surface_weight[sur_id->at(j)] = 0;
0218         // Use type and new ranges in Initialise_info
0219         Initialise_info(surface_info[sur_id->at(j)], surface_name, sur_id->at(j), type, pos, range_min, range_max);
0220         Initialise_hist(surface_hist[sur_id->at(j)], surface_info[sur_id->at(j)]);
0221       }
0222       // Weight for each surface = number of hit associated to it.
0223       surface_weight[sur_id->at(j)]++;
0224     }
0225 
0226     // loop over all the material hit to fill the histogram
0227     for(int j=0; j<mat_X0->size(); j++ ){
0228 
0229       int type = surface_info[sur_id->at(j)].type;
0230 
0231       // Ignore surface of incorrect type
0232       if(type == -1) continue;
0233 
0234       if(type == 1){
0235         surface_hist[sur_id->at(j)][0]->Fill(v_eta, v_phi, (mat_step_length->at(j)/mat_X0->at(j)));
0236         surface_hist[sur_id->at(j)][1]->Fill(v_eta, v_phi, (mat_step_length->at(j)/mat_L0->at(j)));
0237         surface_hist[sur_id->at(j)][2]->Fill(v_eta, v_phi, (1/surface_weight[sur_id->at(j)]));
0238       }
0239       if(type == 2 || type == 4){
0240         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)));
0241         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)));
0242         surface_hist[sur_id->at(j)][2]->Fill(sur_x->at(j), sur_y->at(j), (1/surface_weight[sur_id->at(j)]));
0243       }
0244     }
0245   }
0246   // Normalise the histograms
0247   for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
0248     hist_it->second[0]->Divide(hist_it->second[2]);
0249     hist_it->second[1]->Divide(hist_it->second[2]);
0250   }
0251 }
0252 
0253 /// Plot the material on each surface.
0254 /// nbprocess : number of parameter to be processed.
0255 /// name : name of the output directory.
0256 
0257 void Mat_map_surface_plot(std::string input_file = "", int nbprocess = -1, std::string name = "", std::string geometry_file = "geometry-map.json"){
0258 
0259   gStyle->SetOptStat(0);
0260   gStyle->SetOptTitle(0);
0261 
0262   std::map<std::uint64_t,std::vector<TH2F*>> surface_hist;
0263   std::map<std::uint64_t,sinfo> surface_info;
0264 
0265   Fill(surface_hist, surface_info, input_file, geometry_file, nbprocess);
0266   for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
0267     plot(hist_it->second, surface_info[hist_it->first], name);
0268     for (auto hist : hist_it->second){
0269       delete hist;
0270     }
0271     hist_it->second.clear();
0272   }
0273 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!