// ESAF : Euso Simulation and Analysis Framework
// $Id: TrackDetectorPlaneModule.cc,v 1.11 2005/01/18 09:51:09 pesce Exp $
// R. Pesce created May, 10 2004
#include "RecoPixelData.hh"
#include "RecoEvent.hh"
#include "TrackDetectorPlaneModule.hh"
#include "LeastSquaresFit.hh"
#include "MedianFit.hh"
#include "RecoRootEvent.hh"
#include "Config.hh"
#include <TVector3.h>
#include <TObject.h>
#include <TH2F.h>
#include <TH1F.h>
#include <TLine.h>
#include <TGraph2D.h>
#include <TDirectory.h>
#include <TKey.h>
ClassImp(TDPData)
ClassImp(TrackDetectorPlaneModule)
//_____________________________________________________________________________
TrackDetectorPlaneModule::TrackDetectorPlaneModule() : RecoModule("TrackDetectorPlane") {
// ctor
}
//_____________________________________________________________________________
TrackDetectorPlaneModule::~TrackDetectorPlaneModule() {
// dtor
}
//_____________________________________________________________________________
Bool_t TrackDetectorPlaneModule::Init() {
Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDetectorPlaneModule::PreProcess() {
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDetectorPlaneModule::Process(RecoEvent *ev) {
fEv = ev;
const RecoModuleData *gcm = fEv->GetModuleData("GTUClustering");
if ( gcm==NULL ) {
Msg(EsafMsg::Panic) << "GTUClusteringModule is not found." << MsgDispatch;
} else {
if ( gcm->GetObj("CluPixels") == NULL ) {
Msg(EsafMsg::Warning) << "No data in event" << MsgDispatch;
return kFALSE;
}
fPointsId = *(vector<Int_t>*)gcm->GetObj("CluPixels");
if (fPointsId.size() == 0) {
fNumPoints = 0;
Msg(EsafMsg::Warning) << "No pixels selected: terminated." << MsgDispatch;
return kFALSE;
} else {
fNumPoints = fPointsId.size();
Msg(EsafMsg::Debug) << "Processing " << fNumPoints << " points" << MsgDispatch;
}
}
FillPlaneData();
UseVelocities();
//save data
MyData()->Add("NormX",fNorm.X());
MyData()->Add("NormY",fNorm.Y());
MyData()->Add("NormZ",fNorm.Z());
MyData()->Add("WAxisX",fWAxis.X());
MyData()->Add("WAxisY",fWAxis.Y());
MyData()->Add("WAxisZ",fWAxis.Z());
MyData()->Add("UAxisX",fUAxis.X());
MyData()->Add("UAxisY",fUAxis.Y());
MyData()->Add("UAxisZ",fUAxis.Z());
MyData()->Add("XAxisX",fXAxis.X());
MyData()->Add("XAxisY",fXAxis.Y());
MyData()->Add("XAxisZ",fXAxis.Z());
MyData()->Add("YAxisX",fYAxis.X());
MyData()->Add("YAxisY",fYAxis.Y());
MyData()->Add("YAxisZ",fYAxis.Z());
vector<TVector3> *spPoints = &fPlaneData.fSpPos;
MyData()->Add("SpPoints",spPoints);
MyData()->Add("NumHits",fNumHits);
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDetectorPlaneModule::PostProcess() {
TVector3 centerPos = fEv->GetHeader().GetTrueShowerMaxPos();
if ( fNumPoints !=0 ) fPlaneErrors.push_back(fTestPlane);
if ( fNumPoints !=0 ) fEASCenterX.push_back(centerPos.X());
if ( fNumPoints !=0 ) fEASCenterY.push_back(centerPos.Y());
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDetectorPlaneModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
return kTRUE;
}
//_____________________________________________________________________________
Bool_t TrackDetectorPlaneModule::Done() {
Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
Int_t numEv = (Int_t)fPlaneErrors.size();
if (numEv < 2) return true;
Float_t *xc, *yc, *pErr, *cpErr;
xc = new Float_t[numEv];
yc = new Float_t[numEv];
pErr = new Float_t[numEv];
cpErr = new Float_t[numEv];
for( Int_t i=0; i<numEv; i++ ) {
xc[i] = fEASCenterX[i];
yc[i] = fEASCenterY[i];
cpErr[i] = fPlaneErrors[i];
pErr[i] = TMath::ACos(fPlaneErrors[i])*TMath::RadToDeg();
}
delete [] xc;
delete [] yc;
delete [] pErr;
delete [] cpErr;
return kTRUE;
}
//_____________________________________________________________________________
void TrackDetectorPlaneModule::UserMemoryClean() {
if (fPlaneData.fSpPos.size()==0) return;
vector<TVector3> *spPoints = (vector<TVector3>*)MyData()->GetObj("SpPoints");
(*spPoints).clear();
MyData()->RemoveObj("SpPoints");
}
//_____________________________________________________________________________
void TrackDetectorPlaneModule::FillPlaneData() {
fPlaneData.fNumPoints = fNumPoints;
TH2F *uSphere = new TH2F("USphere", "Points on the unitary sphere", 100, -0.5, 0.5, 100, -0.5, 0.5);
// retrieve data from the event
for(Int_t i=0; i<fNumPoints; i++) {
RecoPixelData *pix = fEv->GetRecoPixelData(fPointsId[i]);
//FIXME
Double_t xx = TMath::Sin(pix->GetTheta())*TMath::Cos(pix->GetPhi());
Double_t yy = TMath::Sin(pix->GetTheta())*TMath::Sin(pix->GetPhi());
Double_t zz = TMath::Cos(pix->GetTheta());
Double_t dxx = TMath::Abs(TMath::Cos(pix->GetTheta()))*pix->GetThetaSigma() +
TMath::Abs(TMath::Sin(pix->GetPhi()))*pix->GetPhiSigma();
Double_t dyy = TMath::Abs(TMath::Cos(pix->GetTheta()))*pix->GetThetaSigma() +
TMath::Abs(TMath::Cos(pix->GetPhi()))*pix->GetPhiSigma();
Double_t dzz = TMath::Abs(TMath::Sin(pix->GetTheta()))*pix->GetThetaSigma();
TVector3 pos(-xx,-yy,-zz); //FIXME
TVector3 err(dxx,dyy,dzz); //FIXME
fPlaneData.fSpPos.push_back(pos);
fPlaneData.fSpErr.push_back(err); //FIXME
fPlaneData.fTime.push_back((Double_t)2500.*pix->GetGtu()*ns);
fPlaneData.fHits.push_back(pix->GetCounts());
fPlaneData.fCentroid += pos*pix->GetCounts(); //FIXME
uSphere->Fill(pos[0],pos[1]); //FIXME
}
// uSphere->Write();
delete uSphere;
fNumHits = 0;
for( Int_t i=0; i<fNumPoints; i++ )
fNumHits += fPlaneData.fHits[i];
fPlaneData.fNumHits = fNumHits;
fPlaneData.fCentroid.SetMag(1.);
}
//____________________________________________________________________________
void TrackDetectorPlaneModule::UseVelocities() {
// points for fitting velocities
vector<Double_t> x;
vector<Double_t> y;
vector<Double_t> t;
vector<Double_t> errx;
vector<Double_t> erry;
TH2F *tProj = new TH2F("TProj", "Proj on z axis tg Plane", 100, -0.5, 0.5, 100, -0.5, 0.5);
for( Int_t i=0; i<fNumPoints; i++ ) {
RecoPixelData *pix = fEv->GetRecoPixelData(fPointsId[i]);
TVector3 dummy = fPlaneData.fSpPos[i];
Double_t dummyZ = dummy.Z();
dummy.SetZ(dummyZ);
dummy *= 1/dummy.CosTheta();
tProj->Fill(dummy.X(),dummy.Y());
Double_t dx = TMath::Sqrt(
TMath::Power(TMath::Sin(dummy.Phi())*pix->GetPhiSigma(),2.)
+TMath::Power(TMath::Power(1/TMath::Cos(dummy.Theta()),2)*pix->GetThetaSigma(),2.)
);
Double_t dy = TMath::Sqrt(
TMath::Power(TMath::Cos(dummy.Phi())*pix->GetPhiSigma(),2.)
+TMath::Power(TMath::Power(1/TMath::Cos(dummy.Theta()),2)*pix->GetThetaSigma(),2.)
);
for( Int_t j=0; j<fPlaneData.fHits[i]; j++ ) {
x.push_back(dummy.X());
y.push_back(dummy.Y());
t.push_back(fPlaneData.fTime[i]);
errx.push_back(dx);
erry.push_back(dy);
}
}
// tProj->Write();
delete tProj;
// Median Fit
// fit in x-t
MedianFit *xMf = new MedianFit(fNumHits, t, x);
Double_t sX = xMf->GetSlope();
Double_t oX = xMf->GetOffset();
// Double_t aX = xMf->GetAbsoluteDeviation();
// fit in x-t
MedianFit *yMf = new MedianFit(fNumHits, t, y);
Double_t sY = yMf->GetSlope();
Double_t oY = yMf->GetOffset();
// normal to TDP
//fNorm.SetXYZ(vY,-vX,wY*vX-wX*vY); //with least squares
fNorm.SetXYZ(sY,-sX,oY*sX-oX*sY); // with median
fNorm.SetMag(1.);
if ( fNorm.Z() < 0 ) fNorm *= -1;
// calculate TDP axis
// x axis is the centroid direction
fXAxis = fPlaneData.fCentroid.Unit();
//
fYAxis = (fNorm.Cross(fXAxis)).Unit();
fWAxis = (TVector3(0,0,1).Cross(fNorm)).Unit();
fUAxis = fNorm.Cross(fWAxis).Unit();
// save a TLine for cross check
TVector3 lVers = fYAxis;
TVector3 zVers;
zVers.SetXYZ(0.,0.,-1.);
TLine *line = new TLine(zVers.X()-lVers.X(),
zVers.Y()-lVers.Y(),
zVers.X()+lVers.X(),
zVers.Y()+lVers.Y());
delete line;
// MC truth
Double_t fTrueTheta = TMath::Pi()-fEv->GetHeader().GetTrueTheta();
Double_t fTruePhi = fEv->GetHeader().GetTruePhi();
fTrueDir.SetMagThetaPhi(1.,fTrueTheta,fTruePhi);
fTrueTheta *= TMath::RadToDeg();
fTruePhi *= TMath::RadToDeg();
fTestPlane = fTrueDir*fNorm;
}
//___________________________________________________
void TDPData::Clear() {
// TDPData clear
fPos.clear();
fErr.clear();
fHits.clear();
fSpPos.clear();
fSpErr.clear();
fTime.clear();
fNumPoints = 0;
fNumHits = 0;
fCentroid.SetXYZ(0,0,0);
}