Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TrackDetectorPlaneModule - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: TrackDetectorPlaneModule.cc,v 1.13 2005/10/02 14:18:08 thea 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>

using namespace sou;

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");
    // check if GTUClusteringModule is not called and try to call HoughTransformModule
    if ( gcm == NULL ) {
         gcm = fEv->GetModuleData("HoughTransform");
    }
    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);
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005