Skip to content

Commit

Permalink
migrate to podio::Frame based datamodel
Browse files Browse the repository at this point in the history
  • Loading branch information
@SanghyunKo
SanghyunKo committed Jul 6, 2023
1 parent 0a229d4 commit 73c571f
Showing 1 changed file with 81 additions and 54 deletions.
135 changes: 81 additions & 54 deletions analysis/analysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@
#include "edm4hep/TimeSeriesCollection.h"
#include "edm4hep/RawTimeSeriesCollection.h"

#include "podio/ROOTReader.h"
#include "podio/EventStore.h"
#include "podio/ROOTFrameReader.h"
#include "podio/Frame.h"

#include "TROOT.h"
#include "TStyle.h"
Expand All @@ -21,70 +21,91 @@
#include "TLorentzVector.h"
#include "TGraph.h"

#include "Math/Vector3Dfwd.h"
#include "Math/Vector3D.h"

#include <iostream>
#include <string>

int main(int , char* argv[]) {
TString filename = argv[1];
float low = std::stof(argv[2]);
float high = std::stof(argv[3]);
int main(int argc, char* argv[]) {
float low = std::stof(argv[1]);
float high = std::stof(argv[2]);
std::vector<std::string> filenames;

for (int idx = 3; idx < argc; idx++)
filenames.push_back( std::string(argv[idx]) );

gStyle->SetOptFit(1);

TH1F* tEdep = new TH1F("totEdep","Total Energy deposit;GeV;Evt",100,low,high);
tEdep->Sumw2(); tEdep->SetLineColor(kRed); tEdep->SetLineWidth(2);
TH1F* tE_C = new TH1F("E_C","Energy of Cerenkov ch.;GeV;Evt",100,low,high);
TH1F* tE_C = new TH1F("E_C","Energy of Cerenkov ch.;GeV;Evt",80,low,high);
tE_C->Sumw2(); tE_C->SetLineColor(kBlue); tE_C->SetLineWidth(2);
TH1F* tE_S = new TH1F("E_S","Energy of Scintillation ch.;GeV;Evt",100,low,high);
TH1F* tE_S = new TH1F("E_S","Energy of Scintillation ch.;GeV;Evt",80,low,high);
tE_S->Sumw2(); tE_S->SetLineColor(kRed); tE_S->SetLineWidth(2);
TH1F* tE_SC = new TH1F("E_SC","E_{S}+E_{C};GeV;Evt",100,2.*low,2.*high);
TH1F* tE_SC = new TH1F("E_SC","E_{S}+E_{C};GeV;Evt",80,2.*low,2.*high);
tE_SC->Sumw2(); tE_SC->SetLineColor(kBlack); tE_SC->SetLineWidth(2);
TH1F* tE_DR = new TH1F("E_DR","Dual-readout corrected Energy;GeV;Evt",100,low,high);
TH1F* tE_DR = new TH1F("E_DR","Dual-readout corrected Energy;GeV;Evt",80,low,high);
tE_DR->Sumw2(); tE_DR->SetLineColor(kBlack); tE_DR->SetLineWidth(2);

functions::dualhist1D* tT_dual = new functions::dualhist1D("time","time;ns;p.e.",150,10.,70.);
functions::dualhist1D* tWav_dual = new functions::dualhist1D("wavlen","wavelength;nm;p.e.",120,300.,900.);
functions::dualhist1D* tNhit_dual = new functions::dualhist1D("nHits","of Scint p.e./SiPM;p.e.;n",100,0.,100.);
functions::dualhist1D* tD_dual = new functions::dualhist1D("digi","digi waveform;ns;a.u.",300,10.,100.);
functions::dualhist1D* tToa_dual = new functions::dualhist1D("ToA","Time of Arrival;ns;a.u.",150,10.,70.);
functions::dualhist1D* tInt_dual = new functions::dualhist1D("Int","ADC integration;ADC counts;a.u.",400,0.,2000.);
functions::dualhist1D* tProc_dual = new functions::dualhist1D("proc","processed waveform;ns;a.u.",300,10.,100.);
functions::dualhist1D* tInt_dual = new functions::dualhist1D("Int","ADC integration;ADC counts;a.u.",500,0.,500.);
functions::dualhist1D* tProc_dual = new functions::dualhist1D("proc","processed waveform;ns;a.u.",475,5.,100.);

auto pReader = std::make_unique<podio::ROOTReader>();
pReader->openFile(static_cast<std::string>(filename));
std::vector<float> E_Ss,E_Cs;

auto pStore = std::make_unique<podio::EventStore>();
pStore->setReader(pReader.get());
auto reader = podio::ROOTFrameReader();
reader.openFiles(filenames);

std::vector<float> E_Ss,E_Cs;
unsigned int entries = reader.getEntries("events");

unsigned int entries = pReader->getEntries();
for (unsigned int iEvt = 0; iEvt < entries; iEvt++) {
if (iEvt % 100 == 0) printf("Analyzing %dth event ...\n", iEvt);

auto& edepHits = pStore->get<edm4hep::SimCalorimeterHitCollection>("SimCalorimeterHits");
auto& rawTimeStructs = pStore->get<edm4hep::RawTimeSeriesCollection>("RawTimeStructs");
auto& rawWavlenStructs = pStore->get<edm4hep::RawTimeSeriesCollection>("RawWavlenStructs");
auto& digiWaveforms = pStore->get<edm4hep::TimeSeriesCollection>("DigiWaveforms");
auto& caloHits = pStore->get<edm4hep::CalorimeterHitCollection>("DRcalo2dHits");
auto& rawHits = pStore->get<edm4hep::RawCalorimeterHitCollection>("RawCalorimeterHits");
auto& digiHits = pStore->get<edm4hep::RawCalorimeterHitCollection>("DigiCalorimeterHits");
auto& procTimes = pStore->get<edm4hep::TimeSeriesCollection>("DRpostprocTime");
auto aframe = podio::Frame(reader.readEntry("events",iEvt));

const auto* edepHits = static_cast<const edm4hep::SimCalorimeterHitCollection*>(aframe.get("SimCalorimeterHits"));
const auto* rawTimeStructs = static_cast<const edm4hep::RawTimeSeriesCollection*>(aframe.get("RawTimeStructs"));
const auto* rawWavlenStructs = static_cast<const edm4hep::RawTimeSeriesCollection*>(aframe.get("RawWavlenStructs"));
const auto* digiWaveforms = static_cast<const edm4hep::TimeSeriesCollection*>(aframe.get("DigiWaveforms"));
const auto* caloHits = static_cast<const edm4hep::CalorimeterHitCollection*>(aframe.get("DRcalo2dHits"));
const auto* rawHits = static_cast<const edm4hep::RawCalorimeterHitCollection*>(aframe.get("RawCalorimeterHits"));
const auto* digiHits = static_cast<const edm4hep::CalorimeterHitCollection*>(aframe.get("DigiCalorimeterHits"));
const auto* procTimes = static_cast<const edm4hep::TimeSeriesCollection*>(aframe.get("DRpostprocTime"));
const auto* leakages = static_cast<const edm4hep::MCParticleCollection*>(aframe.get("Leakages"));

float leakage = 0.;

for (unsigned int iLeak = 0; iLeak < leakages->size(); iLeak++) {
auto& momentum = (*leakages)[iLeak].getMomentum();
ROOT::Math::XYZVectorF vec(momentum.x, momentum.y, momentum.z);
leakage += std::sqrt(vec.Mag2());
}

if (leakage > 5.)
continue;

float Edep = 0.;
for (unsigned int iEdep = 0; iEdep < edepHits.size(); iEdep++)
Edep += edepHits[iEdep].getEnergy();

for (unsigned int iEdep = 0; iEdep < edepHits->size(); iEdep++)
Edep += (*edepHits)[iEdep].getEnergy();

tEdep->Fill(Edep);

float en_S = 0.; float en_C = 0.;
for (unsigned int idx = 0; idx < caloHits.size(); idx++) {
const auto& caloHit = caloHits.at(idx);
const auto& rawHit = rawHits.at(idx);
const auto& digiHit = digiHits.at(idx);
const auto& timeStruct = rawTimeStructs.at(idx);
const auto& wavlenStruct = rawWavlenStructs.at(idx);
const auto& waveform = digiWaveforms.at(idx);
const auto& procTime = procTimes.at(idx);

for (unsigned int idx = 0; idx < caloHits->size(); idx++) {
const auto& caloHit = caloHits->at(idx);
const auto& rawHit = rawHits->at(idx);
const auto& digiHit = digiHits->at(idx);
const auto& timeStruct = rawTimeStructs->at(idx);
const auto& wavlenStruct = rawWavlenStructs->at(idx);
const auto& waveform = digiWaveforms->at(idx);
const auto& procTime = procTimes->at(idx);

int type = caloHit.getType();
float en = caloHit.getEnergy();
Expand All @@ -103,27 +124,27 @@ int main(int , char* argv[]) {
tNhit->Fill(nhits);
tToa->Fill(caloHit.getTime());

if ( digiHit.getAmplitude() > 0 )
tInt->Fill(static_cast<double>(digiHit.getAmplitude()));
if ( digiHit.getEnergy() > 0 )
tInt->Fill(static_cast<double>(digiHit.getEnergy()));

for (unsigned int bin = 0; bin < timeStruct.adcCounts_size(); bin++) {
int con = timeStruct.getAdcCounts(bin);
tT->Fill(timeStruct.getTime()+static_cast<float>(bin)*timeStruct.getInterval(),static_cast<float>(con));
tT->Fill(timeStruct.getTime()+(static_cast<float>(bin)+0.5)*timeStruct.getInterval(),static_cast<float>(con));
}

for (unsigned int bin = 0; bin < wavlenStruct.adcCounts_size(); bin++) {
int con = wavlenStruct.getAdcCounts(bin);
tWav->Fill(wavlenStruct.getTime()+static_cast<float>(bin)*wavlenStruct.getInterval(),static_cast<float>(con));
tWav->Fill(wavlenStruct.getTime()+(static_cast<float>(bin)+0.5)*wavlenStruct.getInterval(),static_cast<float>(con));
}

for (unsigned int bin = 0; bin < waveform.amplitude_size(); bin++) {
float con = waveform.getAmplitude(bin);
tD->Fill(waveform.getTime()+static_cast<float>(bin)*waveform.getInterval(),con);
tD->Fill(waveform.getTime()+(static_cast<float>(bin)+0.5)*waveform.getInterval(),con);
}

for (unsigned int bin = 0; bin < procTime.amplitude_size(); bin++) {
float con = procTime.getAmplitude(bin);
tProc->Fill(procTime.getTime()+static_cast<float>(bin)*procTime.getInterval(),con);
tProc->Fill(procTime.getTime()+(static_cast<float>(bin)+0.5)*procTime.getInterval(),con);
}
}

Expand All @@ -134,37 +155,46 @@ int main(int , char* argv[]) {
tE_C->Fill(en_C);
tE_SC->Fill(en_C+en_S);
tE_DR->Fill(functions::E_DR291(en_C,en_S));

pStore->clear();
pReader->endOfEvent();
} // event loop

std::string filenameStd = static_cast<std::string>(filename);
std::string filenameStd = static_cast<std::string>(filenames.front());
std::string extension = ".root";
auto where = filenameStd.find(extension);
if (where != std::string::npos) {
filenameStd.erase(where, extension.length());
}

filename = static_cast<TString>(filenameStd);
TString filename = static_cast<TString>(filenameStd);

TCanvas* c = new TCanvas("c","");

tEdep->Draw("Hist"); c->SaveAs(filename+"_Edep.png");

TF1* grE_DR = new TF1("Efit","gaus",low,high); grE_DR->SetLineColor(kBlack);
tE_DR->SetOption("p"); tE_DR->Fit(grE_DR,"R+&same");

c->cd();
tE_DR->SetTitle("");
tE_DR->Draw(""); c->Update();
TPaveStats* statsE_DR = (TPaveStats*)c->GetPrimitive("stats");
statsE_DR->SetName("DR");
statsE_DR->SetTextColor(kBlack);
statsE_DR->SetX1NDC(.7);
statsE_DR->SetY1NDC(.7); statsE_DR->SetY2NDC(1.);

tE_S->SetTitle("");
tE_S->Draw("Hist"); c->Update();
tE_S->Draw("Hist&sames"); c->Update();
TPaveStats* statsE_S = (TPaveStats*)c->GetPrimitive("stats");
statsE_S->SetName("Scint");
statsE_S->SetTextColor(kRed);
statsE_S->SetY1NDC(.6); statsE_S->SetY2NDC(.8);
statsE_S->SetY1NDC(.5); statsE_S->SetY2NDC(.7);

tE_C->Draw("Hist&sames"); c->Update();
TPaveStats* statsE_C = (TPaveStats*)c->GetPrimitive("stats");
statsE_C->SetName("Cerenkov");
statsE_C->SetTextColor(kBlue);
statsE_C->SetY1NDC(.8); statsE_C->SetY2NDC(1.);
statsE_C->SetY1NDC(.3); statsE_C->SetY2NDC(.5);

c->SaveAs(filename+"_EcsHist.png");

TF1* grE_C = new TF1("Cfit","gaus",low,high); grE_C->SetLineColor(kBlue);
Expand All @@ -189,14 +219,11 @@ int main(int , char* argv[]) {
c->SaveAs(filename+"_Ecs.png");

TF1* grE_SC = new TF1("S+Cfit","gaus",2.*low,2.*high); grE_SC->SetLineColor(kBlack);
TF1* grE_DR = new TF1("Efit","gaus",low,high); grE_DR->SetLineColor(kBlack);
tE_SC->SetOption("p"); tE_SC->Fit(grE_SC,"R+&same");
tE_DR->SetOption("p"); tE_DR->Fit(grE_DR,"R+&same");

tE_SC->Draw(""); c->SaveAs(filename+"_Esum.png");
tE_DR->Draw(""); c->SaveAs(filename+"_Ecorr.png");

TGraph* grSvsC = new TGraph(entries,&(E_Ss[0]),&(E_Cs[0]));
TGraph* grSvsC = new TGraph(E_Ss.size(),&(E_Ss[0]),&(E_Cs[0]));
grSvsC->SetTitle("SvsC;E_S;E_C");
grSvsC->SetMarkerSize(0.5); grSvsC->SetMarkerStyle(20);
grSvsC->GetXaxis()->SetLimits(0.,high);
Expand Down

0 comments on commit 73c571f

Please sign in to comment.