// ESAF : Euso Simulation and Analysis Framework
// $Id: TrackDirectionModule.cc,v 1.18 2005/03/25 18:32:42 naumov Exp $
// R. Pesce created May, 11 2004
#include "RecoPixelData.hh"
#include "RecoEvent.hh"
#include "TrackDirectionModule.hh"
#include "EConst.hh"
#include "MedianFit.hh"
#include "LeastSquaresFit.hh"
#include "RecoRootEvent.hh"
#include "Config.hh"
#include <TVector3.h>
#include <TObject.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TGraph2D.h>
#include <TDirectory.h>
#include <TKey.h>
#include <TMath.h>
#include <TMinuit.h>
#include <TTree.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TGraph2D.h>
ClassImp(TrackData)
ClassImp(TrackDirectionModule)
// ChiSquare function for beta fitting
void ChiSquareTrack(Int_t &npar, Double_t *deriv, Double_t &f, Double_t *par, Int_t flag);
void MyChiSquareTrack(Int_t &npar, Double_t *deriv, Double_t &f, Double_t *par, Int_t flag);
Double_t TimeAtEuso(Double_t *x, Double_t *par);
//_____________________________________________________________________________
TrackDirectionModule::TrackDirectionModule() : RecoModule("TrackDirection"),
fTimeAtEuso(0) {
// ctor
fAll = NULL;
}
//_____________________________________________________________________________
TrackDirectionModule::~TrackDirectionModule() {
// dtor
delete fAll;
}
//_____________________________________________________________________________
Bool_t TrackDirectionModule::Init() {
fMinuitOutputLevel = (Int_t)Conf()->GetNum("TrackDirectionModule.fMinuitOutputLevel");
if ( fMinuitOutputLevel < -1 || fMinuitOutputLevel > 3 )
Msg(EsafMsg::Panic) << "Wrong config value TrackDirectionModule.fMinuitOutputLevel" << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirectionModule::PreProcess() {
Double_t kHUGE = 1.e20;
fAll = NULL;
fTheta = -kHUGE;
fPhi = -kHUGE;
fErrorTheta = -kHUGE;
fErrorPhi = -kHUGE;
fErrorDirection = -kHUGE;
fNumPoints = 0;
fNumHits = 0;
fNum = -1;
fBetaData.Clear();
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirectionModule::Process(RecoEvent *ev) {
fEv = ev;
const RecoModuleData *gcm = fEv->GetModuleData("GTUClustering");
if ( gcm == NULL ) {
Msg(EsafMsg::Panic) << "TrackDirectionModule: not found GTUClusteringModule." << MsgDispatch;
}
else {
if ( gcm->GetObj("CluPixels") == NULL ) {
Msg(EsafMsg::Warning) << "TrackDirectionModule: no data in event" << MsgDispatch;
return kFALSE;
}
fPointsId = *(vector<Int_t>*)gcm->GetObj("CluPixels");
if (fPointsId.size() == 0) {
fNumPoints = 0;
Msg(EsafMsg::Warning) << "TrackDirectionModule: no pixels selected: terminated." << MsgDispatch;
return kFALSE;
} else {
fNumPoints = fPointsId.size();
Msg(EsafMsg::Debug) << "TrackDirectionModule: processing " << fNumPoints << " points" << MsgDispatch;
}
}
const RecoModuleData *tdpm = fEv->GetModuleData("TrackDetectorPlane");
if ( gcm==NULL ) {
Msg(EsafMsg::Panic) << "TrackDirectionModule: not found TrackDetectorModule." << MsgDispatch;
} else {
Double_t nx = tdpm->GetDouble("NormX");
Double_t ny = tdpm->GetDouble("NormY");
Double_t nz = tdpm->GetDouble("NormZ");
Double_t wx = tdpm->GetDouble("WAxisX");
Double_t wy = tdpm->GetDouble("WAxisY");
Double_t wz = tdpm->GetDouble("WAxisZ");
Double_t ux = tdpm->GetDouble("UAxisX");
Double_t uy = tdpm->GetDouble("UAxisY");
Double_t uz = tdpm->GetDouble("UAxisZ");
Double_t xx = tdpm->GetDouble("XAxisX");
Double_t xy = tdpm->GetDouble("XAxisY");
Double_t xz = tdpm->GetDouble("XAxisZ");
Double_t yx = tdpm->GetDouble("YAxisX");
Double_t yy = tdpm->GetDouble("YAxisY");
Double_t yz = tdpm->GetDouble("YAxisZ");
fNormTDP.SetXYZ(nx,ny,nz);
fWAxisTDP.SetXYZ(wx,wy,wz);
fUAxisTDP.SetXYZ(ux,uy,uz);
fXAxisTDP.SetXYZ(xx,xy,xz);
fYAxisTDP.SetXYZ(yx,yy,yz);
vector<TVector3> *spPoints = (vector<TVector3>*)tdpm->GetObj("SpPoints");
fSpPos = *spPoints;
fNumHits = tdpm->GetInt("NumHits");
}
fNum = fEv->GetHeader().GetNum();
ConfigFileParser *pConfig = Config::Get()->GetCF("General","Euso");
fHISS = pConfig->GetNum("Euso.fAltitude")*km;
fEarthRadius = EConst::EarthRadius();
FillBetaData();
if (Conf()->GetStr("TrackDirectionModule.fMethod")=="approx") {
UseApprox();
}
else if (Conf()->GetStr("TrackDirectionModule.fMethod")=="exact") {
UseApprox();
UseExact();
}
else {
Msg(EsafMsg::Panic) << "Wrong config value TrackDirectionModule.fMethod" << MsgDispatch;
}
PrintBeta();
fTheta = TMath::Pi() - fTheta; // DN FIXME! Roberto uses a definition of theta opposite to default
MyData()->Add("Theta",fTheta);
MyData()->Add("Phi",fPhi);
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirectionModule::PostProcess() {
if ( fNumPoints == 0 ) return kTRUE;
fCenterPos = fEv->GetHeader().GetTrueShowerMaxPos();
Double_t dummyZ = fCenterPos.Z();
fCenterPos.SetZ(-dummyZ+fHISS);
fCenterPos = -fCenterPos;
fCenterX = fCenterPos.X();
fCenterY = fCenterPos.Y();
fDist = TMath::Sqrt(TMath::Power(fCenterPos.X(),2)+TMath::Power(fCenterPos.Y(),2));
fErrorPlane = TMath::ACos(fNormTDP*fTrueDir);
fVisualAngle = TMath::ACos(fTrueDir*(fCenterPos.Unit()));
fEASCenterX.push_back(fCenterPos.X());
fEASCenterY.push_back(fCenterPos.Y());
fDirectionErrors.push_back(fErrorDirection*TMath::RadToDeg());
fBetaErrors.push_back(fErrorBeta*TMath::RadToDeg());
fThetaErrors.push_back(fErrorTheta*TMath::RadToDeg());
fPhiErrors.push_back(fErrorPhi*TMath::RadToDeg());
fHMaxErrors.push_back(fErrorHMax/km);
fNumFastOr.push_back(fEv->GetHeader().GetTrueNumFastOR());
fNumPhEl.push_back(fNumHits);
fPlaneErrors.push_back(fErrorPlane*TMath::RadToDeg());
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirectionModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
fRecoRootEvent->GetRecoTrackDirection().SetTheta(fTheta);
fRecoRootEvent->GetRecoTrackDirection().SetErrorTheta(fErrorTheta);
fRecoRootEvent->GetRecoTrackDirection().SetPhi(fPhi);
fRecoRootEvent->GetRecoTrackDirection().SetErrorPhi(fErrorPhi);
fRecoRootEvent->GetRecoTrackDirection().SetErrorDirection(fErrorDirection);
fRecoRootEvent->GetRecoTrackDirection().SetQuality(0); // Put here a quality variable.
fRecoRootEvent->GetRecoTrackDirection().SetCheckInfo(fAll);
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDirectionModule::Done() {
Int_t numEv = (Int_t)fDirectionErrors.size();
if (numEv < 2) return kTRUE;
Float_t *xc, *yc, *dErr, *tErr, *pErr, *hErr, *bErr;
Float_t *fOr, *r, *plErr, *pEl;
xc = new Float_t[numEv];
yc = new Float_t[numEv];
hErr = new Float_t[numEv];
dErr = new Float_t[numEv];
bErr = new Float_t[numEv];
tErr = new Float_t[numEv];
pErr = new Float_t[numEv];
fOr = new Float_t[numEv];
pEl = new Float_t[numEv];
r = new Float_t[numEv];
plErr = new Float_t[numEv];
for( Int_t i=0; i<numEv; i++ ) {
xc[i] = fEASCenterX[i];
yc[i] = fEASCenterY[i];
dErr[i] = fDirectionErrors[i];
bErr[i] = fBetaErrors[i];
tErr[i] = fThetaErrors[i];
pErr[i] = fPhiErrors[i];
hErr[i] = fHMaxErrors[i];
fOr[i] = (Float_t)fNumFastOr[i];
pEl[i] = (Float_t)fNumPhEl[i];
r[i] = TMath::Sqrt(fEASCenterX[i]*fEASCenterX[i]+fEASCenterY[i]*fEASCenterY[i])/km;
plErr[i] = fPlaneErrors[i];
}
// events statistics
Int_t bestCounter(0);
Int_t medCounter(0);
Int_t worstCounter(0);
for( Int_t i=0; i<numEv; i++ ) {
if ( dErr[i] < 5.0 ) bestCounter++;
else if ( dErr[i] >= 5.0 && dErr[i] <= 10.0 ) medCounter++;
else
worstCounter++;
}
Msg(EsafMsg::Debug) << Form("STATISTICS on %d events",numEv) << MsgDispatch;
Msg(EsafMsg::Debug) << Form("Events with reco error under 5 deg = %d (%.3f %)",
bestCounter,(Float_t)bestCounter/numEv*100.) << MsgDispatch;
Msg(EsafMsg::Debug) << Form("Events with reco error between 5 and 10 deg = %d (%.3f %)",
medCounter,(Float_t)medCounter/numEv*100.) << MsgDispatch;
Msg(EsafMsg::Debug) << Form("Events with reco error above 10 deg = %d (%.3f %)",
worstCounter,(Float_t)worstCounter/numEv*100.) << MsgDispatch;
/*cout << "n-----Updating Error Histogram-----n" << endl;
if ( !fGeneralDir ) {
gDirectory->cd("/");
cout << "GeneralDir undefined. Looking for the key..." << flush;
TKey *kGeneralDir = gDirectory->FindKey("general");
if ( !kGeneralDir ) {
char pathname[] = "general";
cout << " search failed.nCreating new "<< pathname << " directory..." << flush;
fGeneralDir = new TDirectory(pathname,pathname);
} else {
cout << " succeded.nLoading object..." << flush;
fGeneralDir = (TDirectory*)kGeneralDir->ReadObj();
}
cout << ( fGeneralDir ? " done" : " failed" ) << endl;
}
fGeneralDir->cd();
cout << "Plotting " << numEv << " events..." << flush;
TGraph2D *recoerr = new TGraph2D(numEv, xc, yc, dErr);
recoerr->SetNameTitle("recoerr","Error in Reco Direction");
recoerr->Write();
delete recoerr;
TGraph2D *betaerr = new TGraph2D(numEv, xc, yc, bErr);
betaerr->SetNameTitle("betaerr","Error in Reco Beta");
betaerr->Write();
delete betaerr;
TGraph2D *thetaerr = new TGraph2D(numEv, xc, yc, tErr);
thetaerr->SetNameTitle("thetaerr","Error in Reco Theta");
thetaerr->Write();
delete thetaerr;
TGraph2D *phierr = new TGraph2D(numEv, xc, yc, pErr);
phierr->SetNameTitle("phierr","Error in Reco Phi");
phierr->Write();
delete phierr;
TGraph2D *hmaxerr = new TGraph2D(numEv, xc, yc, hErr);
hmaxerr->SetNameTitle("hmaxerr","Error in Reco HMax");
hmaxerr->Write();
delete hmaxerr;
TGraph *fastor = new TGraph(numEv,fOr,dErr);
fastor->SetNameTitle("fastor","RecoError vs NumFastOr");
fastor->SetMarkerSize(.5);
fastor->SetMarkerColor(kRed);
fastor->SetMarkerStyle(23);
fastor->Write();
delete fastor;
TGraph *phel = new TGraph(numEv,pEl,dErr);
phel->SetNameTitle("phel","RecoError vs Num P.E.");
phel->SetMarkerSize(.5);
phel->SetMarkerColor(kRed);
phel->SetMarkerStyle(23);
phel->Write();
delete phel;
TGraph *dist = new TGraph(numEv,r,dErr);
dist->SetNameTitle("dist","RecoError vs Dist from centre of FOV");
dist->SetMarkerSize(.5);
dist->SetMarkerColor(kRed);
dist->SetMarkerStyle(23);
dist->Write();
delete dist;
TGraph *plane = new TGraph(numEv,plErr,dErr);
plane->SetNameTitle("plane","RecoError vs PlaneError");
plane->SetMarkerSize(.5);
plane->SetMarkerColor(kRed);
plane->SetMarkerStyle(23);
plane->Write();
delete plane;
TGraph *planefastor = new TGraph(numEv,fOr,plErr);
planefastor->SetNameTitle("planefastor","PlaneError vs Num FastOr");
planefastor->SetMarkerSize(.5);
planefastor->SetMarkerColor(kRed);
planefastor->SetMarkerStyle(23);
planefastor->Write();
delete plane;
TGraph *planedist = new TGraph(numEv,r,plErr);
planedist->SetNameTitle("planedist","PlaneError vs Dist from center");
planedist->SetMarkerSize(.5);
planedist->SetMarkerColor(kRed);
planedist->SetMarkerStyle(23);
planedist->Write();
delete planedist;
TGraph *planephel = new TGraph(numEv,pEl,plErr);
planephel->SetNameTitle("planephel","PlaneError vs Num P.E.");
planephel->SetMarkerSize(.5);
planephel->SetMarkerColor(kRed);
planephel->SetMarkerStyle(23);
planephel->Write();
delete planephel;
TH1F *numErr = new TH1F("numErr", "Num Events vs RecoError", 100, 0, 180 );
TH1F *plnumErr = new TH1F("plnumErr", "Num Events vs PlaneError", 100, 0, 180 );
TH1F *nThetaErr = new TH1F("nThetaErr", "Num Events vs ThetaError", 100, -180, 180 );
TH1F *nPhiErr = new TH1F("nPhiErr", "Num Events vs PhiError", 100, -180, 180 );
for( Int_t i=0; i<numEv; i++ ) {
numErr->Fill(dErr[i]);
plnumErr->Fill(plErr[i]);
nThetaErr->Fill(tErr[i]);
nPhiErr->Fill(pErr[i]);
}
numErr->Write();
plnumErr->Write();
nThetaErr->Write();
nPhiErr->Write();
delete numErr;
delete plnumErr;
delete nThetaErr;
delete nPhiErr;
*/
delete [] xc;
delete [] yc;
delete [] dErr;
delete [] tErr;
delete [] pErr;
delete [] hErr;
delete [] bErr;
delete [] fOr;
delete [] r;
delete [] pEl;
delete [] plErr;
return kTRUE;
}
//_____________________________________________________________________________
void TrackDirectionModule::UserMemoryClean() {
fAll = NULL;
}
//_____________________________________________________________________________
void TrackDirectionModule::FillBetaData() {
fBetaData.fNumPoints = fNumPoints;
fBetaData.fNumHits = fNumHits;
for( Int_t i=0; i<fNumPoints; i++ ) {
RecoPixelData *pix = fEv->GetRecoPixelData(fPointsId[i]);
fBetaData.fAlpha.push_back(TMath::ATan2(fSpPos[i]*fUAxisTDP,fSpPos[i]*fWAxisTDP));
fBetaData.fTime.push_back((Double_t)2500.*pix->GetGtu()*ns);
fBetaData.fHits.push_back(pix->GetCounts());
fBetaData.fErrors.push_back(2500.*ns);
}
// search gtu min and max and with max nuber of hits
Double_t tmin = fBetaData.fTime[0];
Double_t tmax = fBetaData.fTime[0];
Int_t maxhits = fBetaData.fHits[0];
Int_t minid = 0;
Int_t maxid = 0;
Int_t mostpopid = 0;
for( Int_t i=1; i<fNumPoints; i++ ) {
if ( fBetaData.fTime[i] < tmin ) {
tmin = fBetaData.fTime[i];
minid = i;
}
if ( fBetaData.fTime[i] > tmax ) {
tmax = fBetaData.fTime[i];
maxid = i;
}
if ( fBetaData.fHits[i] > maxhits ) {
maxhits = fBetaData.fHits[i];
mostpopid = i;
}
}
vector<Int_t> maxCounts;
for( Int_t i=0; i<fNumPoints; i++ ) {
if ( fBetaData.fHits[i] == maxhits ) maxCounts.push_back(i);
}
Double_t sum(0);
for( size_t i=0; i<maxCounts.size(); i++ ) {
sum += fBetaData.fAlpha[maxCounts[i]];
}
sum /= maxCounts.size();
vector<Double_t> diffAlpha;
for( size_t i=0; i<maxCounts.size(); i++ ) {
diffAlpha.push_back(fBetaData.fAlpha[maxCounts[i]]-sum);
}
Int_t mpid(0);
for( size_t i=1; i<maxCounts.size(); i++ ) {
if ( diffAlpha[i] < diffAlpha[mpid] ) mpid=i;
}
mostpopid = maxCounts[mpid];
Double_t alphamin = TMath::Abs(fBetaData.fAlpha[minid]);
Double_t alphamax = TMath::Abs(fBetaData.fAlpha[maxid]);
if ( alphamin < alphamax) {
for( Int_t i=0; i<fNumPoints; i++ ) {
if ( fBetaData.fTime[i] == tmin ) {
if ( fBetaData.fAlpha[i] < alphamin ) {
alphamin = fBetaData.fAlpha[i];
minid = i;
}
}
if ( fBetaData.fTime[i] == tmax ) {
if ( fBetaData.fAlpha[i] > alphamax ) {
alphamax = fBetaData.fAlpha[i];
maxid = i;
}
}
}
} else if ( alphamax > alphamin ) {
for( Int_t i=0; i<fNumPoints; i++ ) {
if ( fBetaData.fTime[i] == tmin ) {
if ( fBetaData.fAlpha[i] > alphamin ) {
alphamin = fBetaData.fAlpha[i];
minid = i;
}
}
if ( fBetaData.fTime[i] == tmax ) {
if ( fBetaData.fAlpha[i] < alphamax ) {
alphamax = fBetaData.fAlpha[i];
maxid = i;
}
}
}
}
TVector3 pMax = fSpPos[mostpopid];
fBetaData.fIdMostPop = mostpopid;
fAlphaMax = fBetaData.fAlpha[fBetaData.fIdMostPop];
fTMax = fBetaData.fTime[fBetaData.fIdMostPop];
// MC truth
fTrueEnergy = fEv->GetHeader().GetTrueEnergy();
//Double_t zInit = fEv->GetHeader().GetTrueInitPos().Z();
fTrueHMax = fEv->GetHeader().GetTrueShowerMaxPos().Z();
fTrueThetaInc = fEv->GetHeader().GetTrueTheta();
fTrueTheta = TMath::Pi()-fEv->GetHeader().GetTrueTheta();
fTruePhi = fEv->GetHeader().GetTruePhi();
fTrueDir.SetMagThetaPhi(1.,fTrueTheta,fTruePhi);
fTrueBeta = TMath::ATan2(fTrueDir*fUAxisTDP, fTrueDir*fWAxisTDP);
//TVector3 perp = fNormTDP.Cross(fTrueDir);
//TVector3 initpos = fEv->GetHeader().GetTrueInitPos();
//fRMax = CalculateRMax(fTrueHMax);
fHMax = 5.*km;
fRMax = CalculateRMax(fHMax);
// create the function used to fit
fTimeAtEuso = new TF1("timeteo", TimeAtEuso, -TMath::Pi(), TMath::Pi(), 5);
// range's too wide. Get the correct interval
if ( fBetaData.fAlpha[mostpopid] > 0)
fTimeAtEuso->SetRange(TMath::Pi()/4., 3./4.*TMath::Pi());
else
fTimeAtEuso->SetRange(-3./4.*TMath::Pi(), -TMath::Pi()/4.);
// set parameters
fTimeAtEuso->SetParName(0,"h_max");
fTimeAtEuso->SetParName(1,"beta");
fTimeAtEuso->SetParName(2,"alpha_max");
fTimeAtEuso->SetParName(3,"t_max");
fTimeAtEuso->SetParName(4,"r_max");
// fix const params
fTimeAtEuso->SetParameter(0, 5.*km);
fTimeAtEuso->SetParameter(1, fTrueBeta);
fTimeAtEuso->FixParameter(2, fAlphaMax);
fTimeAtEuso->FixParameter(3, fTMax);
fTimeAtEuso->FixParameter(4, fRMax);
Msg(EsafMsg::Debug) << Form("fTimeAtEuso range: alpha_max = %.3f, min = %.3f, max = %.3f",
fBetaData.fAlpha[mostpopid], fTimeAtEuso->GetXmin(), fTimeAtEuso->GetXmax()) << MsgDispatch;
fBetaData.fIdMin = minid;
fBetaData.fIdMax = maxid;
fBetaData.fNorm = fNormTDP;
fBetaData.fWAxis = fWAxisTDP;
fBetaData.fUAxis = fUAxisTDP;
fBetaData.fTimeAtEuso = fTimeAtEuso;
fTimeAtEuso->Eval(fBetaData.fAlpha[mostpopid]);
fBetaData.fCsi = TMath::ACos((pMax.Unit()).Z());
//calculate omega mean
fOmegaMean = 0;
Int_t used = 0;
for( Int_t i(0); i<fNumPoints; i++ ) {
if (fBetaData.fTime[i]!=fTMax) {
fOmegaMean += (fBetaData.fAlpha[i]-fAlphaMax)/(fBetaData.fTime[i]-fTMax);
used++;
}
}
fOmegaMean /= used;
//calculate omega first and last
fOmegaFirst = (fBetaData.fAlpha[minid]-fAlphaMax)/(fBetaData.fTime[minid]-fTMax);
fOmegaLast = (fBetaData.fAlpha[maxid]-fAlphaMax)/(fBetaData.fTime[maxid]-fTMax);
}
//_____________________________________________________________________________
void TrackDirectionModule::UseApprox() {
// median fit in order to find an approximate value of beta to initialize fit
LeastSquaresFit *lBeta = new LeastSquaresFit(fNumPoints, fBetaData.fAlpha, fBetaData.fTime, fBetaData.fErrors);
MedianFit *mBeta = new MedianFit(fNumPoints, fBetaData.fTime, fBetaData.fAlpha);
Double_t visual = 2*TMath::ATan(1/(mBeta->GetSlope()*fRMax/EConst::C()));
fBeta = fBetaData.fAlpha[fBetaData.fIdMostPop] - visual - TMath::Pi();
if ( fBeta > 1.e11 ) fBeta -= TMath::TwoPi() *1.e10;
Int_t bPi = (Int_t) (fBeta/TMath::TwoPi());
fBeta -= TMath::TwoPi()*bPi;
fBetaInit = fBeta;
fOmega = mBeta->GetSlope();
fOmegaLeast = 1./lBeta->GetSlope();
MsgForm(EsafMsg::Debug, "__Beta Pre-Fit_____________________________________________________");
MsgForm(EsafMsg::Debug, "Slope = %.3e", mBeta->GetSlope());
MsgForm(EsafMsg::Debug, "Offset = %.3e", mBeta->GetOffset());
MsgForm(EsafMsg::Debug, "Absolute Deviation = %.3f", mBeta->GetAbsoluteDeviation());
MsgForm(EsafMsg::Debug, "RMax init = %.3f km", fRMax/km);
MsgForm(EsafMsg::Debug, "Visual Angle = %.3f deg", visual*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, "Beta Init = %.3f deg", fBeta*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, "____________________________________________________________________");
fDir = TMath::Cos(fBeta) * fWAxisTDP + TMath::Sin(fBeta) * fUAxisTDP;
fTheta = fDir.Theta();
fPhi = fDir.Phi();
if ( fPhi < - TMath::Pi() ) fPhi += TMath::Pi();
if ( fPhi > TMath::Pi() ) fPhi -= TMath::Pi();
fThetaInit = fTheta;
fPhiInit = fPhi;
// Error in reco direction
fErrorDirection = TMath::ACos(TMath::Sin(fTheta)*TMath::Sin(fTrueTheta)
*TMath::Cos(fPhi)*TMath::Cos(fTruePhi) +
TMath::Sin(fTheta)*TMath::Sin(fTrueTheta)
*TMath::Sin(fPhi)*TMath::Sin(fTruePhi) +
TMath::Cos(fTheta)*TMath::Cos(fTrueTheta));
fErrorDirectionInit = fErrorDirection;
fErrorBetaInit = fBetaInit - fTrueBeta;
fErrorThetaInit = fThetaInit - fTrueTheta;
fErrorPhiInit = fPhiInit - fTruePhi;
MsgForm(EsafMsg::Debug, "__Angular Velocity__________________________________________________");
MsgForm(EsafMsg::Debug, "Fitted (median) = %.3e rad/ns", fOmega);
MsgForm(EsafMsg::Debug, "Fitted (least squares) = %.3e rad/ns", fOmegaLeast);
MsgForm(EsafMsg::Debug, "Mean = %.3e rad/ns", fOmegaMean);
MsgForm(EsafMsg::Debug, "First = %.3e rad/ns", fOmegaFirst);
MsgForm(EsafMsg::Debug, "Last = %.3e rad/ns", fOmegaLast);
MsgForm(EsafMsg::Debug, "____________________________________________________________________");
}
//_____________________________________________________________________________
void TrackDirectionModule::UseExact() {
// correction of HMax
fHMax = 7.5*km * TMath::Log(1033./850./TMath::Abs(TMath::Cos(fTheta)));
fRMax = CalculateRMax(fHMax);
fTimeAtEuso->FixParameter(4, fRMax);
// initialize first fit
TMinuit *minbeta = new TMinuit(2);
minbeta->SetPrintLevel(fMinuitOutputLevel);
minbeta->SetObjectFit(&fBetaData);
minbeta->SetFCN(MyChiSquareTrack);
minbeta->SetErrorDef(1.);
// initialize parameters
Double_t par[2];
// For debug purpose: initialize parameters with true values
// par[0] = fTrueHMax;
// par[1] = fTrueBeta;
par[0] = fHMax;
par[1] = fBeta;
Double_t step[2] = {0.000, 0.001};
Double_t min[2] = {0.0, 0.0};
Double_t max[2] = {0.0, 0.0};
string cpar[2] = {"H_max","Beta"};
for( Int_t i=0; i<2; i++ ) {
minbeta->DefineParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
}
//Double_t plist[1] = {2.};
//Int_t ierflg;
//minbeta->mnexcm("SET STRategy",plist,1,ierflg);
//plist[0] = 1;
//minbeta->mnexcm("FIX",plist,1,ierflg);
minbeta->Migrad();
// retrieve results
Double_t outpar[2],err[2];
for( Int_t i=0; i<2; i++ ) {
minbeta->GetParameter(i,outpar[i],err[i]);
}
fBeta = outpar[1];
fHMax = outpar[0];
if ( fBeta > 1.e15 ) fBeta = fBetaInit;
Int_t bPi = (Int_t) (fBeta/TMath::TwoPi());
fBeta -= TMath::TwoPi()*bPi;
MsgForm(EsafMsg::Debug, "n___First beta fit: hmax = %.3f km", fHMax/km);
MsgForm(EsafMsg::Debug, "___First beta fit: beta = %.3f degn", fBeta*TMath::RadToDeg());
fDir = TMath::Cos(fBeta) * fWAxisTDP + TMath::Sin(fBeta) * fUAxisTDP;
fTheta = fDir.Theta();
fPhi = fDir.Phi();
if ( fPhi < - TMath::Pi() ) fPhi += TMath::Pi();
if ( fPhi > TMath::Pi() ) fPhi -= TMath::Pi();
// correction of HMax
fHMax = 7.5*km * TMath::Log(1033./850./TMath::Abs(TMath::Cos(fTheta)));
fRMax = CalculateRMax(fHMax);
fTimeAtEuso->FixParameter(4, fRMax);
// second fit
par[0] = fHMax;
par[1] = fBeta;
for( Int_t i=0; i<2; i++ ) {
minbeta->DefineParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
}
minbeta->Migrad();
for( Int_t i=0; i<2; i++ ) {
minbeta->GetParameter(i,outpar[i],err[i]);
}
fBeta = outpar[1];
fHMax = outpar[0];
if ( fBeta > 1.e11 ) fBeta = fBetaInit;
bPi = (Int_t) (fBeta/TMath::TwoPi());
fBeta -= TMath::TwoPi()*bPi;
fDir = TMath::Cos(fBeta) * fWAxisTDP + TMath::Sin(fBeta) * fUAxisTDP;
fTheta = fDir.Theta();
fPhi = fDir.Phi();
if ( fPhi < - TMath::Pi() ) fPhi += TMath::Pi();
if ( fPhi > TMath::Pi() ) fPhi -= TMath::Pi();
// Error in reco direction
fErrorDirection = TMath::ACos(TMath::Sin(fTheta)*TMath::Sin(fTrueTheta)
*TMath::Cos(fPhi)*TMath::Cos(fTruePhi) +
TMath::Sin(fTheta)*TMath::Sin(fTrueTheta)
*TMath::Sin(fPhi)*TMath::Sin(fTruePhi) +
TMath::Cos(fTheta)*TMath::Cos(fTrueTheta));
}
//______________________________________________________________________________
Double_t TrackDirectionModule::CalculateRMax(Double_t hmax) {
Double_t rmax;
Double_t thmax = fEv->GetRecoPixelData(fPointsId[fBetaData.fIdMostPop])->GetTheta();
rmax = (fEarthRadius+fHISS)*TMath::Cos(thmax)-
TMath::Sqrt(TMath::Power(fEarthRadius+hmax,2)
-TMath::Power((fEarthRadius+fHISS)*TMath::Sin(thmax),2));
return rmax;
}
//______________________________________________________________________________
void TrackDirectionModule::PrintBeta() {
MsgForm(EsafMsg::Debug, "_EAS Direction Fit Results_____________________________________________________");
MsgForm(EsafMsg::Debug, " TrueEnergy = %.3e eV ", fTrueEnergy);
MsgForm(EsafMsg::Debug, " TrueDir = (%.3f, %.3f, %.3f)", fTrueDir[0],fTrueDir[1],fTrueDir[2]);
MsgForm(EsafMsg::Debug, " Dir = (%.3f, %.3f, %.3f)", fDir[0],fDir[1],fDir[2]);
MsgForm(EsafMsg::Debug, " Theta = %.3f deg TrueTheta = %.3f deg",
fTheta*TMath::RadToDeg(), fTrueTheta*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " Phi = %.3f deg TruePhi = %.3f deg",
fPhi*TMath::RadToDeg(), fTruePhi*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " Beta = %.3f deg TrueBeta = %.3f deg",
fBeta*TMath::RadToDeg(), fTrueBeta*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " Error in Direction = %f deg", fErrorDirection*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " HMax = %.3f km TrueHMax = %.3f km", fHMax/km, fTrueHMax/km);
MsgForm(EsafMsg::Debug, "_______________________________________________________________________________");
fThetaInc = TMath::Pi()-fTheta;
fErrorTheta = fTheta - fTrueTheta;
fErrorPhi = fPhi - fTruePhi;
if ( fErrorPhi < -TMath::Pi() ) fErrorPhi += TMath::TwoPi();
fErrorHMax = fHMax - fTrueHMax;
fErrorBeta = fBeta - fTrueBeta;
Double_t t_ref = fBetaData.fTime[fBetaData.fIdMostPop];
Double_t t_prop, Dt, t_prop_tr, Dt_tr;
TVector3 initpos = fEv->GetHeader().GetTrueInitPos();
Double_t dummyZ = initpos.Z();
initpos.SetZ(-dummyZ+fHISS);
TVector3 trueMaxPos = fEv->GetHeader().GetTrueShowerMaxPos();
dummyZ = trueMaxPos.Z();
trueMaxPos.SetZ(-dummyZ+fHISS);
trueMaxPos = -trueMaxPos;
TVector3 perp = (fBetaData.fNorm.Cross(fDir)).Unit();
TVector3 truenorm = (fTrueDir.Cross(initpos)).Unit();
if ( truenorm.Z() < 0 ) truenorm *= -1;
TVector3 trueperp = (truenorm.Cross(fTrueDir)).Unit();
Double_t trueCsi = TMath::ACos((trueMaxPos.Unit()).Z());
fTrueAlphaMax = TMath::ATan2(trueMaxPos*fUAxisTDP, trueMaxPos*fWAxisTDP);
fTrueRMax = (-(fEarthRadius+fHISS)*TMath::Cos(trueCsi)-
TMath::Sqrt(TMath::Power(fEarthRadius+fTrueHMax,2)
-TMath::Power((fEarthRadius+fHISS)*TMath::Sin(trueCsi),2)));
fTrueR0 = fTrueRMax*TMath::Sin(fTrueAlphaMax-fTrueBeta);
fR0 = fRMax*TMath::Sin(fAlphaMax-fBeta);
fAlpha0 = TMath::ATan2(perp*fBetaData.fUAxis, perp*fBetaData.fWAxis);
fTrueAlpha0 = TMath::ATan2(trueperp*fBetaData.fUAxis, trueperp*fBetaData.fWAxis);
Double_t dAlpha = fAlphaMax-fAlpha0;
Double_t dAlpha_tr = fTrueAlphaMax-fTrueAlpha0;
fTau0 = t_ref-TMath::Abs(fR0/EConst::C()*(1/TMath::Cos(dAlpha)))
+fR0/EConst::C()*TMath::Tan(dAlpha);
fTrueTau0 = t_ref-TMath::Abs(fTrueR0/EConst::C()*(1/TMath::Cos(dAlpha_tr)))
+fTrueR0/EConst::C()*TMath::Tan(dAlpha_tr);
//TVector3 testdir;
//testdir.SetMagThetaPhi(1.,TMath::Pi()-fTrueDir.Theta(),fTrueDir.Phi());
MsgForm(EsafMsg::Debug, "_______________________________________________________________________________");
MsgForm(EsafMsg::Debug, " Norm = (%.3f, %.3f, %.3f)", fNormTDP[0],fNormTDP[1],fNormTDP[2]);
MsgForm(EsafMsg::Debug, " TrueNorm = (%.3f, %.3f, %.3f)", truenorm[0],truenorm[1],truenorm[2]);
MsgForm(EsafMsg::Debug, " Perp = (%.3f, %.3f, %.3f)", perp[0],perp[1],perp[2]);
MsgForm(EsafMsg::Debug, " TruePerp = (%.3f, %.3f, %.3f)", trueperp[0],trueperp[1],trueperp[2]);
MsgForm(EsafMsg::Debug, " TrueMax = (%.3f, %.3f, %.3f)", trueMaxPos[0],trueMaxPos[1],trueMaxPos[2]);
MsgForm(EsafMsg::Debug, " RMax = %.3f km", fRMax/km);
MsgForm(EsafMsg::Debug, " True RMax = %.3f km", fTrueRMax/km);
MsgForm(EsafMsg::Debug, " R0 = %.3f km", fR0/km);
MsgForm(EsafMsg::Debug, " True R0 = %.3f km", fTrueR0/km);
MsgForm(EsafMsg::Debug, " test: t_max = %.3f microsecond", t_ref/microsecond);
MsgForm(EsafMsg::Debug, " Alpha_max = %.3f deg", fAlphaMax*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " True Alpha_max = %.3f deg", fTrueAlphaMax*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " Csi = %.3f deg", fBetaData.fCsi*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " True Csi = %.3f deg", trueCsi*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " Alpha_0 = %.3f deg", fAlpha0*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " True Alpha_0 = %.3f deg", fTrueAlpha0*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " dAlpha = %.3f deg", dAlpha*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " True dAlpha = %.3f deg", dAlpha_tr*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " tau0 = %.3f microsecond", fTau0/microsecond);
MsgForm(EsafMsg::Debug, " True tau0 = %.3f microsecond", fTrueTau0/microsecond);
MsgForm(EsafMsg::Debug, " Visual Angle = %.3f deg",(fAlphaMax-fBeta-TMath::Pi())*TMath::RadToDeg() );
MsgForm(EsafMsg::Debug, " True Visual Angle 1 = %.3f deg",(fTrueAlphaMax-fTrueBeta-TMath::Pi())*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, " True Visual Angle 2 = %.3f deg",(TMath::ACos((-trueMaxPos.Unit())*fTrueDir)-TMath::TwoPi())*TMath::RadToDeg());
MsgForm(EsafMsg::Debug, "_______________________________________________________________________________");
Double_t *alpha, *t_teo, *t_sper, *t_trueteo, *hits;
alpha = new Double_t[fNumPoints];
t_teo = new Double_t[fNumPoints];
t_sper = new Double_t[fNumPoints];
t_trueteo = new Double_t[fNumPoints];
hits = new Double_t[fNumPoints];
TF1* trueTimeAtEuso = new TF1(*fTimeAtEuso);
trueTimeAtEuso->SetParameter(0, fTrueHMax);
trueTimeAtEuso->SetParameter(1, fTrueBeta);
trueTimeAtEuso->SetParameter(2, fTrueAlphaMax);
trueTimeAtEuso->SetParameter(3, t_ref);
trueTimeAtEuso->SetParameter(4, fTrueRMax);
fTimeAtEuso->SetParameter(0, fHMax);
fTimeAtEuso->SetParameter(1, fBeta);
for (Int_t i(0); i<fNumPoints ; i++) {
dAlpha = fBetaData.fAlpha[i]-fAlpha0;
dAlpha_tr = fBetaData.fAlpha[i]-fTrueAlpha0;
t_prop = TMath::Abs(fR0/EConst::C()*1/TMath::Cos(dAlpha));
Dt = fR0/EConst::C()*(-TMath::Sin(dAlpha))/TMath::Cos(dAlpha);
t_prop_tr = TMath::Abs(fTrueR0/EConst::C()*1/TMath::Cos(dAlpha_tr));
Dt_tr = fTrueR0/EConst::C()*(-TMath::Sin(dAlpha_tr))/TMath::Cos(dAlpha_tr);
hits[i] = fBetaData.fHits[i];
alpha[i] = fBetaData.fAlpha[i]*TMath::RadToDeg();
t_sper[i] = fBetaData.fTime[i];
// t_teo[i] = tau0 + Dt + t_prop;
// t_trueteo[i] = tau0_tr + Dt_tr + t_prop_tr;
t_teo[i] = fTimeAtEuso->Eval(fBetaData.fAlpha[i]);
t_trueteo[i] = trueTimeAtEuso->Eval(fBetaData.fAlpha[i]);
}
// sort arrays
Int_t j,k;
Double_t a, tt, ts, ttt;
for(k=1; k<fNumPoints; k++) {
a = alpha[k];
tt = t_teo[k];
ts = t_sper[k];
ttt = t_trueteo[k];
j=k-1;
while( j>=0 && alpha[j]>a ){
alpha[j+1]=alpha[j];
t_teo[j+1]=t_teo[j];
t_sper[j+1]=t_sper[j];
t_trueteo[j+1]=t_trueteo[j];
j--;
}
alpha[j+1]=a;
t_teo[j+1]=tt;
t_sper[j+1]=ts;
t_trueteo[j+1]=ttt;
}
TGraph *gTimeTeo = new TGraph(fNumPoints, alpha, t_teo);
gTimeTeo->SetLineWidth(3);
gTimeTeo->SetLineColor(kBlue);
// gTimeTeo->SetMarkerStyle(20);
TGraph *gTimeSper = new TGraph(fNumPoints, alpha, t_sper);
gTimeSper->SetMarkerSize(.5);
gTimeSper->SetMarkerColor(kRed);
gTimeSper->SetMarkerStyle(23);
TGraph *gTimeTrueTeo = new TGraph(fNumPoints, alpha, t_trueteo);
gTimeTrueTeo->SetLineWidth(1);
gTimeTrueTeo->SetLineColor(kGreen);
if (!fAll)
fAll = new TMultiGraph("BetaCheck","Comparison between t_sper and t_teo");
else
fAll->Clear();
fAll->Add(gTimeSper,"p");
fAll->Add(gTimeTeo,"l");
fAll->Add(gTimeTrueTeo,"l");
/* delete gTimeTeo;
delete gTimeSper;
delete gTimeTrueTeo;
*/
delete [] t_teo;
delete [] t_trueteo;
// TGraph2D *gHits = new TGraph2D(fNumPoints, alpha, t_sper, hits);
// gHits->SetNameTitle("hits", "Hits of Pixels");
// gHits->Write();
delete [] alpha;
delete [] t_sper;
delete [] hits;
// delete gHits;
}
//______________________________________________________________________________
/*void ChiSquareTrack(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
f = 0;
TrackData *tData = (TrackData*) gMinuit->GetObjectFit();
Double_t t_teo;
Double_t t_ref = tData->fTime[tData->fIdMostPop];
Double_t alpha_ref = tData->fAlpha[tData->fIdMostPop];
for(Int_t i(0); i<tData->fNumPoints; i++) {
t_teo = t_ref - (tData->fHISS - par[0])/(TMath::Cos(tData->fCsi)*EConst::C()) *
TMath::Sin(alpha_ref-par[1])/TMath::Tan(par[1]/2.) * (tData->fAlpha[i] - alpha_ref);
f += tData->fHits[i]*TMath::Power(((tData->fTime[i]-t_teo)/tData->fErrors[i]),2);
}
f = (Double_t) f/(tData->fNumHits - 2.);
}*/
//______________________________________________________________________________
void ChiSquareTrack(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
f = 0;
TrackData *bData = (TrackData*) gMinuit->GetObjectFit();
Double_t t_teo;
//TVector3 dir = TMath::Cos(par[1])*bData->fWAxis+TMath::Sin(par[1])*bData->fUAxis;
//TVector3 perp = bData->fNorm.Cross(dir);
//Double_t alpha0 = TMath::ATan2(perp*bData.fUAxis, perp*bData.fWAxis);
Double_t alpha_ref = bData->fAlpha[bData->fIdMostPop];
Double_t alpha0 = par[1] + TMath::PiOver2();
Double_t dAlpha = alpha_ref - alpha0;
Double_t t_ref = bData->fTime[bData->fIdMostPop];
Double_t R0 = (430.*km - par[0])/TMath::Cos(bData->fCsi)*TMath::Sin(alpha_ref-par[1]);
Double_t Dt = R0/EConst::C()*(-TMath::Sin(dAlpha))/TMath::Cos(dAlpha);
Double_t t_prop = TMath::Abs(R0/EConst::C()*(1)/TMath::Cos(dAlpha));
Double_t tau0 = t_ref-(Dt+t_prop);
//Double_t tau0 = t_ref-R0/EConst::C()*(1-TMath::Sin(dAlpha))/TMath::Cos(dAlpha);
// Double_t tau0 = t_ref-par[0]/EConst::C()*(1-TMath::Sin(dAlpha))/TMath::Cos(dAlpha);
for (Int_t i(0); i<bData->fNumPoints ; i++) {
dAlpha = bData->fAlpha[i]-alpha0;
t_prop = TMath::Abs(R0/EConst::C()*(1)/TMath::Cos(dAlpha));
Dt = R0/EConst::C()*(-TMath::Sin(dAlpha))/TMath::Cos(dAlpha);
t_teo = tau0+Dt+t_prop;
f+=bData->fHits[i]*TMath::Power(((bData->fTime[i]-t_teo)/bData->fErrors[i]),2);
}
f = (Double_t) f/(bData->fNumHits - 2.);
}
//______________________________________________________________________________
void MyChiSquareTrack(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
f = 0;
TrackData *bData = (TrackData*) gMinuit->GetObjectFit();
Double_t t_teo;
bData->fTimeAtEuso->SetParameter("h_max", par[0]);
bData->fTimeAtEuso->SetParameter("beta", par[1]);
for (Int_t i(0); i<bData->fNumPoints ; i++) {
t_teo = bData->fTimeAtEuso->Eval(bData->fAlpha[i]) ;
f+=bData->fHits[i]*TMath::Power(((bData->fTime[i]-t_teo)/bData->fErrors[i]),2);
}
f = (Double_t) f/(bData->fNumHits - 2.);
}
//_____________________________________________________________________________
void TrackData::Clear() {
// TrackData clear
fTime.clear();
fAlpha.clear();
fHits.clear();
fErrors.clear();
fNumPoints = 0;
fNumHits = 0;
fIdMin = 0;
fIdMax = 0;
fNorm.SetXYZ(0,0,0);
fWAxis.SetXYZ(0,0,0);
fUAxis.SetXYZ(0,0,0);
}
//_____________________________________________________________________________
Double_t TimeAtEuso(Double_t *x, Double_t *par){
// par[0] = H_max
// par[1] = Beta
// par[2] = alpha_max (alpha_ref)
// par[3] = t_max
// par[4] = RMax
// x[0] = alpha
Double_t t_teo;
Double_t alpha0 = par[1] + TMath::PiOver2();
Double_t dAlpha = par[2] - alpha0;
Double_t R0 = par[4]*TMath::Sin(par[2] -par[1]);
Double_t t_prop = TMath::Abs(R0/EConst::C()*1/TMath::Cos(dAlpha));
Double_t Dt = -R0/EConst::C()*TMath::Tan(dAlpha);
Double_t tau0 = par[3]-(Dt+t_prop);
/*
cout << Form("_______________________________________________________________________________") << endl;
cout << "dAlpha = " << dAlpha*TMath::RadToDeg() << " deg"<< endl;
cout << "R0 = " << R0/km << " km"<< endl;
cout << "Dt = " << Dt/microsecond << " microsecond"<< endl;
cout << "t_prop = " << t_prop/microsecond<< " microsecond" << endl;
cout << "tau0 " << tau0/microsecond<< " microsecond" << endl;
cout << "par[0] = " << par[0]/km << " km" << endl;
cout << "par[1] = " << par[1]*TMath::RadToDeg() << " deg" << endl;
cout << "par[2] = " << par[2]*TMath::RadToDeg() << " deg" << endl;
cout << "par[3] = " << par[3]/microsecond << " microsecond" << endl;
cout << "par[4] = " << par[4]*TMath::RadToDeg() << " deg" << endl;
cout << "par[5] = " << par[5]/km << " km"<< endl;
cout << "x[0] = " << x[0] << endl;
cout << Form("_______________________________________________________________________________") << endl;
*/
dAlpha = x[0]-alpha0;
t_prop = TMath::Abs(R0/EConst::C()*(1)/TMath::Cos(dAlpha));
Dt = -R0/EConst::C()*TMath::Tan(dAlpha);
t_teo = tau0+Dt+t_prop;
return t_teo;
}