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