Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:12:12 

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 //
0010 //  layerMaterial.C
0011 //
0012 //
0013 //  Created by Julia Hrdinka on 18/08/16.
0014 //
0015 //
0016 
0017 #include <algorithm>
0018 #include <iostream>
0019 #include <vector>
0020 #include "TFile.h"
0021 #include "TH1F.h"
0022 #include "TObject.h"
0023 #include "TROOT.h"
0024 #include "TTree.h"
0025 
0026 // This root script generates material histograms and maps of a layer.
0027 // It is foreseen to be used for output generated by the RootMaterialWriter or
0028 // the RootMaterialStepWriter. For the indicated Layer it produces the
0029 // histograms and maps of the material properties along the local coordinates of
0030 // the layer.
0031 // It also writes the assigned and real global positions into histograms, if
0032 // given.
0033 
0034 void
0035 layerMaterial(std::string inFile,
0036               std::string outFile,
0037               int         binsLoc1,
0038               int         binsLoc2,
0039               int         binsZ,
0040               float       minGlobZ,
0041               float       maxGlobZ,
0042               int         binsR,
0043               float       minGlobR,
0044               float       maxGlobR,
0045               std::string layerName)
0046 {
0047   std::cout << "Opening file: " << inFile << std::endl;
0048   TFile  inputFile(inFile.c_str());
0049   TTree* layer = (TTree*)inputFile.Get(layerName.c_str());
0050   std::cout << "Reading tree: " << layerName << std::endl;
0051 
0052   std::vector<float>* loc0  = new std::vector<float>;
0053   std::vector<float>* loc1  = new std::vector<float>;
0054   std::vector<float>* A     = new std::vector<float>;
0055   std::vector<float>* Z     = new std::vector<float>;
0056   std::vector<float>* x0    = new std::vector<float>;
0057   std::vector<float>* l0    = new std::vector<float>;
0058   std::vector<float>* d     = new std::vector<float>;
0059   std::vector<float>* rho   = new std::vector<float>;
0060   std::vector<float>* dInX0 = new std::vector<float>;
0061   std::vector<float>* dInL0 = new std::vector<float>;
0062   std::vector<float>* t     = new std::vector<float>;
0063 
0064   std::vector<float>* globX = new std::vector<float>;
0065   std::vector<float>* globY = new std::vector<float>;
0066   std::vector<float>* globZ = new std::vector<float>;
0067   std::vector<float>* globR = new std::vector<float>;
0068 
0069   std::vector<float>* assignedGlobX = new std::vector<float>;
0070   std::vector<float>* assignedGlobY = new std::vector<float>;
0071   std::vector<float>* assignedGlobZ = new std::vector<float>;
0072   std::vector<float>* assignedGlobR = new std::vector<float>;
0073 
0074   layer->SetBranchAddress("loc0", &loc0);
0075   layer->SetBranchAddress("loc1", &loc1);
0076   layer->SetBranchAddress("A", &A);
0077   layer->SetBranchAddress("Z", &Z);
0078   layer->SetBranchAddress("x0", &x0);
0079   layer->SetBranchAddress("l0", &l0);
0080   layer->SetBranchAddress("thickness", &d);
0081   layer->SetBranchAddress("rho", &rho);
0082   layer->SetBranchAddress("tInX0", &dInX0);
0083   layer->SetBranchAddress("tInL0", &dInL0);
0084   layer->SetBranchAddress("thickness", &t);
0085 
0086   if (layer->FindBranch("globX") && layer->FindBranch("globY")
0087       && layer->FindBranch("globZ") && layer->FindBranch("globR")) {
0088     layer->SetBranchAddress("globX", &globX);
0089     layer->SetBranchAddress("globY", &globY);
0090     layer->SetBranchAddress("globZ", &globZ);
0091     layer->SetBranchAddress("globR", &globR);
0092   }
0093 
0094   if (layer->FindBranch("assignedGlobX") && layer->FindBranch("assignedGlobY")
0095       && layer->FindBranch("assignedGlobZ")
0096       && layer->FindBranch("assignedGlobR")) {
0097     layer->SetBranchAddress("assignedGlobX", &assignedGlobX);
0098     layer->SetBranchAddress("assignedGlobY", &assignedGlobY);
0099     layer->SetBranchAddress("assignedGlobZ", &assignedGlobZ);
0100     layer->SetBranchAddress("assignedGlobR", &assignedGlobR);
0101   }
0102   layer->GetEntry(0);
0103 
0104   // get minima and maxima for different layers
0105   auto  minmax0 = std::minmax_element(loc0->begin(), loc0->end());
0106   float min0    = *minmax0.first;
0107   float max0    = *minmax0.second;
0108 
0109   auto  minmax1 = std::minmax_element(loc1->begin(), loc1->end());
0110   float min1    = *minmax1.first;
0111   float max1    = *minmax1.second;
0112 
0113   inputFile.Close();
0114   std::cout << "Creating new output file: " << outFile
0115             << " and writing "
0116                "material maps"
0117             << std::endl;
0118   TFile       outputFile(outFile.c_str(), "update");
0119   TDirectory* dir = outputFile.mkdir(layerName.c_str());
0120   dir->cd();
0121   // thickness in X0
0122   TProfile* dInX0_loc1
0123       = new TProfile("dInX0_loc1", "dInX0_loc1", binsLoc1, min1, max1);
0124   TProfile* dInX0_loc0
0125       = new TProfile("dInX0_loc0", "dInX0_loc0", binsLoc1, min0, max0);
0126   TProfile2D* dInX0_map = new TProfile2D(
0127       "dInX0", "dInX0", binsLoc1, min1, max1, binsLoc1, min0, max0);
0128   // thickness in L0
0129   TProfile* dInL0_loc1
0130       = new TProfile("dInL0_loc1", "dInL0_loc1", binsLoc2, min1, max1);
0131   TProfile* dInL0_loc0
0132       = new TProfile("dInL0_loc0", "dInL0_loc0", binsLoc1, min0, max0);
0133   TProfile2D* dInL0_map = new TProfile2D(
0134       "dInL0", "dInL0", binsLoc2, min1, max1, binsLoc1, min0, max0);
0135   // A
0136   TProfile*   A_loc1 = new TProfile("A_loc1", "A_loc1", binsLoc2, min1, max1);
0137   TProfile*   A_loc0 = new TProfile("A_loc0", "A_loc0", binsLoc1, min0, max0);
0138   TProfile2D* A_map
0139       = new TProfile2D("A", "A", binsLoc2, min1, max1, binsLoc1, min0, max0);
0140   // Z
0141   TProfile*   Z_loc1 = new TProfile("Z_loc1", "Z_loc1", binsLoc2, min1, max1);
0142   TProfile*   Z_loc0 = new TProfile("Z_loc0", "Z_loc0", binsLoc1, min0, max0);
0143   TProfile2D* Z_map
0144       = new TProfile2D("Z", "Z", binsLoc2, min1, max1, binsLoc1, min0, max0);
0145   // Rho
0146   TProfile* rho_loc1
0147       = new TProfile("rho_loc1", "rho_loc1", binsLoc2, min1, max1);
0148   TProfile* rho_loc0
0149       = new TProfile("rho_loc0", "rho_loc0", binsLoc1, min0, max0);
0150   TProfile2D* rho_map = new TProfile2D(
0151       "rho", "rho", binsLoc2, min1, max1, binsLoc1, min0, max0);
0152   // x0
0153   TProfile* x0_loc1 = new TProfile("x0_loc1", "x0_loc1", binsLoc2, min1, max1);
0154   TProfile* x0_loc0 = new TProfile("x0_loc0", "x0_loc0", binsLoc1, min0, max0);
0155   TProfile2D* x0_map
0156       = new TProfile2D("x0", "x0", binsLoc2, min1, max1, binsLoc1, min0, max0);
0157   // l0
0158   TProfile* l0_loc1 = new TProfile("l0_loc1", "l0_loc1", binsLoc2, min1, max1);
0159   TProfile* l0_loc0 = new TProfile("l0_loc0", "l0_loc0", binsLoc1, min0, max0);
0160   TProfile2D* l0_map
0161       = new TProfile2D("l0", "l0", binsLoc2, min1, max1, binsLoc1, min0, max0);
0162 
0163   // thickness
0164   TProfile*   t_loc1 = new TProfile("t_loc1", "t_loc1", binsLoc2, min1, max1);
0165   TProfile*   t_loc0 = new TProfile("t_loc0", "t_loc0", binsLoc1, min0, max0);
0166   TProfile2D* t_map
0167       = new TProfile2D("t", "t", binsLoc2, min1, max1, binsLoc1, min0, max0);
0168 
0169   // global
0170   TH2F* glob_r_z = new TH2F(
0171       "r_z", "r_z", binsZ, minGlobZ, maxGlobZ, binsR, minGlobR, maxGlobR);
0172   TH2F* assigned_r_z = new TH2F("r_z_assigned",
0173                                 "r_z_assigned",
0174                                 binsZ,
0175                                 minGlobZ,
0176                                 maxGlobZ,
0177                                 binsR,
0178                                 minGlobR,
0179                                 maxGlobR);
0180 
0181   TH2F* glob_x_y = new TH2F(
0182       "x_y", "x_y", binsR, -maxGlobR, maxGlobR, binsR, -maxGlobR, maxGlobR);
0183   TH2F* assigned_x_y = new TH2F("x_y_assigned",
0184                                 "x_y_assigned",
0185                                 binsR,
0186                                 -maxGlobR,
0187                                 maxGlobR,
0188                                 binsR,
0189                                 -maxGlobR,
0190                                 maxGlobR);
0191 
0192   std::size_t nEntries = loc1->size();
0193   for (int i = 0; i < nEntries; i++) {
0194     // A
0195     A_loc1->Fill(loc1->at(i), A->at(i));
0196     A_loc0->Fill(loc0->at(i), A->at(i));
0197     A_map->Fill(loc1->at(i), loc0->at(i), A->at(i));
0198     // Z
0199     Z_loc1->Fill(loc1->at(i), Z->at(i));
0200     Z_loc0->Fill(loc0->at(i), Z->at(i));
0201     Z_map->Fill(loc1->at(i), loc0->at(i), Z->at(i));
0202     // x0
0203     x0_loc1->Fill(loc1->at(i), x0->at(i));
0204     x0_loc0->Fill(loc0->at(i), x0->at(i));
0205     x0_map->Fill(loc1->at(i), loc0->at(i), x0->at(i));
0206     // l0
0207     l0_loc1->Fill(loc1->at(i), l0->at(i));
0208     l0_loc0->Fill(loc0->at(i), l0->at(i));
0209     l0_map->Fill(loc1->at(i), loc0->at(i), l0->at(i));
0210     // rho
0211     rho_loc1->Fill(loc1->at(i), rho->at(i));
0212     rho_loc0->Fill(loc0->at(i), rho->at(i));
0213     rho_map->Fill(loc1->at(i), loc0->at(i), rho->at(i));
0214     // thickness in X0
0215     dInX0_loc1->Fill(loc1->at(i), dInX0->at(i));
0216     dInX0_loc0->Fill(loc0->at(i), dInX0->at(i));
0217     dInX0_map->Fill(loc1->at(i), loc0->at(i), dInX0->at(i));
0218     // thickness in L0
0219     dInL0_loc1->Fill(loc1->at(i), dInL0->at(i));
0220     dInL0_loc0->Fill(loc0->at(i), dInL0->at(i));
0221     dInL0_map->Fill(loc1->at(i), loc0->at(i), dInL0->at(i));
0222     // thickness
0223     t_loc1->Fill(loc1->at(i), t->at(i));
0224     t_loc0->Fill(loc0->at(i), t->at(i));
0225     t_map->Fill(loc1->at(i), loc0->at(i), t->at(i));
0226 
0227     // fill global r/z
0228     if (globZ->size() && globR->size())
0229       glob_r_z->Fill(globZ->at(i), globR->at(i));
0230 
0231     if (assignedGlobZ->size() && assignedGlobR->size())
0232       assigned_r_z->Fill(assignedGlobZ->at(i), assignedGlobR->at(i));
0233 
0234     // fill global x/y
0235     if (globX->size() && globY->size())
0236       glob_x_y->Fill(globX->at(i), globY->at(i));
0237 
0238     if (assignedGlobX->size() && assignedGlobY->size())
0239       assigned_x_y->Fill(assignedGlobX->at(i), assignedGlobY->at(i));
0240   }
0241 
0242   gStyle->SetOptStat(0);
0243   // A
0244   A_loc1->Write();
0245   A_loc0->Write();
0246   A_map->Write();
0247   delete A_loc1;
0248   delete A_loc0;
0249   delete A_map;
0250   // Z
0251   Z_loc1->Write();
0252   Z_loc0->Write();
0253   Z_map->Write();
0254   delete Z_loc1;
0255   delete Z_loc0;
0256   delete Z_map;
0257   // x0
0258   x0_loc1->Write();
0259   x0_loc0->Write();
0260   x0_map->Write();
0261   delete x0_loc1;
0262   delete x0_loc0;
0263   delete x0_map;
0264   // l0
0265   l0_loc1->Write();
0266   l0_loc0->Write();
0267   l0_map->Write();
0268   delete l0_loc1;
0269   delete l0_loc0;
0270   delete l0_map;
0271   // rho
0272   rho_loc1->Write();
0273   rho_loc0->Write();
0274   rho_map->Write();
0275   delete rho_loc1;
0276   delete rho_loc0;
0277   delete rho_map;
0278   // thickness in X0
0279   dInX0_loc1->Write();
0280   dInX0_loc0->Write();
0281   dInX0_map->Write();
0282   delete dInX0_loc1;
0283   delete dInX0_loc0;
0284   delete dInX0_map;
0285   // thickness in L0
0286   dInL0_loc1->Write();
0287   dInL0_loc0->Write();
0288   dInL0_map->Write();
0289   delete dInL0_loc1;
0290   delete dInL0_loc0;
0291   delete dInL0_map;
0292   // thickness
0293   t_loc1->Write();
0294   t_loc0->Write();
0295   t_map->Write();
0296 
0297   delete loc0;
0298   delete loc1;
0299   delete A;
0300   delete Z;
0301   delete x0;
0302   delete l0;
0303   delete d;
0304   delete rho;
0305   delete dInX0;
0306   delete dInL0;
0307   delete t;
0308 
0309   glob_r_z->Write();
0310   assigned_r_z->Write();
0311   glob_x_y->Write();
0312   assigned_x_y->Write();
0313   delete glob_r_z;
0314   delete assigned_r_z;
0315   delete glob_x_y;
0316   delete assigned_x_y;
0317   outputFile.Close();
0318 }
[ 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!