// ESAF : Euso Simulation and Analysis Framework
// $Id: HmaxByShapeMethodModule.cc,v 1.13 2005/02/18 11:52:50 naumov Exp $
// Dmitry V. Naumov created Jun, 6 2004
#include "HmaxByShapeMethodModule.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "RecoRootEvent.hh"
#include "Atmosphere.hh"
#include "Config.hh"
#include "utils.hh"
#include <TProfile.h>
#include <TCanvas.h>
#include <TH1F.h>
#include <TF1.h>
ClassImp(HmaxByShapeMethodModule)
HmaxByShapeMethodModule* tmpHMaxpointer = NULL; // Points to the HmaxByShapeMethodModule class
// fFluorescence function
Double_t HmaxByShapeMethod_fFluorescence(Double_t *x, Double_t *par) {
return par[0]*exp(-0.5*pow(x[0]-par[1],2)/pow(par[2],2));
}
// Sum of background and peak functions
Double_t HmaxByShapeMethod_fitFunction(Double_t *x, Double_t *par) {
return HmaxByShapeMethod_fFluorescence(x,&par[0]) + tmpHMaxpointer->GetBackground();
}
//_____________________________________________________________________________
float HmaxByShapeMethod_ShowerMaximumDerivative(float t){
float meanlambda = 360; // nm
float xr = 2974/37.15; // g/cm2
float tmax = 1.7 + 0.76*log(tmpHMaxpointer->GetEnergy()/8.1e7);
float s = 2./(1.+tmax/t);
float CosTheta = cos(tmpHMaxpointer->GetTheta());
RecoPixelData* pix =
tmpHMaxpointer->GetClusters()[tmpHMaxpointer->PixelId(tmpHMaxpointer->GetFluorMaxTime())];
float ct = cos(pix->GetTheta());
return 1. - 2*log(s) - tmax/t*s - pow(400/meanlambda,4)/xr*CosTheta/ct;
}
//_____________________________________________________________________________
float HmaxByShapeMethod_DensityToAltitude(float h) {
return Atmosphere::Get()->Air_Density(h) - tmpHMaxpointer->GetDensity();
}
//_____________________________________________________________________________
HmaxByShapeMethodModule::HmaxByShapeMethodModule() : RecoModule("HmaxByShapeMethod") {
// ctor
tmpHMaxpointer = this;
fEnergy = 1.e20;
fTemp = NULL;
fFitFcn = NULL;
fFluor = NULL;
}
//_____________________________________________________________________________
HmaxByShapeMethodModule::~HmaxByShapeMethodModule() {
// dtor
if (!fc)
delete fc;
}
//_____________________________________________________________________________
Bool_t HmaxByShapeMethodModule::Init() {
Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
Atmosphere::Get();
Msg(EsafMsg::Info) << "HmaxByShapeMethodModule Enabled Atmosphere: " << Atmosphere::Get()->GetType() << MsgDispatch;
// Initialize the mean background
// FIXME:need to be read from the Event and not from config file
ConfigFileParser *pConfig = Config::Get()->GetCF("Electronics","EusoElectronics");
fMeanBackground = pConfig->GetNum("EusoElectronics.fNightGlowRate");
return kTRUE;
}
//_____________________________________________________________________________
Bool_t HmaxByShapeMethodModule::PreProcess() {
fc = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("hmaxbyshape");
if (!fc) fc = new TCanvas("hmaxbyshape","hmaxbyshape");
fc->Clear();
fc->cd();
fc->SetFillColor(33);
fc->SetFrameFillColor(41);
fc->SetGrid();
fTemp = NULL;
fFitFcn = NULL;
fFluor = NULL;
fEv = NULL;
fAltitude = -1;
fFluorMaxAt = -1;
fFluorMaxPe = -1;
fFluorMaxPeObs = -1;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t HmaxByShapeMethodModule::Process(RecoEvent *ev) {
fEv = ev;
const RecoModuleData *td = ev->GetModuleData("TrackDirection");
if ( td==NULL ) {
Msg(EsafMsg::Panic) << "TrackDirectionModule is not found." << MsgDispatch;
throw runtime_error( "HmaxByShapeMethodModule failed.");
}
fTheta = td->GetDouble("Theta");
fPhi = td->GetDouble("Phi");
fOmega.SetXYZ(TMath::Sin(fTheta)*TMath::Cos(fPhi), TMath::Sin(fTheta)*TMath::Sin(fPhi),
TMath::Cos(fTheta));
const RecoModuleData *gcm = ev->GetModuleData("GTUClustering");
if ( gcm==NULL ) {
Msg(EsafMsg::Panic) << "GTUClusteringModule is not found." << MsgDispatch;
throw runtime_error("HmaxByShapeMethodModule failed.");
}
else {
if ( gcm->GetObj("CluPixels") == NULL ) {
Msg(EsafMsg::Panic) << "No data in event" << MsgDispatch;
return kTRUE;
}
fPointsId = *(vector<Int_t>*)gcm->GetObj("CluPixels");
if (fPointsId.size() == 0) {
fNumPoints = 0;
Msg(EsafMsg::Panic) << "No pixels selected: terminated." << MsgDispatch;
return kTRUE;
} else {
fNumPoints = fPointsId.size();
Msg(EsafMsg::Debug) << "Processing " << fNumPoints << " points" << MsgDispatch;
}
}
fGtuLength = ev->GetHeader().GetGtuLength()/1000;
fFirstHit = ev->GetStartGtu();
fLastHit = ev->GetEndGtu();
Float_t Nbins = (Int_t) (fLastHit - fFirstHit)/fGtuLength;
fTemp = (TH1F*)gDirectory->FindObject("fTemp_HMaxShapeModule");
if ( fTemp != NULL) delete fTemp;
fTemp = new TH1F("fTemp_HMaxShapeModule","",(Int_t)Nbins,fFirstHit,fLastHit);
fTemp->GetXaxis()->SetTitle("time, #mu sec");
fTemp->GetYaxis()->SetTitle("PhotoElectrons");
// Save event unique tag
fEventTag = Form("hmaxshape.run.%d.event%d",ev->GetHeader().GetRun(),ev->GetHeader().GetNum());
for(Int_t i=0; i<fNumPoints; i++) {
RecoPixelData *pix = ev->GetRecoPixelData(fPointsId[i]);
fTemp->Fill(pix->GetGtu(),pix->GetCounts());
cluster.push_back(pix);
}
Double_t par[3]={1,1,1};
fFitFcn = new TF1("fFitFcn",HmaxByShapeMethod_fitFunction,ev->GetStartGtu(),ev->GetEndGtu(),3);
fFluor = new TF1("fFluor",HmaxByShapeMethod_fFluorescence,ev->GetStartGtu(),ev->GetEndGtu(),3);
fFitFcn->SetNpx(500);
fFitFcn->SetLineWidth(2);
fFitFcn->SetLineColor(kMagenta);
// find center of mass of fTemp
// Try to find the fFluorescent peak
float maxintfl(0);
int FlWidth = 10;
for(int i=0;i<(Int_t)fTemp->GetNbinsX();i++){
if(maxintfl < fTemp->Integral(i,i+FlWidth) ){
maxintfl = fTemp->Integral(i,i+FlWidth);
double ftime_fTemp(0), fstat(0);
for (int j=i;j<i+FlWidth;j++) {
ftime_fTemp += fTemp->GetBinContent(j)*fTemp->GetBinCenter(j);
fstat += fTemp->GetBinContent(j);
}
fFluorMaxTime = ftime_fTemp/fstat;
fFluorMaxPe = fTemp->GetBinContent(fTemp->FindBin(fFluorMaxTime));
}
}
fFitFcn->SetParLimits(0,0,fTemp->GetMaximum());
fFitFcn->SetParLimits(1,0,ev->GetEndGtu());
fFitFcn->SetParLimits(2,0,0.5*ev->GetEndGtu());
fFitFcn->SetParameters(fFluorMaxPe,fFluorMaxTime,fTemp->GetRMS());
fTemp->Fit("fFitFcn","RQ");
fChi2 = fTemp->GetFunction("fFitFcn")->GetChisquare()/fTemp->GetFunction("fFitFcn")->GetNDF();
fTemp->GetFunction("fFitFcn")->GetParameters(&par[0]);
fFitFcn->SetParameters(par);
fTemp->Draw("E");
fFitFcn->Draw("same");
fFluor->SetParameters(par);
fFluorMaxPe = par[0];
fFluorMaxTime = par[1];
if ( PixelId(fFluorMaxTime) == -1) return kFALSE;
GetAltitude();
MyData()->Add("Altitude",fAltitude*km);
MyData()->Add("FluorMaxAt",fFluorMaxAt);
MyData()->Add("FluorMaxPe",fFluorMaxPe);
MyData()->Add("FluorMaxPeObs",fTemp->GetMaximum());
MyData()->Add("DeltaL",DeltaL(fFluorMaxTime)*km);
return kTRUE;
}
//_____________________________________________________________________________
Bool_t HmaxByShapeMethodModule::PostProcess() {
if (fTemp!=(TH1F*)NULL) {
fTemp->SetMarkerStyle(21);
fTemp->SetMarkerSize(0.8);
fTemp->SetStats(0);
}
ConfigFileParser *pConfig = Config::Get()->GetCF("Reco","HmaxByShapeMethodModule");
string saveEps = pConfig->GetStr("HmaxByShapeMethodModule.SaveEps");
if (saveEps == "yes")
fc->Print("output/"+fEventTag+".eps");
return kTRUE;
}
//_____________________________________________________________________________
Bool_t HmaxByShapeMethodModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
fRecoRootEvent->GetRecoHmaxByShapeMethod().SetQuality(fChi2);
fRecoRootEvent->GetRecoHmaxByShapeMethod().SetHmax(fAltitude);
fRecoRootEvent->GetRecoHmaxByShapeMethod().SetTempHisto(fTemp);
fRecoRootEvent->GetRecoHmaxByShapeMethod().SetErrorHmax(fAltitude -
fEv->GetHeader().GetTrueShowerMaxPos().Z()/km);
return kTRUE;
}
//_____________________________________________________________________________
Bool_t HmaxByShapeMethodModule::Done() {
Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
void HmaxByShapeMethodModule::UserMemoryClean() {
//
//
//
SafeDelete(fFitFcn);
SafeDelete(fFluor);
SafeDelete(fTemp);
cluster.clear();
}
//_____________________________________________________________________________
Int_t HmaxByShapeMethodModule::PixelId(Double_t time) {
fFluorMaxAt = -1;
// This methods returns the pixel id number corresponding to the given time moment
// returns -1 if no pixel found (WARNING)
for(UInt_t i = 0; i < cluster.size() - 1; i++)
if (cluster[i]->GetGtu() < time && time <= cluster[i+1]->GetGtu()) {
fFluorMaxAt = i;
return i;
}
Msg(EsafMsg::Warning) << "PixelId can not determine pixel id corresponding to time " << time << MsgDispatch;
return fFluorMaxAt;
}
//_____________________________________________________________________________
Double_t HmaxByShapeMethodModule::nOmega(Double_t time) {
// This method computes the scalar product of the shower unit vector
// Omega (sin(Theta)*cos(Phi), sin(Theta)*sin(Phi), cos(Theta))
// and the given pixel unit vector aling its field of view direction
Int_t GtuHit = PixelId(time);
if (GtuHit == -1)
return -100;
TVector3 pixelUV;
RecoPixelData *pix = cluster[GtuHit];
Double_t st = TMath::Sin(pix->GetTheta());
Double_t ct = TMath::Cos(pix->GetTheta());
Double_t sp = TMath::Sin(pix->GetPhi());
Double_t cp = TMath::Cos(pix->GetPhi());
pixelUV.SetXYZ(st*cp,st*sp,ct);
return pixelUV*fOmega;
}
//_____________________________________________________________________________
Double_t HmaxByShapeMethodModule::DeltaL(Double_t time)
{
// This method returns the track length inclined with Theta and Phi zenith and azimuth angles
// respectively seen by the given pixel in one GTU time interval
return 0.3*fGtuLength/(1.+nOmega(time));
}
//_____________________________________________________________________________
void HmaxByShapeMethodModule::GetAltitude()
{
Float_t x0 = 37.15*g/cm2;
Float_t t1 = 10, t2 = 100;
Float_t (*f)(float)= &HmaxByShapeMethod_ShowerMaximumDerivative;
zbrac(f,&t1,&t2);
Float_t t0 = rtbis(f,t1,t2,0.01);
Double_t FSecondDerivative = -1/(2*t0);
Double_t FluorEntries = fFluor->Integral(fFirstHit,fLastHit)/fGtuLength;
Double_t erfLeft = TMath::Erf(sqrt(log(fFluorMaxPe/fFluor->Eval(fFirstHit))));
Double_t erfRight = TMath::Erf(sqrt(log(fFluorMaxPe/fFluor->Eval(fLastHit))));
fDensity = x0*fFluorMaxPe/FluorEntries*(erfLeft+erfRight)*
TMath::Sqrt(-TMath::Pi()/FSecondDerivative/2)/(DeltaL(fFluorMaxTime)*km);
DensityToAltitude();
}
//_____________________________________________________________________________
void HmaxByShapeMethodModule::DensityToAltitude() {
Float_t h1 = 0*km, h2 = 100*km;
Float_t (*f)(float)= &HmaxByShapeMethod_DensityToAltitude;
zbrac(f,&h1,&h2);
fAltitude = rtbis(f,h1,h2,0.01*km)/km;
}