// ESAF : Euso Simulation and Analysis Framework
// $Id: EnergyModule.cc,v 1.24 2005/04/19 14:27:31 naumov Exp $
// Naumov Dmitry created Jul, 4 2004
//////////////////////////////////////////////////////////////////////////////
// //
// EnergyModule //
// //
// //
//////////////////////////////////////////////////////////////////////////////
#include "EnergyModule.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "FluoCalculator.hh"
#include "Config.hh"
#include "EsafSpectrum.hh"
#include "LightSourceFactory.hh"
#include "SinglePhoton.hh"
#include "RadiativeProcessesCalculator.hh"
#include "O1_ClearSkyPropagator.hh"
#include "RadiativeFactory.hh"
#include "EOpticsResponse.hh"
#include "EORSample.hh"
#include "EORHeader.hh"
#include "ERunParameters.hh"
#include "EDetector.hh"
#include "EConst.hh"
#include "ShowerTrack.hh"
#include "RecoRootEvent.hh"
#include "SlastShowerGenerator.hh"
#include <TH2F.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TView.h>
#include <TPolyLine.h>
#include <TStyle.h>
#include <TGraphErrors.h>
#include <TText.h>
#include <TCutG.h>
#include <TSpline.h>
#include <TMinuit.h>
#include <TFumili.h>
#include <TVirtualFitter.h>
#include <TPaveText.h>
#include <TLegend.h>
#include <map>
#include <vector>
using namespace TMath;
using namespace EConst;
ClassImp(EnergyModule)
ClassImp(EnergyModuleData)
//______________________________________________________________________________
Double_t PhotoElectronsTotal(Double_t *x, Double_t *par) {
// time (in GTU units)
// par[0] is fitted energy
// par[1] is fitted time at the shower maximum
// par[2] is fitted time sigma
EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
if (!fData) return 0;
Float_t time = x[0];
Float_t energy = par[0];
Float_t timemax = par[1];
Float_t sigma = par[2];
Float_t fQuantumEff = fData->fQuantumEff;
Float_t fAttenuation = fData->fAttenuation;
Float_t fYield = fData->fYield;
Float_t fShowerLength = fData->fShowerLength;
Float_t fRmax = fData->fRmax;
TSpline3 *fPsfIdeal = (TSpline3*)fData->fPsfIdeal;
// make the energy to scale around 10**20 eV
Float_t Theory = energy*1.e14/(1450)*fQuantumEff*fAttenuation*fYield*fShowerLength*fPsfIdeal->Eval(time);
Theory = Theory/( 4.e0*Pi()*fRmax*fRmax);
Theory *= Exp(-0.5*Power(time-timemax,2)/Power(sigma,2));
return Theory;
}
//______________________________________________________________________________
Double_t PhotoElectronsRecovery(Double_t *x, Double_t *par) {
// time (in GTU units)
// par[0] is fitted energy
// par[1] is fitted time at the shower maximum
// par[2] is fitted time sigma
EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
if (!fData) return 0;
Float_t time = x[0];
Float_t energy = par[0];
Float_t timemax = par[1];
Float_t sigma = par[2];
Float_t fQuantumEff = fData->fQuantumEff;
Float_t fAttenuation = fData->fAttenuation;
Float_t fYield = fData->fYield;
Float_t fShowerLength = fData->fShowerLength;
Float_t fRmax = fData->fRmax;
TSpline3 *fRecovery = (TSpline3*)fData->fRecovery;
// make the energy to scale around 10**20 eV
Float_t Theory = energy*1.e14/(1450)*fQuantumEff*fAttenuation*fYield*fShowerLength*fRecovery->Eval(time);
Theory = Theory/( 4.e0*Pi()*fRmax*fRmax);
Theory *= Exp(-0.5*Power(time-timemax,2)/Power(sigma,2));
return Theory;
}
//______________________________________________________________________________
Double_t PhotoElectrons2Total(Double_t *x, Double_t *par) {
// time (in GTU units)
// par[0] is fitted energy
// par[1] is fitted time at the shower maximum
EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
if (!fData) return 0;
Float_t time = x[0];
Float_t energy = par[0];
Float_t timeshift = par[1];
Float_t fQuantumEff = fData->fQuantumEff;
TSpline3 *fPsfIdeal = (TSpline3*)fData->fPsfIdeal;
TH1F *fFlux = (TH1F*)fData->fFlux;
Int_t bin = fFlux->GetXaxis()->FindBin(time-timeshift);
Float_t Theory = energy*fPsfIdeal->Eval(time)*fQuantumEff;
Theory *= fFlux->GetBinContent(bin);
return Theory;
}
//______________________________________________________________________________
Double_t PhotoElectrons2Recovery(Double_t *x, Double_t *par) {
// time (in GTU units)
// par[0] is fitted energy
// par[1] is fitted time at the shower maximum
EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
if (!fData) return 0;
Float_t time = x[0];
Float_t energy = par[0];
Float_t timeshift = par[1];
Float_t fQuantumEff = fData->fQuantumEff;
TSpline3 *fRecovery = (TSpline3*)fData->fRecovery;
TH1F *fFlux = (TH1F*)fData->fFlux;
Int_t bin = fFlux->GetXaxis()->FindBin(time-timeshift);
Float_t Theory = energy*fRecovery->Eval(time)*fQuantumEff;
Theory *= fFlux->GetBinContent(bin);
return Theory;
}
//______________________________________________________________________________
void PhotoElectronsTimeFit(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
// time (in GTU units)
// par[0] is fitted energy
// par[1] is fitted time at the shower maximum
// par[2] is fitted time sigma
EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
if (!fData) return;
Float_t energy = par[0];
Float_t timemax = par[1];
Float_t sigma = par[2];
Float_t fQuantumEff = fData->fQuantumEff;
Float_t fAttenuation = fData->fAttenuation;
Float_t fYield = fData->fYield;
Float_t fShowerLength = fData->fShowerLength;
Float_t fRmax = fData->fRmax;
TGraphErrors *fPhotoElectrons = (TGraphErrors*)fData->fMyPhotoElectrons;
TSpline3 *fRecovery = (TSpline3*)fData->fRecovery;
Double_t Chi2(0);
for (Int_t i(0); i < fPhotoElectrons->GetN(); i++) {
Double_t time,pe;
fPhotoElectrons->GetPoint(i,time,pe);
// find efficacy integrated over pixels
Float_t Theory = energy*1.e14/1450*fRecovery->Eval(time)*fQuantumEff*fAttenuation*fYield*fShowerLength;
Theory = Theory/( 4.e0*Pi()*fRmax*fRmax);
Theory *= Exp(-0.5*Power(time-timemax,2)/Power(sigma,2));
Chi2 += Power(pe-Theory,2)/pe;
}
f = (Double_t)Chi2/(fPhotoElectrons->GetN() - 3);
}
//______________________________________________________________________________
void PhotoElectronsTimeFit2(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
// time (in GTU units)
// par[0] is fitted energy
// par[1] is fitted time at the shower maximum
EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
if (!fData) return;
Float_t energy = par[0];
Float_t timeshift = par[1];
Float_t fQuantumEff = fData->fQuantumEff;
TGraphErrors *fPhotoElectrons = (TGraphErrors*)fData->fMyPhotoElectrons;
TSpline3 *fRecovery = (TSpline3*)fData->fRecovery;
TH1F *fFlux = (TH1F*)fData->fFlux;
Double_t Chi2(0);
for (Int_t i(0); i < fPhotoElectrons->GetN(); i++) {
Double_t time,pe;
fPhotoElectrons->GetPoint(i,time,pe);
Int_t bin = fFlux->GetXaxis()->FindBin(time-timeshift);
Float_t Theory = energy*fRecovery->Eval(time)*fQuantumEff;
Theory *= fFlux->GetBinContent(bin);
Chi2 += Power(pe-Theory,2)/pe;
}
f = (Double_t)Chi2/(fPhotoElectrons->GetN() - 2);
}
//_____________________________________________________________________________
void EnergyModuleData::Clear(){
fQuantumEff = 0;
fAttenuation = 0;
fYield = 0;
fShowerLength = 0;
fRmax = 0;
SafeDelete(fRecovery);
SafeDelete(fPsfIdeal);
SafeDelete(fFlux);
SafeDelete(fPhotoElectrons);
SafeDelete(fMyPhotoElectrons);
}
//_____________________________________________________________________________
EnergyModule::EnergyModule() : RecoModule("Energy") {
// ctor
const string& ORdbname = Conf()->GetStr("EnergyModule.fORname");
fOpticsResponse = new EOpticsResponse(ORdbname.c_str());
fEsafSpectrum = new EsafSpectrum();
fShowerGenerator = new SlastShowerGenerator;
fEPspectrum = NULL;
fEvent = NULL;
fEPphotons = NULL;
fShowerTrack = NULL;
fCanvas = NULL;
fPsfOverlap = NULL;
fPsfIdealOverlap = NULL;
fRecovery = NULL;
fPhotoElectrons = NULL;
fMyPhotoElectrons = NULL;
fBgGraph = NULL;
//
}
//_____________________________________________________________________________
EnergyModule::~EnergyModule() {
// dtor
SafeDelete(fEsafSpectrum);
SafeDelete(fEPspectrum);
SafeDelete(fEPphotons);
SafeDelete(fOpticsResponse);
SafeDelete(fShowerTrack);
SafeDelete(fShowerGenerator);
SafeDelete(fCanvas);
ClearGraphs();
}
//_____________________________________________________________________________
Bool_t EnergyModule::Init() {
// Get fluorescence calculator
Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
string fluoname = Conf()->GetStr("EnergyModule.fFluoCalculator");
fFluoCalculator = LightSourceFactory::Get()->GetFluoCalculator( fluoname );
fRadCalculator = RadiativeFactory::Get()->GetRadiativeProcessesCalculator();
fEusoAlt = Conf()->GetNum("EnergyModule.fEusoAltitude")*km;
fEusoRadius = Conf()->GetNum("EnergyModule.fEusoRadius")*mm;
fSaveEps = Conf()->GetStr("EnergyModule.fSaveEps");
fEnergyDistributionName = Conf()->GetStr("EnergyModule.fEnergyDistributionName");
fViewThreshold = (Int_t)Conf()->GetNum("EnergyModule.fViewThreshold");
fThreshold = (Int_t)Conf()->GetNum("EnergyModule.fThreshold");
fMinuitOutputLevel = (Int_t)Conf()->GetNum("EnergyModule.fMinuitOutputLevel");
fFitMode = (Int_t)Conf()->GetNum("EnergyModule.fFitMode");
if ( fMinuitOutputLevel < -1 || fMinuitOutputLevel > 3 ) {
Msg(EsafMsg::Warning) << "Wrong config value EnergyModule.fMinuitOutputLevel "
<< fMinuitOutputLevel << " Setted to 1" << MsgDispatch;
fMinuitOutputLevel = 1;
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t EnergyModule::PreProcess() {
fAltitude = -1;
fFluorMaxPe = -1;
fDeltaL = -1;
fEnergyChi2 = -1;
fGtuMin = 1.e21;
fGtuMax = -1.21;
fGtuAtMax = -1;
fGtuTimeAtMax = -1;
fXmax = -1.e20;
fXmin = 1.e20;
fYmax = -1.e20;
fYmin = 1.e20;
fData.Clear();
fGtuPixels.clear();
fGtuMyPixels.clear();
fPixelLines.clear();
fPmtLines.clear();
SafeDelete(fEPspectrum);
return kTRUE;
}
//_____________________________________________________________________________
void EnergyModule::CreateAuxilaryObjects() {
ClearGraphs();
fCanvas = new TCanvas("EnergyModule");
fGraphTheta = new TGraphErrors();
fGraphPhi = new TGraphErrors();
fGraphXY = new TGraphErrors();
fPsfOverlap = new TGraph();
fPsfIdealOverlap = new TGraph();
fBgGraph = new TGraph();
fRecovery = new TSpline3();
fPhotoElectrons = new TGraphErrors();
fMyPhotoElectrons = new TGraphErrors();
fEPspectrum = new TH1F("EPspectrum","EP",100,300*nm,400*nm);
}
//_____________________________________________________________________________
void EnergyModule::ClearGraphs() {
SafeDelete(fCanvas);
SafeDelete(fGraphTheta);
SafeDelete(fGraphPhi);
SafeDelete(fGraphXY);
SafeDelete(fPsfOverlap);
SafeDelete(fPsfIdealOverlap);
SafeDelete(fRecovery);
SafeDelete(fPhotoElectrons);
SafeDelete(fMyPhotoElectrons);
SafeDelete(fBgGraph);
}
//_____________________________________________________________________________
void EnergyModule::LoadModules() {
// Load information from previous modules essential for the energy reconstruction
fTDM = fEvent->GetModuleData("TrackDirection");
if ( fTDM == NULL ) Msg(EsafMsg::Panic) << "TrackDirectionModule is not found." << MsgDispatch;
fTheta = fTDM->GetDouble("Theta");
fPhi = fTDM->GetDouble("Phi");
fShowerDir.SetXYZ(Sin(fTheta)*Cos(fPhi), Sin(fTheta)*Sin(fPhi), Cos(fTheta));
fGtuCM = fEvent->GetModuleData("GTUClustering");
if ( fGtuCM == NULL ) Msg(EsafMsg::Panic) << "GTUClusteringModule is not found." << MsgDispatch;
if ( fGtuCM->GetObj("CluPixels") == NULL )
Msg(EsafMsg::Panic) << "Gtu Cluster Module has NULL CluPixels" << MsgDispatch;
fPixelsCluster.clear();
fAllActivePixels.clear();
fPixelsCluster = *(vector<Int_t>*) fGtuCM->GetObj("CluPixels");
fAllActivePixels = (vector<RecoPixelData*>)fEvent->GetRecoPixels();
if ( (fNumPoints = fPixelsCluster.size()) == 0) Msg(EsafMsg::Panic) << "CluPixels has no pixels"
<< MsgDispatch;
Msg(EsafMsg::Debug) << "Processing " << fNumPoints << " 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");
}
//_____________________________________________________________________________
Float_t EnergyModule::MakeFluorYield() {
// Estimate fluorescence yield for the shower maximum
TF2* EnergyDistribution = (TF2*)ShowerTrack::GetEnergyDistribution(fEnergyDistributionName);
TF12 EnergyDistributionMax("EDS_name",EnergyDistribution,1.e0,"X");
fTotalYield = fFluoCalculator->GetFluoYield(fAltitude,&EnergyDistributionMax,fEsafSpectrum);
// Estimate the fluorescence spectrum at the Entrace Pupil
Int_t NumTrials(0);
Float_t H1 = fEusoAlt-fAltitude;
Float_t ct2 = Power(Cos(fThetaFoVAtMax),2);
Float_t st2 = Power(Sin(fThetaFoVAtMax),2);
Float_t tmp1 = fEusoAlt*st2 - EarthRadius()*ct2;
Float_t tmp2 = Power(fEusoAlt,2)*st2 - (2*EarthRadius()*fAltitude + fAltitude*fAltitude);
Float_t z = tmp1 + Sqrt(tmp1*tmp1 - tmp2);
Float_t rho = H1*Tan(fThetaFoVAtMax);
Float_t x = rho*Cos(fPhiFoVAtMax);
Float_t y = rho*Sin(fPhiFoVAtMax);
EarthVector MaximumPosition(x,y,z);
fMaxPos = MaximumPosition;
EarthVector ISS(0,0,fEusoAlt);
EarthVector dir = (ISS - MaximumPosition).Unit();
fRmax = (ISS - MaximumPosition).Mag();
fDeltaL = EConst::C()*fEvent->GetHeader().GetGtuLength()/(1 + dir*fShowerDir);
fMeanAttenuation = 0;
fMeanLambdaOnEntrancePupil = 0;
// O1_ClearSkyPropagator::GoToDetector(EsafSpectrum& spec,const EarthVector& pos)
while (kTRUE) {
// FIXME An optimization can be done if the whole spectrum at once can be transported.
// FIXME This is what Lowtran does but this is not fully here
NumTrials++;
Float_t lambda = fEsafSpectrum->GetLambda();
SinglePhoton s(0, lambda, MaximumPosition, dir, (PhotonType)1, (PhotonStatus)1);
fMeanAttenuation += fRadCalculator->Trans(s,ISS);
fMeanLambdaOnEntrancePupil += lambda;
fEPspectrum->Fill(lambda);
if (NumTrials > 100) break;
}
fMeanAttenuation /= NumTrials;
fMeanLambdaOnEntrancePupil /= NumTrials;
Msg(EsafMsg::Debug) << Form("MakeFluorYield(): Begin Dump") << MsgDispatch;
Msg(EsafMsg::Debug) << Form(" Estimated ShowerMax(x,y,z): (%2.3f, %2.3f, %2.3f) km",x/km,y/km,z/km)
<< MsgDispatch;
Msg(EsafMsg::Debug) << Form(" Estimated Dir To Euso (Theta,Phi): (%2.3f, %2.3f) deg",
dir.Theta()*RadToDeg(),dir.Phi()*RadToDeg()) << MsgDispatch;
Msg(EsafMsg::Debug) << Form(" Estimated Mean Attenuation: (%2.3f)",fMeanAttenuation) << MsgDispatch;
Msg(EsafMsg::Debug) << Form(" Estimated Mean Lambda on EP: (%2.3f) nm",
fMeanLambdaOnEntrancePupil/nm) << MsgDispatch;
Msg(EsafMsg::Debug) << Form(" Estimated Mean Yield: (%2.3f) ph/m",fTotalYield*m)
<< MsgDispatch;
Msg(EsafMsg::Debug) << Form("MakeFluorYield(): End Dump") << MsgDispatch;
return fTotalYield;
}
//_____________________________________________________________________________
void EnergyModule::DefinePixelsView() {
// returnes vector of RecoPixelData which should be considered in both Viewer and Energy reconstruction
PMTs_t pmts;
//PMTs_t::iterator itPmt;
//Pixels_t::iterator itPixel;
// make list of PMTs which containes found pixels
// find box with found cluster
for (UInt_t i (0); i < fPixelsCluster.size(); i++) {
RecoPixelData *pix = fEvent->GetRecoPixelData(fPixelsCluster[i]);
Int_t pmtId = fRunpars->Pmt(pix->GetPixelId());
for(Int_t iCorner(0); iCorner<4; iCorner++){
TVector3 dummy = fRunpars->PmtCorner(pmtId,iCorner);
if (fXmin > dummy.X()) fXmin = dummy.X();
if (fYmin > dummy.Y()) fYmin = dummy.Y();
if (fXmax < dummy.X()) fXmax = dummy.X();
if (fYmax < dummy.Y()) fYmax = dummy.Y();
}
}
// increase size of the box by OffSet
Float_t OffSet = 10; // mm
fXmin -= OffSet;
fYmin -= OffSet;
fXmax += OffSet;
fYmax += OffSet;
// find pmts in this box
for( Int_t i(1); i <= fRunpars->GetNumPmts(); i++) {
TVector3 dummy = fRunpars->GetPmtGeo(i)->GetCenter();
if ( dummy.X() < fXmin || dummy.X() > fXmax || dummy.Y() < fYmin || dummy.Y() > fYmax ) continue;
// check if this pmt is not yet there
UInt_t pmtsize = pmts.size();
Bool_t IsPmtThere = kFALSE;
if (pmtsize) {
for (UInt_t j(0); j<pmtsize; j++) {
if (pmts[i] == i) {
IsPmtThere = kTRUE;
break;
}
}
}
if (!IsPmtThere) pmts.push_back(i);
}
// make list of lines to draw arond each PMT
fPmtLines.clear();
fPixelLines.clear();
for (UInt_t i(0); i < pmts.size(); i++) {
Float_t x[5], y[5];
Int_t pmtUid = pmts[i];
for(Int_t iCorner(0); iCorner<4; iCorner++){
TVector3 dummy = fRunpars->PmtCorner(pmtUid,iCorner);
x[iCorner] = dummy.X();
y[iCorner] = dummy.Y();
if (fXmin > x[iCorner]) fXmin = x[iCorner];
if (fYmin > y[iCorner]) fYmin = y[iCorner];
if (fXmax < x[iCorner]) fXmax = x[iCorner];
if (fYmax < y[iCorner]) fYmax = y[iCorner];
}
x[4] = x[0];
y[4] = y[0];
TPolyLine *p = new TPolyLine(5,x,y);
p->SetLineColor(20);
fPmtLines[pmtUid] = p;
// for each pmt to draw prepare the list of pixels inside the given pmt
const EPmtGeo* pmt = fRunpars->GetPmtGeo(pmtUid);
for (Int_t PixUid = pmt->GetStartUniqueId(); PixUid <= pmt->GetLastUniqueId(); PixUid++) {
Float_t x[5], y[5];
for(Int_t iCorner(0); iCorner<4; iCorner++){
TVector3 dummy = fRunpars->PixelCorner(PixUid,iCorner);
x[iCorner] = dummy.X();
y[iCorner] = dummy.Y();
}
x[4] = x[0];
y[4] = y[0];
TPolyLine *p = new TPolyLine(5,x,y);
p->SetLineColor(kWhite);
Color_t color = 20;
Style_t style = 3490; // special style for defaut
// search if this pixel belongs to active pixels
for (UInt_t j(0); j<fAllActivePixels.size(); j++){
if (PixUid == fAllActivePixels[j]->GetPixelId()) {
// make yellow color for active pixel with number of counts above the ViewThreshold
if (fAllActivePixels[j]->GetCounts() > fViewThreshold) color = 20;
// make grey color for active pixel with number of counts below the ViewThreshold
else color = 20;
// define fill style
for (UInt_t k(0); k<fPixelsCluster.size(); k++){
RecoPixelData* pixel = fEvent->GetRecoPixelData(fPixelsCluster[k]);
if (pixel->GetPixelId() == PixUid) {
style = 1001; // solid fill for the cluster
break;
}
}
break;
}
else
// make black color for not active pixel
color = 1;
}
p->SetFillColor(color);
fPixelLines[PixUid] = p;
}
}
}
//_____________________________________________________________________________
void EnergyModule::DisplayCluster() {
DefinePixelsView();
if (fCanvas) fCanvas->Clear();
if (gPad) {
Float_t OffSet = 5;
gPad->Range(fXmin-OffSet,fYmin-OffSet,fXmax+OffSet,fYmax+OffSet);
}
// Draw PMTs
map <Int_t, TPolyLine* >::const_iterator itPmt;
for (itPmt = fPmtLines.begin(); itPmt != fPmtLines.end(); itPmt++) {
itPmt->second->Draw("f");
itPmt->second->Draw();
}
// compute number of p.e. integrated over each pixel and draw the text in the center of each pixel.
map <Int_t,Int_t> mapPixelsIntegral;
for (UInt_t i(0) ; i<fAllActivePixels.size(); i++) {
RecoPixelData *pix = fAllActivePixels[i];
mapPixelsIntegral[pix->GetPixelId()] += pix->GetCounts();
}
// Draw Pixels
for (itPmt = fPixelLines.begin(); itPmt != fPixelLines.end(); itPmt++) {
itPmt->second->Draw("f");
itPmt->second->Draw();
Int_t counts = mapPixelsIntegral[itPmt->first];
TVector3 PixelCenter = fRunpars->PixelCenter(itPmt->first);
Float_t x = PixelCenter.X() + 0.2*fRunpars->GetPmtData().GetPadSide();
Float_t y = PixelCenter.Y() + 0.2*fRunpars->GetPmtData().GetPadSide();
TText *text = new TText(x,y,Form("%d",counts));
Float_t padMin = Min( gPad->GetWh(),gPad->GetWw() );
Float_t textsize = 3.e-5*padMin;
text->SetTextSize(textsize);
text->Draw();
}
// Draw fitted line track
if (fGraphXY)
fGraphXY->Draw("PX");
if (fSaveEps == "yes")
if(fCanvas) fCanvas->Print(Form("output/check_focalplane.run%d.event%d.eps",
fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
}
//_____________________________________________________________________________
void EnergyModule::SortGtu() {
// Prepare map <GtuHit_t,Pixels_t>
fGtuPixels.clear();
Pixels_t OneGtuPixels;
for (UInt_t i(0); i < fPixelsCluster.size(); i++) {
RecoPixelData *pix = fEvent->GetRecoPixelData(fPixelsCluster[i]);
if (fGtuMin > pix->GetGtu()) fGtuMin = pix->GetGtu();
if (fGtuMax < pix->GetGtu()) fGtuMax = pix->GetGtu();
}
Int_t previousGtu = fEvent->GetRecoPixelData(fPixelsCluster[0])->GetGtu();
for (UInt_t i(0); i < fPixelsCluster.size(); i++) {
RecoPixelData *pix = fEvent->GetRecoPixelData(fPixelsCluster[i]);
if (pix->GetGtu() == previousGtu) {
OneGtuPixels.push_back(i);
}
else {
fGtuPixels[previousGtu] = OneGtuPixels;
previousGtu = pix->GetGtu();
OneGtuPixels.clear();
OneGtuPixels.push_back(i);
}
}
// make another map of pixels from all active pixels
fGtuMyPixels.clear();
previousGtu = fEvent->GetRecoPixelData(fPixelsCluster[0])->GetGtu();
for (UInt_t i(0); i < fAllActivePixels.size(); i++) {
RecoPixelData *pix = fAllActivePixels[i];
if (pix->GetGtu() == previousGtu) {
OneGtuPixels.push_back(i);
}
else {
fGtuMyPixels[previousGtu] = OneGtuPixels;
previousGtu = pix->GetGtu();
OneGtuPixels.clear();
OneGtuPixels.push_back(i);
}
}
DumpGtuSorted();
}
//_____________________________________________________________________________
void EnergyModule::DumpGtuSorted(){
// dump on screen
Msg(EsafMsg::Info) << Form("DumpGtuSorted(): Dumping maps <GtuHit_t,Pixels_t> into Output File")
<< MsgDispatch;
Msg(EsafMsg::Debug) << Form("DumpGtuSorted(): Dumping maps <GtuHit_t,Pixels_t>") << MsgDispatch;
Msg(EsafMsg::Debug) << Form("DumpGtuSorted(): Dumping TDM map <GtuHit_t,Pixels_t>") << MsgDispatch;
GtuPixels_t::iterator it;
Pixels_t OneGtuPixels;
Int_t MaxCounts(0);
for (it = fGtuPixels.begin(); it != fGtuPixels.end(); it++) {
OneGtuPixels.clear();
OneGtuPixels = it->second;
Msg(EsafMsg::Debug) << Form("GTU Entry: %d vector of %d pixels:",it->first,OneGtuPixels.size())
<< MsgDispatch;
Int_t sum(0);
for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[OneGtuPixels[i]]);
Msg(EsafMsg::Debug)
<< Form(" GTU %d PixelUid %d Counts %d",pix->GetGtu(),pix->GetPixelId(),pix->GetCounts())
<< MsgDispatch;
sum += pix->GetCounts();
}
Msg(EsafMsg::Debug) << Form("GTU Entry: Total counts %d",sum) << MsgDispatch;
if (MaxCounts < sum) {
MaxCounts = sum;
fGtuTimeAtMax = it->first;
}
}
Msg(EsafMsg::Debug) << Form("Most numerous GTU %d with %d count",Int_t(fGtuTimeAtMax),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 (it = fGtuMyPixels.begin(); it != fGtuMyPixels.end(); it++) {
OneGtuPixels.clear();
OneGtuPixels = it->second;
Msg(EsafMsg::Debug) << Form("GTU Entry: %d vector of %d pixels:",it->first,OneGtuPixels.size())
<< MsgDispatch;
Int_t sum(0);
for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
RecoPixelData* pix = fAllActivePixels[OneGtuPixels[i]];
Msg(EsafMsg::Debug)
<< Form(" GTU %d PixelUid %d Counts %d",pix->GetGtu(),pix->GetPixelId(),pix->GetCounts())
<< MsgDispatch;
sum += pix->GetCounts();
}
Msg(EsafMsg::Debug) << Form("GTU Entry: Total counts %d",sum) << MsgDispatch;
if (MaxCounts < sum) {
MaxCounts = sum;
fGtuTimeAtMax = it->first;
}
}
Msg(EsafMsg::Debug) << Form("Most numerous GTU %d with %d count",Int_t(fGtuTimeAtMax),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;
}
//_____________________________________________________________________________
void EnergyModule::BuildFoVAngleFunctions() {
//
// make two maps
Int_t pixels(0);
GtuPixels_t::iterator it;
Pixels_t OneGtuPixels;
for (it = fGtuPixels.begin(); it != fGtuPixels.end(); it++) {
OneGtuPixels.clear();
OneGtuPixels = it->second;
for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[OneGtuPixels[i]]);
fGraphTheta->SetPoint(pixels,it->first,pix->GetTheta()*RadToDeg());
fGraphTheta->SetPointError(pixels,0,Sqrt(1.e0/pix->GetCounts()));
fGraphPhi->SetPoint(pixels,it->first,pix->GetPhi()*RadToDeg());
fGraphPhi->SetPointError(pixels,0,Sqrt(1.e0/pix->GetCounts()));
TVector3 PixelCenter = fRunpars->PixelCenter(pix->GetPixelId());
fGraphXY->SetPoint(pixels,PixelCenter.X(), PixelCenter.Y());
fGraphXY->SetPointError(pixels,0,Sqrt(1.e0/pix->GetCounts()));
pixels++;
}
}
fCanvas->Clear();
fCanvas->Divide(1,2);
fCanvas->cd(1);
fGraphTheta->Fit("pol1","Q");
fGraphTheta->Draw("A*X");
fGraphTheta->GetYaxis()->SetTitle("#Theta, [deg]");
fCanvas->cd(2);
fGraphPhi->Fit("pol1","Q");
fGraphPhi->Draw("A*X");
fGraphPhi->GetYaxis()->SetTitle("#Phi, [deg]");
fGraphPhi->GetXaxis()->SetTitle("Gtu Hit");
fGraphXY->Fit("pol1","Q");
fGraphXY->SetMarkerColor(kRed);
fGraphXY->SetMarkerStyle(kMultiply);
fGraphXY->SetMarkerSize(1.e-3*Min( gPad->GetWh(),gPad->GetWw() ));
fGraphXY->GetFunction("pol1")->SetLineStyle(2);
fGraphXY->GetFunction("pol1")->SetLineWidth(2);
fGraphXY->GetFunction("pol1")->SetLineColor(kBlue);
if (fSaveEps == "yes")
if(fCanvas) fCanvas->Print(Form("output/check_fits.run%d.event%d.eps",
fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
}
//_____________________________________________________________________________
void EnergyModule::FindThetaPhiAtMaximum() {
Int_t MaxCounts(0);
GtuPixels_t::iterator it;
Pixels_t OneGtuPixels;
for (it = fGtuMyPixels.begin(); it != fGtuMyPixels.end(); it++) {
OneGtuPixels.clear();
OneGtuPixels = it->second;
Int_t CurrentGtuCounts(0);
Float_t CurrentGtuTheta(0);
Float_t CurrentGtuPhi(0);
for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
// RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[OneGtuPixels[i]]);
RecoPixelData* pix = fAllActivePixels[OneGtuPixels[i]];
CurrentGtuCounts += pix->GetCounts();
CurrentGtuTheta += pix->GetTheta()*pix->GetCounts();
CurrentGtuPhi += pix->GetPhi()*pix->GetCounts();
}
if (MaxCounts < CurrentGtuCounts) {
MaxCounts = CurrentGtuCounts;
fThetaFoVAtMax = CurrentGtuTheta/CurrentGtuCounts;
fPhiFoVAtMax = CurrentGtuPhi/CurrentGtuCounts;
}
}
Msg(EsafMsg::Debug) << Form("FindThetaPhiAtMaximum(): Theta %3.2f deg, Phi %3.2f deg",
fThetaFoVAtMax*RadToDeg(),fPhiFoVAtMax*RadToDeg()) << MsgDispatch;
}
//_____________________________________________________________________________
Float_t EnergyModule::ThetaFoV(Float_t time) {
//
if (!fGraphTheta) return 0;
Float_t fThetaFoVP0 = fGraphTheta->GetFunction("pol1")->GetParameter(0);
Float_t fThetaFoVP1 = fGraphTheta->GetFunction("pol1")->GetParameter(1);
return DegToRad()*(fThetaFoVP0 + time*fThetaFoVP1);
}
//_____________________________________________________________________________
Float_t EnergyModule::PhiFoV(Float_t time) {
//
if (!fGraphPhi) return 0;
// Get fit parameters to set up the functions
Float_t fPhiFoVP0 = fGraphPhi->GetFunction("pol1")->GetParameter(0);
Float_t fPhiFoVP1 = fGraphPhi->GetFunction("pol1")->GetParameter(1);
return DegToRad()*(fPhiFoVP0 + time*fPhiFoVP1);
}
//_____________________________________________________________________________
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.
Int_t nSteps = 4*Int_t(fGtuMax - fGtuMin);
Int_t CurrentStep = 0;
Int_t GtuOffSet = 2; // extend the track time to (GtuMin - DeltaGtu*GtuOffSet, GtuMax + DeltaGtu*GtuOffSet)
Float_t TimeStep = (fGtuMax - fGtuMin + 2*GtuOffSet*fEvent->GetHeader().GetGtuLength()/microsecond)/nSteps;
Float_t time = fGtuMin - GtuOffSet;
map <Int_t, TPolyLine* >::const_iterator itPixels;
Msg(EsafMsg::Info) << "BuildRecoveryFunction(): 0%" << MsgFlush;
Bool_t next(0),subnext(0);
Float_t left(0),right(10),subleft(0),subright(2.5);
while (CurrentStep<=nSteps) {
Float_t fThetaFoV = ThetaFoV(time);
Float_t fPhiFoV = PhiFoV(time);
// TH2F* Psf = fOpticsResponse->MakePsfHist(fThetaFoV,fPhiFoV,fEPspectrum);
TH2F* Psf = fOpticsResponse->MakePsfHist(fThetaFoV,fPhiFoV,fMeanLambdaOnEntrancePupil);
Float_t IntegratedEfficacy(0);
for (itPixels = fPixelLines.begin(); itPixels!=fPixelLines.end(); itPixels++) {
/* Int_t entry=-1;
for (UInt_t j(0); j<fPixelsCluster.size(); j++){
RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[j]);
if (itPixels->first == pix->GetPixelId()) {
entry = j;
break;
}
}
if (entry == -1) continue;
RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[entry]);
if (pix->GetCounts() < fThreshold) continue;*/
TPolyLine *p = itPixels->second;
Double_t *x = p->GetX();
Double_t *y = p->GetY();
TCutG *PixelCut = new TCutG("PixelCut",p->GetN(),p->GetX(),p->GetY());
// find min and max of x and y to save time needed for integration
Float_t xmin(1.e20), xmax(-1.e20), ymin(1.e20), ymax(-1.e20);
for (Int_t i(0); i < p->GetN(); i++) {
if (xmin > x[i]) xmin = x[i];
if (xmax < x[i]) xmax = x[i];
if (ymin > y[i]) ymin = y[i];
if (ymax < y[i]) ymax = y[i];
}
Int_t xbinFirst = Psf->GetXaxis()->FindBin(xmin);
Int_t xbinLast = Psf->GetXaxis()->FindBin(xmax);
Int_t ybinFirst = Psf->GetYaxis()->FindBin(ymin);
Int_t ybinLast = Psf->GetYaxis()->FindBin(ymax);
for (Int_t i = xbinFirst; i <= xbinLast; i++) {
Float_t x = Psf->GetXaxis()->GetBinCenter(i);
for (Int_t j = ybinFirst; j <=ybinLast; j++ ) {
Float_t y = Psf->GetYaxis()->GetBinCenter(j);
if (PixelCut->IsInside(x,y)) {
IntegratedEfficacy += Psf->GetBinContent(i,j);
}
}
}
}
// end of integration
fPsfOverlap->SetPoint(CurrentStep,time,IntegratedEfficacy);
fPsfIdealOverlap->SetPoint(CurrentStep,time,Psf->Integral());
SafeDelete(Psf);
time += TimeStep;
CurrentStep++;
Float_t fraction = 1.e2*CurrentStep/nSteps;
if (left<fraction && fraction <=right)
next = kFALSE;
else if (fraction>right) {
next = kTRUE;
left = right;
right += 10;
}
if (subleft<fraction && fraction <=subright) {
subnext = kFALSE;
}
if (fraction > subright) {
subnext = kTRUE;
subleft = subright;
subright +=2;
}
if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;
}
Msg(EsafMsg::Info) << " -done. " << MsgDispatch;
if(fCanvas) fCanvas->Clear();
fPsfOverlap->Draw("A*X");
fPsfOverlap->GetYaxis()->SetTitle("Psf Overlap with pixels, [mm#^2]");
fPsfOverlap->GetXaxis()->SetTitle("Gtu Hit");
SafeDelete(fRecovery);
SafeDelete(fEffIdeal);
fRecovery = new TSpline3("Recovery",fPsfOverlap->GetX(),fPsfOverlap->GetY(),fPsfOverlap->GetN());
fRecovery->Draw("same");
fEffIdeal = new TSpline3("Recovery",fPsfIdealOverlap->GetX(),fPsfIdealOverlap->GetY(),fPsfIdealOverlap->GetN());
fEffIdeal->SetLineColor(kRed);
fEffIdeal->Draw("same");
if (fSaveEps == "yes")
if(fCanvas) fCanvas->Print(Form("output/check_recovery.run%d.event%d.eps",
fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
return kTRUE;
}
//_____________________________________________________________________________
void EnergyModule::FillDataStructure() {
// Fill data container object fData of type EnergyModuleData
fData.fQuantumEff = fQuantumEff;
fData.fAttenuation = fMeanAttenuation;
fData.fYield = fTotalYield;
fData.fRmax = fRmax;
fData.fShowerLength = fDeltaL;
fData.fRecovery = (TSpline3*)fRecovery->Clone();
fData.fPhotoElectrons = (TGraphErrors*)fPhotoElectrons->Clone();
fData.fMyPhotoElectrons = (TGraphErrors*)fMyPhotoElectrons->Clone();
fData.fPsfIdeal = (TSpline3*)fEffIdeal->Clone();
fData.fFlux = (TH1F*)fEPphotons->Clone();
}
//_____________________________________________________________________________
Bool_t EnergyModule::ShowerTrack3D() {
//
// make the first point and get the shower track for energy of 1.e21 and reconstructed angles
Float_t Length = 30*km;
EarthVector FirstPoint = fMaxPos + fShowerDir*Length;
fShowerGenerator->SetShower(1.e21,fTheta,fPhi,FirstPoint);
fShowerGenerator->SetQuiet(kTRUE);
fShowerGenerator->SetDepthStep(0.5);
fShowerTrack = (ShowerTrack*)fShowerGenerator->Get();
// make estimated time (in GTU) distribution of photons on Entrance Pupil
EarthVector ISS(0,0,fEusoAlt);
Float_t tMin(1.e20), tMax(-1.e20);
Double_t tempData[fShowerTrack->Size()][2];
Double_t FluxMax(-1), tFluxMax(-1);
for (UInt_t i(0); i < fShowerTrack->Size(); i++) {
const ShowerStep& s = fShowerTrack->GetStep(i);
// get position and compute for it: fluo yield + attenuation + time at EUSO
EarthVector pos = 0.5*(s.GetXYZi() + s.GetXYZf());
Float_t step_length = (s.GetXYZi() - s.GetXYZf()).Mag();
Float_t time = 0.5*(s.GetTimei() + s.GetTimef());
EarthVector dir = (ISS - pos);
Double_t SolidAngle = 1/dir.Mag2()/(4*Pi());
Double_t flux = SolidAngle*s.GetNelectrons()*fMeanAttenuation*fTotalYield*step_length;
Float_t EusoTime = time + dir.Mag()/EConst::C();
tempData[i][0] = EusoTime;
tempData[i][1] = flux;
if (tMin > EusoTime) tMin = EusoTime;
if (tMax < EusoTime) tMax = EusoTime;
if (FluxMax < flux) {
FluxMax = flux;
tFluxMax = EusoTime;
}
}
SafeDelete(fEPphotons);
// shift time distribution to have maximum close to that by p.e. distribution
Float_t fGtuLength = fEvent->GetHeader().GetDetector()->GetElectronics().GetGtuLength();
tMax -= (tFluxMax - fGtuLength*fGtuTimeAtMax);
tMin -= (tFluxMax - fGtuLength*fGtuTimeAtMax);
Float_t bins = (tMax - tMin)/fGtuLength;
fEPphotons = new TH1F("EPphotons","",Int_t(bins),tMin/fGtuLength,tMax/fGtuLength);
for (UInt_t i(0); i < fShowerTrack->Size(); i++) {
Float_t GtuTime = (tempData[i][0]-(tFluxMax - fGtuLength*fGtuTimeAtMax))/fGtuLength;
Float_t photons = tempData[i][1];
fEPphotons->Fill(GtuTime,photons);
}
if (fSaveEps == "yes") {
if(fCanvas) {
fCanvas->Clear();
fEPphotons->Draw();
fCanvas->Print(Form("output/check_shower_photons.run%d.event%d.eps",
fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
}
}
fShowerTrack->Reset();
fShowerGenerator->Reset();
return kTRUE;
}
//_____________________________________________________________________________
void EnergyModule::FitTimeDistribution() {
// Fit temporal distribution
// define the data point and plot them for the check
GtuPixels_t::iterator it;
Pixels_t OneGtuPixels;
Int_t point(0);
for (it = fGtuPixels.begin(); it != fGtuPixels.end(); it++) {
OneGtuPixels.clear();
OneGtuPixels = it->second;
Int_t CurrentGtuPixels(0);
Float_t CurrentGtuCounts(0);
for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[OneGtuPixels[i]]);
if (pix->GetCounts() >= fThreshold) {
CurrentGtuCounts += pix->GetCounts();
CurrentGtuPixels++;
}
}
// substruct the background
// CurrentGtuCounts -= fBackground*CurrentGtuPixels;
fPhotoElectrons->SetPoint(point,it->first,CurrentGtuCounts);
fPhotoElectrons->SetPointError(point,0,Sqrt(1.e0/CurrentGtuCounts));
point++;
}
// make my cluster
point = 0;
for (it = fGtuMyPixels.begin(); it != fGtuMyPixels.end(); it++) {
OneGtuPixels.clear();
OneGtuPixels = it->second;
Int_t CurrentGtuPixels(0);
Float_t CurrentGtuCounts(0);
for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
RecoPixelData* pix = fAllActivePixels[OneGtuPixels[i]];
if (pix->GetCounts() >= fThreshold) {
CurrentGtuCounts += pix->GetCounts();
CurrentGtuPixels++;
}
}
// substruct the background
// CurrentGtuCounts -= fBackground*CurrentGtuPixels;
fMyPhotoElectrons->SetPoint(point,it->first,CurrentGtuCounts);
fMyPhotoElectrons->SetPointError(point,0,Sqrt(1.e0/CurrentGtuCounts));
point++;
}
if(fCanvas) fCanvas->Clear();
fMyPhotoElectrons->SetMarkerColor(kBlue);
fMyPhotoElectrons->SetMarkerStyle(22);
fPhotoElectrons->SetMarkerStyle(21);
fPhotoElectrons->SetMarkerColor(kBlack);
fMyPhotoElectrons->Draw("APE");
fPhotoElectrons->Draw("PE");
fPhotoElectrons->GetYaxis()->SetTitle("p.e.");
fPhotoElectrons->GetXaxis()->SetTitle("Gtu Hit");
TLegend *legend=new TLegend(0.1,0.85,0.3,1);
legend->SetTextFont(72);
legend->SetTextSize(0.02);
legend->AddEntry(fMyPhotoElectrons,"all pixels","lp");
legend->AddEntry(fPhotoElectrons,"TDM cluster","lp");
legend->Draw();
// make the fit
FillDataStructure();
if (fFitMode == 1) FitMode1();
if (fFitMode == 2) FitMode2();
if (fFitMode == 3) {FitMode1(); FitMode2();}
if (fFitMode == 4) FitMode2fumili();
}
//_____________________________________________________________________________
void EnergyModule::FitMode1() {
TMinuit *minuit = new TMinuit(3);
minuit->SetPrintLevel(fMinuitOutputLevel);
minuit->SetObjectFit(&fData);
minuit->SetFCN(PhotoElectronsTimeFit);
minuit->SetErrorDef(1.);
Double_t par[3] = {1.e0,10,10};
Double_t step[3] = {1.e-1,1.e-1,1.e-1};
Double_t min[3] = {0.00, 0.0, 0.0};
Double_t max[3] = {0.00, 1.e3, 1.e3};
string cpar[3] = {"Energy","tMax","Sigma"};
for( Int_t i=0; i<3; i++ ) {
minuit->DefineParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
}
minuit->Migrad();
// retrieve results
Double_t outpar[3],err[3];
for( Int_t i=0; i<3; i++ ) {
minuit->GetParameter(i,outpar[i],err[i]);
}
fEnergy = outpar[0]*1e14;
fEnergyError = err[0]*1e14;
Double_t amin,edm,errdef;
Int_t nvpar,nparx,icstat;
minuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
fEnergyChi2 = amin/nvpar;
TF1 *f1 = new TF1("f1",PhotoElectronsTotal,fEvent->GetStartGtu(),fEvent->GetEndGtu(),3);
TF1 *f2 = new TF1("f2",PhotoElectronsRecovery,fEvent->GetStartGtu(),fEvent->GetEndGtu(),3);
f1->SetParameters(outpar);
f2->SetParameters(outpar);
f1->Draw("same");
f2->Draw("same");
f1->SetLineStyle(1); f1->SetLineColor(1);
f2->SetLineStyle(2); f2->SetLineColor(2);
SaveFitInfo();
SafeDelete(f1);
SafeDelete(f2);
SafeDelete(minuit);
}
//_____________________________________________________________________________
void EnergyModule::FitMode2() {
TMinuit *minuit = new TMinuit(2);
minuit->SetPrintLevel(fMinuitOutputLevel);
minuit->SetObjectFit(&fData);
minuit->SetFCN(PhotoElectronsTimeFit2);
minuit->SetErrorDef(1.);
Double_t par[2] = {1.e-1,0.0};
Double_t step[3] = {1.e-1,1.e0};
Double_t min[3] = {0.00, 0.00};
Double_t max[3] = {0.00, 0.00};
string cpar[3] = {"Energy","tMax"};
for( Int_t i=0; i<2; i++ ) {
minuit->DefineParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
}
// minuit->FixParameter(1);
minuit->Migrad();
// retrieve results
Double_t outpar[2],err[2];
for( Int_t i=0; i<2; i++ ) {
minuit->GetParameter(i,outpar[i],err[i]);
}
fEnergy = outpar[0]*1e15;
fEnergyError = err[0]*1e15;
Double_t amin,edm,errdef;
Int_t nvpar,nparx,icstat;
minuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
fEnergyChi2 = amin/nvpar;
TF1 *f1 = new TF1("f1",PhotoElectrons2Total,fEvent->GetStartGtu(),fEvent->GetEndGtu(),2);
TF1 *f2 = new TF1("f2",PhotoElectrons2Recovery,fEvent->GetStartGtu(),fEvent->GetEndGtu(),2);
f1->SetParameters(outpar);
f2->SetParameters(outpar);
f1->Draw("same");
f2->Draw("same");
f1->SetLineStyle(1); f1->SetLineColor(1);
f2->SetLineStyle(2); f2->SetLineColor(2);
SaveFitInfo();
SafeDelete(f1);
SafeDelete(f2);
SafeDelete(minuit);
}
//_____________________________________________________________________________
void EnergyModule::FitMode2fumili() {
TFumili *fumili = new TFumili(2);
fumili->SetObjectFit(&fData);
fumili->SetFCN(PhotoElectronsTimeFit2);
fumili->SetErrorDef(1.);
Double_t par[2] = {1.e-1,0.0};
Double_t step[3] = {1.e-1,1.e0};
Double_t min[3] = {0.00, 0.00};
Double_t max[3] = {0.00, 0.00};
string cpar[3] = {"Energy","tMax"};
for( Int_t i=0; i<2; i++ ) {
fumili->SetParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
}
// fumili->FixParameter(1);
// fumili->Migrad();
fumili->Minimize();
// retrieve results
Double_t outpar[2],err[2];
char cname[100];
Double_t vlow[2], vhigh[2];
for( Int_t i=0; i<2; i++ ) {
fumili->GetParameter(i,cname,outpar[i],err[i],vlow[i],vhigh[i]);
}
fEnergy = outpar[0]*1e15;
fEnergyError = err[0]*1e15;
Double_t amin,edm,errdef;
Int_t nvpar,nparx;
fumili->GetStats(amin,edm,errdef,nvpar,nparx);
fEnergyChi2 = amin/nvpar;
TF1 *f1 = new TF1("f1",PhotoElectrons2Total,fEvent->GetStartGtu(),fEvent->GetEndGtu(),2);
TF1 *f2 = new TF1("f2",PhotoElectrons2Recovery,fEvent->GetStartGtu(),fEvent->GetEndGtu(),2);
f1->SetParameters(outpar);
f2->SetParameters(outpar);
f1->Draw("same");
f2->Draw("same");
f1->SetLineStyle(1); f1->SetLineColor(1);
f2->SetLineStyle(2); f2->SetLineColor(2);
SaveFitInfo();
SafeDelete(f1);
SafeDelete(f2);
SafeDelete(fumili);
}
//_____________________________________________________________________________
void EnergyModule::SaveFitInfo() {
TPaveText *cl = new TPaveText(0.78,0.6,1,1,"brNDC");
cl->SetTextFont(72); cl->SetTextSize(0.02); cl->SetFillColor(18); cl->SetTextAlign(12);
// Event Information
cl->AddText("Truth Information");
cl->AddText(Form("Energy = %.2E MeV",fEvent->GetHeader().GetTrueEnergy()));
cl->AddText(Form("#Theta = %.2f deg",fEvent->GetHeader().GetTrueTheta()*RadToDeg()));
cl->AddText(Form("#phi = %.2f deg",fEvent->GetHeader().GetTruePhi()*RadToDeg()));
cl->AddText("Reco Information");
cl->AddText(Form("Energy = %.2E #pm %2.E MeV",fEnergy,fEnergyError));
cl->AddText(Form("#Theta = %.2f deg",fTheta*RadToDeg()));
cl->AddText(Form("#phi = %.2f deg",fPhi*RadToDeg()));
cl->Draw();
if (fSaveEps == "yes")
if(fCanvas) fCanvas->Print(Form("output/check_timedistribution.run%d.event%d.eps",
fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
SafeDelete(cl);
}
//_____________________________________________________________________________
void EnergyModule::EstimateBackground() {
map <Int_t,Int_t> mapGtuCounts;
map <Int_t,Int_t> mapGtuPixels;
vector<RecoPixelData*> Bg;
for (UInt_t i(0); i<fAllActivePixels.size(); i++) {
// look if this pixel at the given time was taken as signal
RecoPixelData* pix = fAllActivePixels[i];
Bool_t IsPixelFound = kFALSE;
for (UInt_t j(0); j<fPixelsCluster.size(); j++) {
RecoPixelData* pix2 = fEvent->GetRecoPixelData(fPixelsCluster[j]);
if (pix->GetPixelId() == pix2->GetPixelId() && pix->GetGtu() == pix2->GetGtu()) {
if (pix->GetCounts() != pix2->GetCounts())
Msg(EsafMsg::Warning) << "EstimateBackground() two pixels with same Gtu and same Uid have different number of Counts. Is it possible?!" << MsgDispatch;
IsPixelFound = kTRUE;
break;
}
}
if (!IsPixelFound) Bg.push_back(pix);
}
for (UInt_t i(0); i<Bg.size(); i++) {
RecoPixelData* pix = Bg[i];
mapGtuCounts[pix->GetGtu()] += pix->GetCounts();
mapGtuPixels[pix->GetGtu()] += 1;
}
map <Int_t,Int_t>::iterator it;
Int_t point(0);
for (it = mapGtuCounts.begin(); it != mapGtuCounts.end(); it++) {
Float_t noise = 1.e0*it->second/mapGtuPixels[it->first];
fBgGraph->SetPoint(point++,it->first,noise);
}
if(fCanvas) fCanvas->Clear();
fCanvas->Divide(1,2);
fCanvas->cd(1);
fBgGraph->Fit("pol0","Q");
fBackground = fBgGraph->GetFunction("pol0")->GetParameter(0);
fBgGraph->Draw("A*");
fBgGraph->GetXaxis()->SetTitle("Gtu Hit");
fBgGraph->GetYaxis()->SetTitle("Noise/pixel/GTU");
fCanvas->cd(2);
TH1F* hNoise = new TH1F("hNoise","Noise/pixel/GTU",100,1,2);
Double_t *noise = fBgGraph->GetY();
for (Int_t i(0); i < fBgGraph->GetN(); i++) {
hNoise->Fill(noise[i]);
}
hNoise->Draw();
hNoise->GetXaxis()->SetTitle("Noise/pixel/GTU");
fCanvas->cd(0);
if (fSaveEps == "yes")
if(fCanvas) fCanvas->Print(Form("output/check_noisedistribution.run%d.event%d.eps",
fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
SafeDelete(hNoise);
}
//_____________________________________________________________________________
Bool_t EnergyModule::Process(RecoEvent *ev) {
// manager of the energy reconstruction chain
fEvent = ev;
fRunpars = fEvent->GetHeader().GetRunPars();
fQuantumEff = fRunpars->GetPmtData().GetQuantum();
CreateAuxilaryObjects();
// Load information from previous modules
LoadModules();
// Make map <GtuHit_t,Pixels_t>
SortGtu();
// Make maps <GtuHit_t,MeanTheta>, <GtuHit_t,MeanPhi>
BuildFoVAngleFunctions();
// Make map of pmts and pixels to display
DisplayCluster();
// Find Theta, Phi at Shower Maximum
FindThetaPhiAtMaximum();
// Estimate Fluorescence yield for the Shower Maximum
MakeFluorYield();
// Estimate Background
EstimateBackground();
// Build the recovery function along the track
BuildRecoveryFunction();
// Generate Shower Track, fill expected histo for flux of photons on EP
ShowerTrack3D();
// Fit the time distribution
FitTimeDistribution();
// Save data
MyData()->Add("Energy",fEnergy);
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
//
Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
void EnergyModule::UserMemoryClean() {
//
// Here we clear any auxilary objects created during one event energy reconstruction
//
fTheta = -1;
fPhi = -1;
fNumPoints = -1;
fAltitude = -1;
fGtuMin = 1.e20;
fGtuMax = -1.e20;
fQuantumEff = -1;
fTotalYield = -1;
fEnergy = -1;
fEnergyError = -1;
fEnergyChi2 = -1;
fFluorMaxPe = -1;
fDeltaL = -1;
fGtuAtMax = -1;
fGtuTimeAtMax = -1;
fThetaFoVAtMax = -1;
fPhiFoVAtMax = -1;
fMeanAttenuation = -1;
fMeanLambdaOnEntrancePupil = -1;
fRmax = -1;
fBackground = -1;
fXmax = -1.e20;
fXmin = 1.e20;
fYmax = -1.e20;
fYmin = 1.e20;
fData.Clear();
fMaxPos.Clear();
SafeDelete(fEPspectrum);
SafeDelete(fEPphotons);
fPixelsCluster.clear();
fAllActivePixels.clear();
fShowerDir.Clear();
fPmtLines.clear();
fPixelLines.clear();
fGtuPixels.clear();
fGtuMyPixels.clear();
}