// ESAF : Euso Simulation and Analysis Framework
// $Id: EnergyModule.cc,v 1.74 2005/11/14 17:02:49 naumov Exp $
// Naumov Dmitry created Jul, 4 2004
//////////////////////////////////////////////////////////////////////////////
// //
// EnergyModule //
// //
// //
//////////////////////////////////////////////////////////////////////////////
//
/*

*/
//
//
/*

*/
//
#include "EnergyModule.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "Config.hh"
#include "EOpticsResponse.hh"
#include "EConst.hh"
#include "RecoRootEvent.hh"
#include "RecoShowerTrack.hh"
// ---> Begin of Energy Module helpers
#include "EnergyViewer.hh"
#include "EnergyFitUtils.hh"
// <--- End of Energy Module helpers
// ---> Begin of Root includes
#include <TCutG.h>
#include <TMinuit.h>
#include <TSpectrum.h>
#include <TProfile.h>
// <--- End of Root includes
#include <map>
#include <vector>
using namespace TMath;
using namespace EConst;
using namespace sou;
ClassImp(EnergyModule)
//_____________________________________________________________________________
EnergyModule::EnergyModule() :
RecoModule("Energy") {
// ctor
const string& ORdbname = Conf()->GetStr("
EnergyModule.fORname");
fOpticsResponse = new
EOpticsResponse(ORdbname.c_str());
Constructor();
}
//_____________________________________________________________________________
void EnergyModule::Constructor() {
//
// Default ctor initialization
fEvent = NULL;
fShowerTrack = NULL;
fGraphTheta = NULL;
fGraphPhi = NULL;
fGraphXY = NULL;
fGraphX = NULL;
fGraphY = NULL;
fHistoOverlap = NULL;
fHistoIdealOverlap = NULL;
fUseTrueAngles = kFALSE;
}
//_____________________________________________________________________________
EnergyModule::~
EnergyModule() {
// dtor
SafeDelete(
fShowerTrack);
ClearGraphs();
}
//_____________________________________________________________________________
Bool_t EnergyModule::Init() {
// Get fluorescence calculator
Msg(
EsafMsg::Info) << "Initializing " << MsgDispatch;
fDebugPs = (
Int_t)Conf()->GetNum("
EnergyModule.
fDebugPs");
fMinuitOutputLevel = (
Int_t)Conf()->GetNum("
EnergyModule.
fMinuitOutputLevel");
if (
fMinuitOutputLevel < -1 ||
fMinuitOutputLevel > 3 ) {
Msg(
EsafMsg::Warning) << "Wrong config value
EnergyModule.
fMinuitOutputLevel "
<<
fMinuitOutputLevel << " Setted to 1" << MsgDispatch;
fMinuitOutputLevel = 1;
}
fUseTrueAngles = (
Int_t)Conf()->GetNum("
EnergyModule.
fUseTrueAngles");
// Initialization of pixels viewer
if (
fDebugPs >= 1)
fEnergyViewer = new
EnergyViewer(this);
fShowerTrack = new
RecoShowerTrack;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t EnergyModule::PreProcess() {
fTheta = -1;
fPhi = -1;
fAltitude = -1;
fEnergyChi2 = -1;
fThetaFoVAtMax = -1;
fPhiFoVAtMax = -1;
fGtuTimeAtMax = -1;
fBackGround = -1;
fFCtime[0] =
fFCtime[1] =
fFCtime[2] =
fFCtime[3] = -1;
fPixelsCluster.clear();
fAllActivePixels.clear();
return kTRUE;
}
//_____________________________________________________________________________
void EnergyModule::CreateAuxilaryObjects() {
ClearGraphs();
fGraphTheta = new
TGraph();
fGraphPhi = new
TGraph();
fGraphXY = new
TGraph();
fGraphX = new
TGraph();
fGraphY = new
TGraph();
if (
fDebugPs >= 1)
fEnergyViewer->
InitFilePS();
}
//_____________________________________________________________________________
void EnergyModule::ClearGraphs() {
SafeDelete(
fGraphTheta);
SafeDelete(
fGraphPhi);
SafeDelete(
fGraphXY);
SafeDelete(
fGraphX);
SafeDelete(
fGraphY);
}
//_____________________________________________________________________________
Bool_t EnergyModule::LoadModules() {
vector <
Int_t> PixelsClusterInt;
// Load information from previous modules essential for the energy reconstruction
fTDM =
fEvent->
GetModuleData("TrackDirection");
// try to get TrackDirection2
if (
fTDM == NULL )
fTDM =
fEvent->
GetModuleData("TrackDirection2");
if (
fTDM == NULL ) Msg(
EsafMsg::Panic) << "None of
TrackDirectionModule and TrackDirectionModule2 is not found." << MsgDispatch;
if (
fUseTrueAngles) {
fTheta =
fEvent->
GetHeader().GetTrueTheta();
fPhi =
fEvent->
GetHeader().GetTruePhi();
}
else {
fTheta =
fTDM->
GetDouble("Theta");
fPhi =
fTDM->
GetDouble("Phi");
if(IsNaN(
fTheta) || IsNaN(
fPhi)) return kFALSE;
}
fShowerDir.
SetXYZ(Sin(
fTheta)*Cos(
fPhi),Sin(
fTheta)*Sin(
fPhi),Cos(
fTheta));
fGtuCM =
fEvent->
GetModuleData("GTUClustering");
if (
fGtuCM == NULL )
fGtuCM =
fEvent->
GetModuleData("HoughTransform");
if (
fGtuCM == NULL ) Msg(
EsafMsg::Panic) << "None of
GTUClusteringModule and
HoughModule is not found." << MsgDispatch;
if (
fGtuCM->
GetObj("CluPixels") == NULL )
Msg(
EsafMsg::Panic) << "Gtu Cluster Module has NULL CluPixels" << MsgDispatch;
PixelsClusterInt.clear();
PixelsClusterInt = *(vector<
Int_t>*)
fGtuCM->
GetObj("CluPixels");
fAllActivePixels = (Pixels_t)
fEvent->
GetRecoPixels();
for (
UInt_t i (0); i < PixelsClusterInt.size(); i++) {
RecoPixelData *pix =
fEvent->
GetRecoPixelData(PixelsClusterInt[i]);
fPixelsCluster.push_back(pix);
}
if (!
fPixelsCluster.size())
Msg(
EsafMsg::Panic) << "CluPixels has no pixels" << MsgDispatch;
Msg(
EsafMsg::Debug) << "Processing " <<
fPixelsCluster.size()
<< " points" << MsgDispatch;
fHmaxShapeM =
fEvent->
GetModuleData("HmaxByShapeMethod");
fHmaxProtonM =
fEvent->
GetModuleData("HmaxForProton");
if (
fHmaxShapeM == NULL &&
fHmaxProtonM == NULL)
fAltitude = 5*km;
if (
fHmaxShapeM)
fAltitude =
fHmaxShapeM->
GetDouble("Altitude");
if (
fHmaxProtonM)
fAltitude =
fHmaxProtonM->
GetDouble("Altitude");
return kTRUE;
}
//_____________________________________________________________________________
void EnergyModule::SortGtu() {
Int_t pmtId_previous = -1;
Float_t OffSet = 10.; // mm
map <
Int_t,
Int_t> EC;
PMTs_t::iterator itCell, itPmt;
Pixels_t::iterator itPix;
map <
Int_t,
Int_t>::iterator itEC;
fGtuMin =
fMyGtuMin = 10000;
fGtuMax =
fMyGtuMax = -10000;
fXmin =
fYmin =
fPhiMin =
fThetaMin = kHuge;
fXmax =
fYmax =
fPhiMax =
fThetaMax = -kHuge;
fGtuPixels.clear();
fGtuMyPixels.clear();
fMyPixelsCluster.clear();
fPixelsBgr.clear();
fECells.clear();
fPmts.clear(); EC.clear();
fActivePixels.clear();
// Sort Gtu and define cluster position
for (itPix =
fPixelsCluster.begin(); itPix !=
fPixelsCluster.end(); itPix++) {
Int_t pmtId =
fRunpars->
Pmt((*itPix)->GetPixelId());
Int_t gtu = (*itPix)->GetGtu();
Float_t theta = (*itPix)->GetTheta()*RadToDeg();
Float_t phi = (*itPix)->GetPhi()*RadToDeg();
fGtuPixels[gtu].push_back(*itPix);
if (
fGtuMin > gtu)
fGtuMin = gtu;
if (
fGtuMax < gtu)
fGtuMax = gtu;
if (
fPhiMin > phi)
fPhiMin = phi;
if (
fPhiMax < phi)
fPhiMax = phi;
if (
fThetaMin > theta)
fThetaMin = theta;
if (
fThetaMax < theta)
fThetaMax = theta;
if (pmtId == pmtId_previous) continue; // Avoid
double calculation
pmtId_previous = pmtId; // for PMT's
for(
Int_t iCorner(0); iCorner < 4; iCorner++){
TVector3 vPmt =
fRunpars->
PmtCorner(pmtId,iCorner);
if (
fXmin > vPmt.X())
fXmin = vPmt.X();
if (
fYmin > vPmt.Y())
fYmin = vPmt.Y();
if (
fXmax < vPmt.X())
fXmax = vPmt.X();
if (
fYmax < vPmt.Y())
fYmax = vPmt.Y();
}
EC[
fRunpars->
GetPmtGeo(pmtId)->GetEC()] = 0; // list of macrocells
}
// Make list of macrocells
for (itEC = EC.begin(); itEC != EC.end(); itEC++)
fECells.push_back(itEC->first);
fXmin -= OffSet;
fYmin -= OffSet;
fXmax += OffSet;
fYmax += OffSet;
// prepeare list of PMTs and pixels
for(
Int_t i = 1; i <=
fRunpars->
GetNumPmts(); i++) {
TVector3 vPmt =
fRunpars->
GetPmtGeo(i)->GetCenter();
if (vPmt.X() <
fXmin || vPmt.X() >
fXmax ||
vPmt.Y() <
fYmin || vPmt.Y() >
fYmax ) continue;
for (itCell =
fECells.begin(); itCell !=
fECells.end();itCell++) {
if (*itCell ==
fRunpars->
GetPmtGeo(i)->GetEC()) {
fPmts.push_back(i);
const
EPmtGeo* pmt =
fRunpars->
GetPmtGeo(i);
for (
Int_t PixUid = pmt->GetStartUniqueId(); PixUid <= pmt->GetLastUniqueId(); PixUid++)
fPixelsBgr.push_back(PixUid);
break;
}
}
}
fMyGtuMin =
fGtuMin;
fMyGtuMax =
fGtuMax;
fMyPixelsCluster =
fPixelsCluster;
fGtuMyPixels =
fGtuPixels;
// Make map of active pixels around cluster
for (itPix =
fAllActivePixels.begin(); itPix !=
fAllActivePixels.end(); itPix++)
if ((*itPix)->GetGtu() >=
fMyGtuMin && (*itPix)->GetGtu() <=
fMyGtuMax)
for (itPmt =
fPixelsBgr.begin(); itPmt !=
fPixelsBgr.end(); itPmt++)
if (*itPmt == (*itPix)->GetPixelId()) {
fActivePixels.push_back(*itPix);
break;
}
InitHisto(); // Histogram initialization
BuildFoVAngleFunctions(); // Build angles and fit distributions
BuildMyPixels(); // Define my pixels around signal
if (
fDebugPs >= 1)
fEnergyViewer->
DisplayFSTrack();
if (
fDebugPs >= 1)
fEnergyViewer->
DisplayPixels();
if (
fDebugPs >= 2)
DumpGtuSorted();
EC.clear();
}
//_____________________________________________________________________________
void EnergyModule::InitHisto() {
Float_t OffSet = 15;
Float_t OffSetDeg = 2;
if (
fPhiMax > 150. &&
fPhiMin < -150.) {
fHistoPhi = new
TH2F("GtuPhi","GtuPhi",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax,3600,-360,360);
}
else {
fHistoPhi = new
TH2F("GtuPhi","GtuPhi",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax,
5*(
Int_t)(
fPhiMax-
fPhiMin+2*OffSetDeg),
fPhiMin-OffSetDeg,
fPhiMax+OffSetDeg);
}
fHistoTheta = new
TH2F("GtuTheta","GtuTheta",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax,
5*(
Int_t)(
fThetaMax-
fThetaMin+2*OffSetDeg),
fThetaMin-OffSetDeg,
fThetaMax+OffSetDeg);
fHistoX = new
TH2F("GtuX","GtuX",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax,
5*(
Int_t)(
fXmax-
fXmin+2*OffSet),
fXmin-OffSet,
fXmax+OffSet);
fHistoY = new
TH2F("GtuY","GtuY",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax,
5*(
Int_t)(
fYmax-
fYmin+2*OffSet),
fYmin-OffSet,
fYmax+OffSet);
fHistoXY = new
TH2F("XY","XY",
5*(
Int_t)(
fXmax-
fXmin+2*OffSet),
fXmin-OffSet,
fXmax+OffSet,
5*(
Int_t)(
fYmax-
fYmin+2*OffSet),
fYmin-OffSet,
fYmax+OffSet);
fHistoOverlap = new
TProfile("Overlap","Overlap",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax);
fHistoIdealOverlap = new
TProfile("IdealOverlap","Overlap",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax);
}
//_____________________________________________________________________________
void EnergyModule::BuildFoVAngleFunctions() {
//
//
//
Float_t binx, biny, bintheta, binphi, phip;
Pixels_t OneGtuPixels;
Pixels_t::iterator itPix;
GtuPixels_t::iterator
Int_t pixel = 0;
itPix = (
fGtuMyPixels.begin()->second).begin();
phip = (*itPix)->GetPhi()*RadToDeg();
for (itGtu =
fGtuMyPixels.begin(); itGtu !=
fGtuMyPixels.end(); itGtu++) {
OneGtuPixels.clear();
OneGtuPixels = itGtu->second;
Int_t Count = 0;
binx = biny = bintheta = binphi = 0;
for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
Int_t nc = (*itPix)->GetCounts();
Float_t phi = (*itPix)->GetPhi()*RadToDeg();
if (phi-phip > 300.) {phi -= 360; phip = phi;}
else if (phi-phip < -300.) {phi += 360; phip = phi;}
fHistoPhi->Fill(itGtu->first,phi,nc);
fHistoTheta->Fill(itGtu->first,(*itPix)->GetTheta()*RadToDeg(),nc);
TVector3 PixelCenter=
fRunpars->
PixelCenter((*itPix)->GetPixelId());
fHistoX->Fill(itGtu->first,PixelCenter.X(),nc);
fHistoY->Fill(itGtu->first,PixelCenter.Y(),nc);
fHistoXY->Fill(PixelCenter.X(),PixelCenter.Y(),nc);
bintheta += (*itPix)->GetTheta()*RadToDeg()*nc;
binphi += phi*nc;
binx += PixelCenter.X()*nc;
biny += PixelCenter.Y()*nc;
Count += nc;
}
for (
Int_t icount = 0; icount < Count; icount++) {
fGraphTheta->
SetPoint(pixel,itGtu->first,bintheta/Count);
fGraphPhi->
SetPoint(pixel,itGtu->first,binphi/Count);
fGraphX->
SetPoint(pixel,itGtu->first,binx/Count);
fGraphY->
SetPoint(pixel,itGtu->first,biny/Count);
fGraphXY->
SetPoint(pixel,binx/Count,biny/Count);
pixel++;
}
}
fGraphTheta->
Fit("pol2","QEM");
fGraphPhi->
Fit("pol2","QEM");
fGraphX->
Fit("pol1","QEM");
fGraphY->
Fit("pol1","QEM");
fGraphXY->
Fit("pol1","QEM");
}
//_____________________________________________________________________________
void EnergyModule::BuildMyPixels() {
// Define the box around fitted line
PMTs_t::iterator itPmt;
Float_t Xnew[5], Ynew[5];
TVector2 VBox1, VBox2;
Float_t HistoSize = 15.*mm; // It is necessary to define it via Histo
TVector2 beginLine =
XYCluster(
fMyGtuMin);
TVector2 endLine =
XYCluster(
fMyGtuMax);
TVector2 Line = endLine-beginLine;
fMyPixels.clear();
VBox1.SetMagPhi(2.*HistoSize,Line.Phi());
// Coordinate in box corners
for (
Int_t Rot = 45, i = 0; Rot <= 315; Rot += 90, i++) {
VBox2 = VBox1.Rotate(Rot*DegToRad());
Xnew[i] = VBox2.X(); Ynew[i] = VBox2.Y();
if (!i || i == 3) {Xnew[i] += endLine.X(); Ynew[i] += endLine.Y();}
else {Xnew[i] += beginLine.X(); Ynew[i] += beginLine.Y();}
}
Xnew[4] = Xnew[0]; Ynew[4] = Ynew[0];
TCutG BoxCut("BoxCut",5,Xnew,Ynew);
// Define pixels within box
for (itPmt =
fPixelsBgr.begin(); itPmt !=
fPixelsBgr.end(); itPmt++) {
TVector3 VCenter =
fRunpars->
PixelCenter(*itPmt);
if (BoxCut.IsInside(VCenter.X(),VCenter.Y()))
fMyPixels.push_back(*itPmt);
}
}
//_____________________________________________________________________________
void EnergyModule::DumpGtuSorted() {
// dump on screen
Int_t MaxCounts = 0, GtuTimeAtMax = 0;
Pixels_t OneGtuPixels;
Pixels_t::iterator itPix;
GtuPixels_t::iterator
Msg(
EsafMsg::Debug) << "Dumped maps of pixels" << MsgDispatch;
for (itGtu =
fGtuPixels.begin(); itGtu !=
fGtuPixels.end(); itGtu++) {
OneGtuPixels.clear();
OneGtuPixels = itGtu->second;
Msg(
EsafMsg::Debug) << Form("GTU Entry: %d vector of %d pixels:",
itGtu->first,OneGtuPixels.size())<<MsgDispatch;
Int_t sum(0);
for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
Msg(
EsafMsg::Debug) << Form(" GTU %d PixelUid %d Counts %d",
(*itPix)->GetGtu(),(*itPix)->GetPixelId(),
(*itPix)->GetCounts()) << MsgDispatch;
sum += (*itPix)->GetCounts();
}
Msg(
EsafMsg::Debug) << Form("GTU Entry: Total counts %d",sum)<<MsgDispatch;
if (MaxCounts < sum) {
MaxCounts = sum;
GtuTimeAtMax = itGtu->first;
}
}
Msg(
EsafMsg::Debug) << Form("Most numerous GTU %d with %d count",
Int_t(GtuTimeAtMax),MaxCounts) << MsgDispatch;
Msg(
EsafMsg::Debug) << Form("Min and Max GTU %d %d ",
Int_t(
fGtuMin),
Int_t(
fGtuMax)) << MsgDispatch;
Msg(
EsafMsg::Debug) << Form("
DumpGtuSorted(): End Dumping") << MsgDispatch;
Msg(
EsafMsg::Debug) << Form("
DumpGtuSorted(): Dumping My map <GtuHit_t,Pixels_t>") << MsgDispatch;
MaxCounts = 0;
for (itGtu =
fGtuMyPixels.begin(); itGtu !=
fGtuMyPixels.end(); itGtu++) {
OneGtuPixels.clear();
OneGtuPixels = itGtu->second;
Msg(
EsafMsg::Debug) << Form("GTU Entry: %d vector of %d pixels:",
itGtu->first,OneGtuPixels.size())<<MsgDispatch;
Int_t sum(0);
for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
Msg(
EsafMsg::Debug) << Form(" GTU %d PixelUid %d Counts %d",
(*itPix)->GetGtu(), (*itPix)->GetPixelId(),
(*itPix)->GetCounts()) << MsgDispatch;
sum += (*itPix)->GetCounts();
}
Msg(
EsafMsg::Debug) << Form("GTU Entry: Total counts %d",sum)<<MsgDispatch;
if (MaxCounts < sum) {
MaxCounts = sum;
GtuTimeAtMax = itGtu->first;
}
}
Msg(
EsafMsg::Debug) << Form("Most numerous GTU %d with %d count",
Int_t(GtuTimeAtMax),MaxCounts) << MsgDispatch;
Msg(
EsafMsg::Debug) << Form("Min and Max GTU %d %d ",
Int_t(
fMyGtuMin),
Int_t(
fMyGtuMax)) << MsgDispatch;
Msg(
EsafMsg::Debug) << Form("
DumpGtuSorted(): End Dumping") << MsgDispatch;
}
//_____________________________________________________________________________
TVector2 EnergyModule::XYCluster(
Float_t time) {
//
if (!
fGraphX || !
fGraphY ) Msg(
EsafMsg::Panic) << "
fGraphX or
fGraphY does not exist! Check your code." << MsgDispatch;
Float_t XP0 =
fGraphX->
GetFunction("pol1")->GetParameter(0);
Float_t XP1 =
fGraphX->
GetFunction("pol1")->GetParameter(1);
Float_t YP0 =
fGraphY->
GetFunction("pol1")->GetParameter(0);
Float_t YP1 =
fGraphY->
GetFunction("pol1")->GetParameter(1);
TVector2 PosCluster(XP0 + XP1*time,YP0 + YP1*time);
return PosCluster;
}
//_____________________________________________________________________________
Float_t EnergyModule::ThetaFoV(
Float_t time) {
//
//
//
if (!
fGraphTheta) Msg(
EsafMsg::Panic) << "
fGraphTheta does not exist! Check your code." << MsgDispatch;
Float_t ThetaP0 =
fGraphTheta->
GetFunction("pol2")->GetParameter(0);
Float_t ThetaP1 =
fGraphTheta->
GetFunction("pol2")->GetParameter(1);
Float_t ThetaP2 =
fGraphTheta->
GetFunction("pol2")->GetParameter(2);
return DegToRad()*(ThetaP0 + ThetaP1*time + ThetaP2*time*time);
}
//_____________________________________________________________________________
Float_t EnergyModule::PhiFoV(
Float_t time) {
//
if (!
fGraphPhi) Msg(
EsafMsg::Panic) << "
fGraphPhi does not exist! Check your code." << MsgDispatch;
Float_t PhiP0 =
fGraphPhi->
GetFunction("pol2")->GetParameter(0);
Float_t PhiP1 =
fGraphPhi->
GetFunction("pol2")->GetParameter(1);
Float_t PhiP2 =
fGraphPhi->
GetFunction("pol2")->GetParameter(2);
return DegToRad()*(PhiP0 + PhiP1*time + PhiP2*time*time);
}
//_____________________________________________________________________________
Bool_t EnergyModule::BuildRecoveryFunction() {
// Build the recovery function along the track
// for a number of steps along the track compute overlap of efficacy with pixels coverage.
Float_t Xnew[5], Ynew[5];
PMTs_t::iterator itPixUid;
vector <
TCutG> Cuts;
vector <
TCutG>::const_iterator itCuts;
Int_t progress = 0;
Float_t MeanLambda = 360*nm;
Float_t GtuOffSet = 0., TimeStep = 0.25;
Float_t InitialTime =
fMyGtuMin - GtuOffSet + 0.5*TimeStep;
Float_t FinalTime =
fMyGtuMax + GtuOffSet - 0.4*TimeStep; // Not error
Float_t TotalTime =
fMyGtuMax -
fMyGtuMin + 2.*GtuOffSet - TimeStep;
// loop along the track
for (
Float_t time = InitialTime; time <= FinalTime; time += TimeStep) {
Float_t fThetaFoV =
ThetaFoV(time);
Float_t fPhiFoV =
PhiFoV(time);
TVector2 XYcluster =
XYCluster(time);
if (fThetaFoV > 30*DegToRad()) continue;
TH2F* htemp =
fOpticsResponse->
MakePsfHist(fThetaFoV,fPhiFoV,MeanLambda);
if (!htemp) continue;
Float_t xminPsf = htemp->GetXaxis()->GetXmin();
Float_t xmaxPsf = htemp->GetXaxis()->GetXmax();
Float_t yminPsf = htemp->GetYaxis()->GetXmin();
Float_t ymaxPsf = htemp->GetYaxis()->GetXmax();
Float_t deltaX = XYcluster.X() - 0.5*(xminPsf+xmaxPsf);
Float_t deltaY = XYcluster.Y() - 0.5*(yminPsf+ymaxPsf);
xminPsf += deltaX; xmaxPsf += deltaX;
yminPsf += deltaY; ymaxPsf += deltaY;
TH2F* Psf = new
TH2F("Psf","",htemp->GetNbinsX(),xminPsf,xmaxPsf,
htemp->GetNbinsY(),yminPsf,ymaxPsf);
for (
Int_t j = 1; j <= htemp->GetNbinsX(); j++)
for (
Int_t k = 1; k <= htemp->GetNbinsY(); k++) {
Float_t x = htemp->GetXaxis()->GetBinCenter(j) + deltaX;
Float_t y = htemp->GetYaxis()->GetBinCenter(k) + deltaY;
Psf->Fill(x,y,htemp->GetBinContent(j,k));
}
SafeDelete(htemp);
// Define the square for signal integration
Float_t PixelSide = 2.*mm;
Float_t histoX[5], histoY[5];
histoX[1] = histoX[2] = xmaxPsf + PixelSide;
histoY[2] = histoY[3] = yminPsf - PixelSide;
histoX[0] = histoX[3] = histoX[4] = xminPsf - PixelSide;
histoY[0] = histoY[1] = histoY[4] = ymaxPsf + PixelSide;
TCutG HistoCut("HistoCut",5,histoX,histoY);
// Define the pixels inside histogram
Cuts.clear();
for (itPixUid =
fMyPixels.begin(); itPixUid !=
fMyPixels.end(); itPixUid++) {
TVector3 VCenter =
fRunpars->
PixelCenter(*itPixUid);
if (HistoCut.IsInside(VCenter.X(),VCenter.Y())) {
for(
Int_t iCorner = 0; iCorner < 4; iCorner++) {
TVector3 dummy =
fRunpars->
PixelCorner(*itPixUid,iCorner);
Xnew[iCorner] = dummy.X(); Ynew[iCorner] = dummy.Y();
}
Xnew[4] = Xnew[0]; Ynew[4] = Ynew[0];
TCutG PixelCut("PixelCut",5,Xnew,Ynew);
Cuts.push_back(PixelCut);
}
}
Float_t IntegratedEfficacy(0);
itCuts = Cuts.begin();
for (
Int_t i = 1; i <= Psf->GetNbinsX(); i++) {
Float_t x = Psf->GetXaxis()->GetBinCenter(i);
for (
Int_t j = 1; j <= Psf->GetNbinsY(); j++) {
Float_t y = Psf->GetYaxis()->GetBinCenter(j);
if (itCuts != Cuts.end() && (*itCuts).IsInside(x,y))
IntegratedEfficacy += Psf->GetBinContent(i,j);
else
for (itCuts = Cuts.begin(); itCuts != Cuts.end(); itCuts++)
if ((*itCuts).IsInside(x,y)) {
IntegratedEfficacy += Psf->GetBinContent(i,j);
break;
}
}
}
if (
fDebugPs >= 2)
fEnergyViewer->
DisplayRecoveryDebug(Psf, IntegratedEfficacy, time+0.5*TimeStep);
fHistoOverlap->
Fill(time,IntegratedEfficacy);
fHistoIdealOverlap->
Fill(time,Psf->Integral());
if (100.*(time-InitialTime)/TotalTime >= progress) {
Msg(
EsafMsg::Info).SetProgress(progress);
Msg(
EsafMsg::Info) << "
BuildRecoveryFunction():" << MsgCount;
progress+=10;
}
SafeDelete(Psf);
}
if (progress <= 100) Msg(
EsafMsg::Info) << MsgDispatch;
if (
fDebugPs >= 1)
fEnergyViewer->
DisplayRecovery();
return kTRUE;
}
//_____________________________________________________________________________
void EnergyModule::FitTimeDistribution() {
// Fit temporal distribution
// define the data point and plot them for the check
Int_t CurrentGtuCounts, CurrentGtuCountsUnfolded;
Pixels_t OneGtuPixels;
Pixels_t::iterator itPix;
GtuPixels_t::iterator
// Currently is not used
// for (itGtu = fGtuPixels.begin(); itGtu != fGtuPixels.end(); itGtu++) {
// CurrentGtuCounts = 0;
// OneGtuPixels.clear();
// OneGtuPixels = itGtu->second;
// for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
// CurrentGtuCounts += (*itPix)->GetCounts();
// if (fBackground)
// CurrentGtuCounts-= fRunpars->GetNightGlowRateByUId((*itPix)->GetPixelId())*
// fEvent->GetHeader().GetDetector()->GetElectronics()->GetGtuLength()/microsecond;
// }
// point++;
// }
// make my cluster
TH1F* htime = new
TH1F("htime","Time distribution",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax);
TH1F* htimeS = new
TH1F("htimeS","Unsaturated time distribution",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax);
TH1F* htimeRec = new
TH1F("Recovery Time distribution","",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax);
for (itGtu =
fGtuMyPixels.begin(); itGtu !=
fGtuMyPixels.end(); itGtu++) {
CurrentGtuCounts = 0;
CurrentGtuCountsUnfolded = 0;
OneGtuPixels.clear();
OneGtuPixels = itGtu->second;
for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end();itPix++) {
CurrentGtuCounts += (*itPix)->GetCounts();
CurrentGtuCountsUnfolded += UnfoldSaturation((*itPix)->GetCounts());
// if (fBackground)
// CurrentGtuCounts-=fRunpars->GetNightGlowRateByUId((*itPix)->GetPixelId()) *
// fEvent->GetHeader().GetDetector()->GetElectronics()->GetGtuLength()/microsecond;
}
htime->Fill(itGtu->first,CurrentGtuCounts);
htimeS->Fill(itGtu->first,CurrentGtuCountsUnfolded);
}
// Find Chrenkov signal
FindCherenkovPeak(htimeS);
htimeRec->Divide(htimeS,
fHistoOverlap);
FindCherenkovPeak(htimeRec);
// Generate Shower Track, fill expected histo for flux of photons on EP
BuildShowerTrack();
// Fitting histogram
TF1 *func_energy = new
TF1("func_energy",Fit_Energy,
fMyGtuMin,
fMyGtuMax,2);
TF1 *func_time = new
TF1("func_time",TimeDistr,
fMyGtuMin,
fMyGtuMax,2);
TF1 *func_bg = new
TF1("func_bg",Background,
fMyGtuMin,
fMyGtuMax,0);
gMinuit->
SetObjectFit(this);
func_energy->SetParLimits(0,0,100);
func_energy->SetParameter(0,5.);
func_energy->SetParameter(1,22.);
htimeS->Fit("func_energy","QEMN");
Double_t *par = func_energy->GetParameters();
Double_t *err = func_energy->GetParErrors();
fEnergy = par[0]*1.e+14;
fEnergyError = err[0]*1.e+14;
fEnergyChi2 = func_energy->GetChisquare()/(
fMyGtuMax-
fMyGtuMin-2.);
if (
fDebugPs >= 1) {
htimeRec->SetMarkerColor(kBlue);
htimeRec->SetMarkerStyle(22);
htime->SetMarkerColor(kBlack);
htime->SetMarkerStyle(21);
htimeS->SetMarkerColor(kBlack);
htimeS->SetMarkerStyle(21);
func_time->SetLineStyle(1);
func_time->SetLineColor(1);
func_time->SetParameters(par);
func_bg->SetLineColor(kBlue);
func_energy->SetLineStyle(2);
func_energy->SetLineColor(2);
fEnergyViewer->
DisplayTimeFit(htimeRec,func_time,func_bg);
fEnergyViewer->
DisplayTimeFit(htime,func_energy,func_bg);
fEnergyViewer->
DisplayTimeFit(htimeS,func_energy,func_bg);
}
Msg(
EsafMsg::Info) << Form("Simulated Reconstructed Information")
<< MsgDispatch;
Msg(
EsafMsg::Info) << Form("Theta %.2E %.2E deg",
fEvent->
GetHeader().GetTrueTheta()*RadToDeg(),
fTheta*RadToDeg()) << MsgDispatch;
Msg(
EsafMsg::Info) << Form("Phi %.2E %.2E deg",
fEvent->
GetHeader().GetTruePhi()*RadToDeg(),
fPhi*RadToDeg()) << MsgDispatch;
Msg(
EsafMsg::Info) << Form("Energy %.2E %.2E MeV",
fEvent->
GetHeader().GetTrueEnergy(),
fEnergy) << MsgDispatch;
SafeDelete(htime);
SafeDelete(htimeRec);
SafeDelete(htimeS);
SafeDelete(func_energy);
SafeDelete(func_time);
SafeDelete(func_bg);
}
//_____________________________________________________________________________
void EnergyModule::FindCherenkovPeak(
TH1F *htime) {
Int_t ShowerTime, TimePeak;
Int_t MaxSignalGtu = 0, npeak = 0, ChBin = 0, ChBinMin = 0;
Int_t ShowerTime2 =
fMyGtuMax -
fMyGtuMin;
Int_t ShowerTime1 = 0;
Float_t *PeakPosition;
TSpectrum *spectrum = new
TSpectrum(10,1);
Int_t PeaksNumber = spectrum->Search(htime,2,"");
if (PeaksNumber) {
// All time unit defines via bin number
// Determine time caracteristics of event
PeakPosition = spectrum->GetPositionX();
// Find a reasonable time of shower evolution
// Searching after last peak
TimePeak = htime->FindBin(PeakPosition[PeaksNumber-1]);
for (
Int_t iTime = TimePeak+1; iTime <
fMyGtuMax; iTime++) {
if (!(htime->GetBinContent(iTime))) {
ShowerTime2 = iTime;
break;
}
}
// Searching before first peak
TimePeak = htime->FindBin(PeakPosition[0]);
for (
Int_t iTime = TimePeak; iTime >
fMyGtuMin; iTime--) {
if (!(htime->GetBinContent(iTime))) {
ShowerTime1 = iTime;
break;
}
}
ShowerTime = ShowerTime2 - ShowerTime1; // Time of shower
// Find maximum signal
for (
Int_t iPeak = 0; iPeak < PeaksNumber; iPeak++) {
Int_t Time = htime->FindBin(PeakPosition[iPeak]);
if (htime->GetBinContent(MaxSignalGtu) < htime->GetBinContent(Time)) {
MaxSignalGtu = Time;
npeak = iPeak;
}
}
ChBin = htime->FindBin(PeakPosition[PeaksNumber-1]);
// Check if the maximum signal is Cherenkov light
if (npeak == PeaksNumber-1) {
if (1.*(MaxSignalGtu-ShowerTime1)/ShowerTime > 0.8) {
if (PeaksNumber > 1) {
// Cherenkov peak founded, try to find fluorescent maximum
MaxSignalGtu = 0;
for (
Int_t iPeak = 0; iPeak < PeaksNumber-1; iPeak++) {
Int_t Time = htime->FindBin(PeakPosition[iPeak]);
if (htime->GetBinContent(MaxSignalGtu) < htime->GetBinContent(Time))
MaxSignalGtu = Time;
}
}
else {MaxSignalGtu = (
Int_t)(0.5*(ShowerTime1+ShowerTime2));}
}
else {ChBin = 0;}
}
// Find a Cherenkov light
if (ChBin) {
// Position of maximum signal from Cherenkov light
for (
Int_t iTime = ChBin-5; iTime < ShowerTime2; iTime++)
if (htime->GetBinContent(ChBin) < htime->GetBinContent(iTime))
ChBin = iTime;
// Time of Cherenkov light signal on focal plane
ChBinMin = ChBin;
for (
Int_t iTime = ChBin; iTime > ChBin - 6; iTime--)
if (htime->GetBinContent(ChBinMin) > htime->GetBinContent(iTime))
ChBinMin = iTime;
}
}
else {
// Put the peak to the point with maximum density
Float_t Sumt = 0;
for (
Int_t iTime =
fMyGtuMin+2; iTime <
fMyGtuMax-3; iTime++) {
Float_t Sum = 0;
for (
Int_t iDh = iTime - 2; iDh <= iTime + 2; iDh ++)
Sum += htime->GetBinContent(iDh);
if (Sum > Sumt) {
MaxSignalGtu = iTime;
Sumt = Sum;
}
}
}
// Theta and Phi at shower maximum
fPhiFoVAtMax =
PhiFoV(MaxSignalGtu);
fThetaFoVAtMax =
ThetaFoV(MaxSignalGtu);
fGtuTimeAtMax = MaxSignalGtu;
Msg(
EsafMsg::Debug) << Form("FindThetaPhiAtMaximum(): Gtu %d, Theta %3.2f deg, Phi %3.2f deg", MaxSignalGtu,
fThetaFoVAtMax*RadToDeg(),
fPhiFoVAtMax*RadToDeg()) << MsgDispatch;
if (
fDebugPs >= 1) {
// Draw polymarker for Cherenkov region and shower maximum
Int_t nPoints = 3;
Float_t x[6], y[6];
x[0] =
fMyGtuMin + ShowerTime1;
y[0] = 0;
x[1] =
fMyGtuMin + MaxSignalGtu-0.5;
y[1] = htime->GetBinContent(MaxSignalGtu);
x[2] =
fMyGtuMin + ShowerTime2;
y[2] = 0;
if (ChBin) {
x[3] =
fMyGtuMin + ChBin-0.5;
y[3] = htime->GetBinContent(ChBin);
nPoints = 4;
if (ChBinMin) {
x[4] =
fMyGtuMin + ChBinMin-0.5;
y[4] = htime->GetBinContent(ChBinMin);
x[5] =
fMyGtuMin + 2*ChBin - ChBinMin - 0.5;
y[5] = htime->GetBinContent(2*ChBin - ChBinMin);
nPoints = 6;
}
}
fEnergyViewer->
DisplayCherenkov(htime,nPoints,x,y);
}
SafeDelete(spectrum);
fFCtime[0] =
fMyGtuMin + MaxSignalGtu-0.5;
fFCtime[1] =
fMyGtuMin + ChBinMin-0.5;
fFCtime[2] =
fMyGtuMin + ChBin-0.5;
fFCtime[3] =
fMyGtuMin + 2*ChBin - ChBinMin - 0.5;
}
//_____________________________________________________________________________
void EnergyModule::EstimateBackground() {
PMTs_t BgrPix;
map <
Int_t,
Float_t> simBgr;
map <
Int_t,
Int_t> Npix, simNpix;
map <
Int_t,
TH1F*> hBgrCount;
PMTs_t::iterator itPmt, itMyPmt;
Pixels_t::iterator itPix;
Npix.clear(); simBgr.clear(); simNpix.clear();
BgrPix.clear(); hBgrCount.clear();
// Make a list of pixels not in signal
for (itPmt =
fPixelsBgr.begin(); itPmt !=
fPixelsBgr.end(); itPmt++) {
for (itMyPmt =
fMyPixels.begin(); itMyPmt !=
fMyPixels.end(); itMyPmt++)
if (*itPmt == *itMyPmt) break;
BgrPix.push_back(*itPmt);
}
// create histogram
for (itPmt =
fECells.begin(); itPmt !=
fECells.end(); itPmt++)
hBgrCount[*itPmt] = new
TH1F(Form("Count%d",*itPmt),"Gtu",
fMyGtuMax-
fMyGtuMin,
fMyGtuMin,
fMyGtuMax);
// fill background histo
for (itPix =
fActivePixels.begin();itPix !=
fActivePixels.end(); itPix++) {
for (itPmt = BgrPix.begin(); itPmt != BgrPix.end(); itPmt++)
if (*itPmt == (*itPix)->GetPixelId()) {
Int_t ECuid =
fRunpars->
ElementaryCell((*itPix)->GetPixelId());
hBgrCount[ECuid]->Fill((*itPix)->GetGtu(),(*itPix)->GetCounts());
break;
}
}
// Calculate number of pixels in each elementary cells
for (itPmt = BgrPix.begin(); itPmt != BgrPix.end(); itPmt++)
Npix[
fRunpars->
ElementaryCell(*itPmt)]++;
// Calculate simulated background
for (itPmt =
fPixelsBgr.begin(); itPmt !=
fPixelsBgr.end(); itPmt++) {
simBgr[
fRunpars->
ElementaryCell(*itPmt)] +=
fRunpars->
GetNightGlowRateByUId(*itPmt);
simNpix[
fRunpars->
ElementaryCell(*itPmt)]++;
}
// calculate time interval excluding zero bin for each elementary cell
Int_t nec = 0;
Float_t rec = 0, sim = 0;
for (itPmt =
fECells.begin(); itPmt !=
fECells.end(); itPmt++) {
Int_t Ngtu = 0;
Float_t integral = 0;
for (
Int_t i = 1; i <= hBgrCount[*itPmt]->GetNbinsX(); i++)
if (hBgrCount[*itPmt]->GetBinContent(i)) {
integral += hBgrCount[*itPmt]->GetBinContent(i);
Ngtu++;
}
// cout << "Cell " << *itPmt << " rec: "
// << hBgrCount[*itPmt]->Integral()/(1.*Npix[*itPmt]*Ngtu)
// << " sim: " << 2.5*NsimBgr[*itPmt]/NpixBgr[*itPmt] << endl;
rec += integral/(1.*Npix[*itPmt]*Ngtu);
sim += 2.5*simBgr[*itPmt]/simNpix[*itPmt];
nec++;
// hBgrCount[*itPmt]->Scale(1./Npix[*itPmt]);
// fEnergyViewer->DrawHisto(hBgrPoisson[*itPmt]);
// fEnergyViewer->DrawHisto(hBgrCount[*itPmt]);
SafeDelete(hBgrCount[*itPmt]);
}
hBgrCount.clear(); BgrPix.clear(); simNpix.clear(); simBgr.clear();
// Int_t np = 0, bgr = 0;
// for (itPix = fPixelsCluster.begin();itPix != fPixelsCluster.end(); itPix++) {
// bgr += (*itPix)->GetCountsNoise(); np++;
// }
// cout << "Mean value sim: " << 1.*sim/nec << endl;
// cout << "Simulated background: " << 1.*bgr/np << endl;
// cout << "Reconstructed: " << 1.*rec/nec << endl;
fBackGround = 1.*rec/nec;
}
//_____________________________________________________________________________
Float_t EnergyModule::BgrPerPixelPerGtu(
Float_t time) {
//
//
//
return
fBackGround;
}
//_____________________________________________________________________________
void EnergyModule::BuildShowerTrack() {
//
//
//
fShowerTrack->
SetEnergyModule(this);
fShowerTrack->
SetHmax(
fAltitude);
fShowerTrack->
SetThetaPhiFoV(
fThetaFoVAtMax,
fPhiFoVAtMax);
fShowerTrack->
SetThetaPhi(
fTheta,
fPhi);
fShowerTrack->
MakeTrack();
}
//_____________________________________________________________________________
Bool_t EnergyModule::Process(
RecoEvent *ev) {
// manager of the energy reconstruction chain
fEvent = ev;
fRunpars =
fEvent->
GetHeader().GetRunPars();
gErrorIgnoreLevel = 1;
CreateAuxilaryObjects();
// Load information from previous modules
if (!
LoadModules()) return kFALSE;
// Make map <GtuHit_t,Pixels_t>
SortGtu();
// Estimate Background
EstimateBackground();
// Build the recovery function along the track
BuildRecoveryFunction();
// Fit the time distribution
FitTimeDistribution();
// Save data
MyData()->Add("Energy",
fEnergy);
if (
fDebugPs >= 1)
fEnergyViewer->
DisplayPixels3D("ClusterSurf");
if (
fDebugPs >= 1)
fEnergyViewer->
DisplayPixels3D("FocalSurf");
if (
fDebugPs >= 1)
fEnergyViewer->
CloseFilePS();
if (
fDebugPs >= 1)
fEnergyViewer->
Clean();
gErrorIgnoreLevel = 0;
fShowerTrack->
Clear();
return kTRUE;
}
//_____________________________________________________________________________
Bool_t EnergyModule::PostProcess() {
return kTRUE;
}
//_____________________________________________________________________________
Bool_t EnergyModule::SaveRootData(
RecoRootEvent *fRecoRootEvent) {
//
// Save desired information into the reco root file
//
fRecoRootEvent->GetRecoEnergy().SetEnergy(
fEnergy,
fEnergyError);
fRecoRootEvent->GetRecoEnergy().SetQuality(
fEnergyChi2);
Double_t MCError = fRecoRootEvent->GetTruth().GetEnergy() -
fEnergy;
fRecoRootEvent->GetRecoEnergy().SetEnergyMCError(MCError);
fRecoRootEvent->GetRecoEnergy().SetThetaFoV(
fThetaFoVAtMax);
fRecoRootEvent->GetRecoEnergy().SetPhiFoV(
fPhiFoVAtMax);
return kTRUE;
}
//_____________________________________________________________________________
Bool_t EnergyModule::Done() {
//
// Here we inform that the module is ended its work
//
if (
fDebugPs >= 1) SafeDelete(
fEnergyViewer);
Msg(
EsafMsg::Info) << "Completed" << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
void EnergyModule::UserMemoryClean() {
//
// Here we clear any auxilary objects created during one event energy reconstruction
//
fPixelsCluster.clear();
fMyPixelsCluster.clear();
fAllActivePixels.clear();
fActivePixels.clear();
fShowerDir.Clear();
fGtuPixels.clear();
fGtuMyPixels.clear();
fMyPixels.clear();
fPixelsBgr.clear();
fECells.clear();
fPmts.clear();
SafeDelete(
fHistoTheta);
SafeDelete(
fHistoPhi);
SafeDelete(
fHistoX);
SafeDelete(
fHistoY);
SafeDelete(
fHistoXY);
SafeDelete(
fHistoOverlap);
SafeDelete(
fHistoIdealOverlap);
}