diff --git a/analysis/analysis.cpp b/analysis/analysis.cpp
index 4ddb451..f4f9087 100644
--- a/analysis/analysis.cpp
+++ b/analysis/analysis.cpp
@@ -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"
@@ -21,25 +21,31 @@
 #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.);
@@ -47,44 +53,59 @@ int main(int , char* argv[]) {
   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();
@@ -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);
       }
     }
 
@@ -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);
@@ -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);