nhcal_basic_distribution/scripts/basic_distribution_analysis.cxx

0001 #include <cmath>
0002 #include <fstream>
0003 #include <iostream>
0004 #include <string>
0005 #include <vector>
0006
0007 #include "ROOT/RDataFrame.hxx"
0008 #include <TH1D.h>
0009 #include <TH2D.h>
0010 #include <TRandom3.h>
0011 #include <TFile.h>
0012 #include <TChain.h>
0013 #include <TTree.h>
0014 #include <TMath.h>
0015 #include <TVector3.h>
0016 #include <TCanvas.h>
0017
0018 #include "TROOT.h"
0019 #include "TRandom.h"
0020 #include "TH3.h"
0021
0022
0023 #include "DD4hep/Detector.h"
0024 #include "DDRec/CellIDPositionConverter.h"
0025
0026 #include <podio/Frame.h>
0027 #include <podio/CollectionBase.h>
0028 #include "podio/ROOTReader.h"
0029 //#include <podio/ROOTFrameReader.h>
0030 #include "podio/CollectionIDTable.h"
0031 #include "podio/ObjectID.h"
0032
0033 #include "edm4hep/MCParticleCollection.h"
0034 #include "edm4hep/MCParticleCollectionData.h"
0035 #include "edm4hep/MCParticle.h"
0036 #include "edm4hep/MCParticleData.h"
0037
0038 #include "edm4hep/SimCalorimeterHitCollectionData.h"
0039 #include "edm4hep/SimCalorimeterHitCollection.h"
0040 #include "edm4hep/SimCalorimeterHitData.h"
0041 #include "edm4hep/SimCalorimeterHit.h"
0042
0043 #include "edm4hep/CalorimeterHit.h"
0044 #include "edm4hep/CalorimeterHitCollectionData.h"
0045 #include "edm4hep/CalorimeterHitCollection.h"
0046 #include "edm4hep/CalorimeterHitData.h"
0047 #include "edm4hep/CalorimeterHit.h"
0048 #include "edm4hep/CalorimeterHitObj.h"
0049
0050 #include "edm4eic/ClusterCollection.h"
0051 #include "edm4eic/Cluster.h"
0052 #include "edm4eic/ClusterData.h"
0053
0054 #include "edm4eic/CalorimeterHit.h"
0055 #include "edm4eic/CalorimeterHitCollectionData.h"
0056 #include "edm4eic/CalorimeterHitCollection.h"
0057 #include "edm4eic/CalorimeterHitData.h"
0058 #include "edm4eic/CalorimeterHit.h"
0059 #include "edm4eic/CalorimeterHitObj.h"
0060
0061
0062 #include <edm4eic/vector_utils_legacy.h>
0063 #include <edm4hep/Vector3f.h>
0064
0065 // #include "FileList.h"
0066 // #include "EICutil.h"
0067 // #include "BasicUtil.h"
0068
0069 // #pragma link C++ class vector<edm4hep::MCParticleData>+;
0070 // #pragma link C++ class vector<podio::ObjectID>+;
0071 // #pragma link C++ class vector<edm4hep::SimCalorimeterHitData>+;
0072
0073 using namespace std;
0074 using namespace ROOT;
0075 using namespace TMath;
0076 using namespace edm4hep;
0077
0078 int basic_distribution_analysis(const string &filename, string outname_pdf, string outname_png)
0079 {
0080     podio::ROOTReader *reader = new podio::ROOTReader();
0081     reader->openFile(filename);
0082     unsigned nEvents = reader->getEntries("events");
0083     cout << "Number of events: " << nEvents << endl;
0084
0085     vector<double> xPosition;
0086     vector<double> yPosition;
0087     vector<double> zPosition;
0088     vector<double> rPosition;
0089     vector<double> energy;
0090     vector<double> totalEventEnergy;
0091     vector<double> totalEventHits;
0092     vector<double> zLayerValues;
0093     vector<int> nHitsPerLayer;
0094     vector<double> nEnergyPerLayer;
0095
0096     for (unsigned ev = 0; ev < nEvents; ev++)
0097     {
0098         auto frameData = reader->readNextEntry(podio::Category::Event);
0099         if (!frameData)
0100         {
0101             cerr << "Invalid FrameData at event " << ev << endl;
0102             continue;
0103         }
0104
0105         podio::Frame frame(std::move(frameData));
0106
0107         auto& hits = frame.get<edm4hep::SimCalorimeterHitCollection>("HcalEndcapNHits");
0108         if (!hits.isValid())
0109         {
0110             cerr << "HcalEndcapNHits collection is not valid!" << endl;
0111         }
0112
0113         map<double, pair<int, double>> layerData;
0114         double totalEnergy = 0;
0115         int totalHits = 0;
0116
0117         for (const auto& hit : hits)
0118         {
0119             totalHits++;
0120             energy.push_back(hit.getEnergy());
0121             totalEnergy += hit.getEnergy();
0122
0123             auto pos = hit.getPosition();
0124             xPosition.push_back(pos.x);
0125             yPosition.push_back(pos.y);
0126             zPosition.push_back(pos.z);
0127             rPosition.push_back(sqrt(pos.x * pos.x + pos.y * pos.y));
0128
0129             double zBin = round(pos.z);
0130             layerData[zBin].first++;
0131             layerData[zBin].second += hit.getEnergy();
0132         }
0133
0134         totalEventEnergy.push_back(totalEnergy);
0135         totalEventHits.push_back(totalHits);
0136
0137         for (const auto& [zValue, stats] : layerData)
0138         {
0139             zLayerValues.push_back(zValue);
0140             nHitsPerLayer.push_back(stats.first);
0141             nEnergyPerLayer.push_back(stats.second);
0142         }
0143
0144     }
0145
0146     auto minmax_xPosition        = minmax_element(xPosition.begin(), xPosition.end());
0147     auto minmax_yPosition        = minmax_element(yPosition.begin(), yPosition.end());
0148     auto minmax_zPosition        = minmax_element(zPosition.begin(), zPosition.end());
0149     auto minmax_rPosition        = minmax_element(rPosition.begin(), rPosition.end());
0150     auto minmax_totalEventEnergy = minmax_element(totalEventEnergy.begin(), totalEventEnergy.end());
0151     auto minmax_totalEventHits   = minmax_element(totalEventHits.begin(), totalEventHits.end());
0152     auto minmax_energy           = minmax_element(energy.begin(), energy.end());
0153     auto minmax_zLayerValues     = minmax_element(zLayerValues.begin(), zLayerValues.end());
0154     auto minmax_nHitsPerLayer    = minmax_element(nHitsPerLayer.begin(), nHitsPerLayer.end());
0155     auto minmax_nEnergyPerLayer  = minmax_element(nEnergyPerLayer.begin(), nEnergyPerLayer.end());
0156
0157     int nbins = nEvents/1000;
0158
0159     TH1D *h_energyTotal = new TH1D("h_energyTotal", "Total Energy per Event; Energy per Event",
0160                                     nbins/2, *minmax_totalEventEnergy.first, *minmax_totalEventEnergy.second);
0161     TH2D *h_layerEnergy = new TH2D("h_layerEnergy", "Energy in Layers (Z); Z; Energy",
0162                                     nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second,
0163                                     nbins/4, *minmax_nEnergyPerLayer.first, *minmax_nEnergyPerLayer.second);
0164     TProfile *p_layerEnergy = new TProfile("p_layerEnergy", "Energy in Layers (Z); Z; Mean Energy",
0165                                     nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second);
0166     TH1D *h_hitCount    = new TH1D("h_hitCount", "Number of Hits per Event; Hits per Event",
0167                                     nbins/2, *minmax_totalEventHits.first, *minmax_totalEventHits.second);
0168     TH2D *h_layerHits   = new TH2D("h_layerHits", "Number of Hits in Layers (Z); Layer(Z); Hits",
0169                                     nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second,
0170                                     *minmax_nHitsPerLayer.second, *minmax_nHitsPerLayer.first, *minmax_nHitsPerLayer.second);
0171     TProfile *p_layerHits = new TProfile("p_layerHits", "Number of Hits in Layers (Z); Z; Mean Hits",
0172                                     nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second);
0173     TH2D *h_XYPos       = new TH2D("h_XYPos", "Hits position X,Y;X [mm];Y [mm]",
0174                                     nbins, *minmax_xPosition.first, *minmax_xPosition.second,
0175                                     nbins, *minmax_yPosition.first, *minmax_yPosition.second);
0176     TH2D *h_ZRPos       = new TH2D("h_ZRPos", "Hits position Z,R;Z [mm];R [mm]",
0177                                     nbins/5, *minmax_zPosition.first, *minmax_zPosition.second,
0178                                     nbins, *minmax_rPosition.first, *minmax_rPosition.second);
0179     TH2D *h_XYEnergy    = new TH2D("h_XYEnergy", "Hits energy X,Y;X [mm];Y [mm]",
0180                                     nbins, *minmax_xPosition.first, *minmax_xPosition.second,
0181                                     nbins, *minmax_yPosition.first, *minmax_yPosition.second);
0182
0183     for(int i = 0; i < xPosition.size(); i++)
0184     {
0185
0186         h_XYPos->Fill(xPosition[i], yPosition[i]);
0187         h_ZRPos->Fill(zPosition[i], rPosition[i]);
0188         h_XYEnergy->Fill(xPosition[i], yPosition[i], energy[i]);
0189     }
0190
0191     for (int i = 0; i < zLayerValues.size(); i++)
0192     {
0193         h_layerHits->Fill(zLayerValues[i], nHitsPerLayer[i]);
0194         p_layerHits->Fill(zLayerValues[i], nHitsPerLayer[i]);
0195
0196         h_layerEnergy->Fill(zLayerValues[i], nEnergyPerLayer[i]);
0197         p_layerEnergy->Fill(zLayerValues[i], nEnergyPerLayer[i]);
0198     }
0199
0200     for(int i = 0; i < nEvents; i++)
0201     {
0202         h_energyTotal->Fill(totalEventEnergy[i]);
0203         h_hitCount->Fill(totalEventHits[i]);
0204     }
0205
0206     gStyle->SetPadLeftMargin(0.13);
0207
0208     TCanvas *canvas = new TCanvas("canvas", "canvas", 1600, 800);
0209     canvas->Divide(4,2);
0210     canvas->cd(1);
0211     h_energyTotal->Draw();
0212     canvas->cd(2);
0213     h_layerEnergy->Draw("COLZ");
0214     p_layerEnergy->Draw("SAME");
0215     canvas->cd(3);
0216     h_hitCount->Draw();
0217     canvas->cd(4);
0218     h_layerHits->Draw("COLZ");
0219     p_layerHits->Draw("SAME");
0220     canvas->cd(5);
0221     h_XYPos->Draw("COLZ");
0222     canvas->cd(6);
0223     h_ZRPos->Draw("COLZ");
0224     canvas->cd(7);
0225     h_XYEnergy->Draw("COLZ");
0226
0227     canvas->SaveAs(outname_pdf.c_str());
0228     canvas->SaveAs(outname_png.c_str());
0229
0230     return 0;
0231 }