source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/reconstruction/modules/shower/fitting/src/TrackDetectorPlaneModule.cc @ 117

Last change on this file since 117 was 117, checked in by moretto, 11 years ago

ESAF version compilable on mac OS
File size: 9.6 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: TrackDetectorPlaneModule.cc 2846 2010-04-18 06:09:05Z fenu $
3// R. Pesce created May, 10 2004
4
5#include "RecoPixelData.hh"
6#include "RecoEvent.hh"
7#include "TrackDetectorPlaneModule.hh"
8#include "LeastSquaresFit.hh"
9#include "MedianFit.hh"
10#include "RecoRootEvent.hh"
11#include "Config.hh"
12
13#include <TVector3.h>
14#include <TObject.h>
15#include <TH2F.h>
16#include <TH1F.h>
17#include <TLine.h>
18#include <TGraph2D.h>
19#include <TDirectory.h>
20#include <TKey.h>
21#include <TMath.h>
22
23using namespace sou;
24
25ClassImp(TDPData)
26ClassImp(TrackDetectorPlaneModule)
27
28//_____________________________________________________________________________
29TrackDetectorPlaneModule::TrackDetectorPlaneModule() : RecoModule("TrackDetectorPlane") {
30    // ctor
31
32}
33
34//_____________________________________________________________________________
35TrackDetectorPlaneModule::~TrackDetectorPlaneModule() {
36    // dtor
37}
38
39
40//_____________________________________________________________________________
41Bool_t TrackDetectorPlaneModule::Init() {
42        ConfigFileParser* pConf = Config::Get()->GetCF("Electronics","MacroCell");
43        fGTU = pConf->GetNum("MacroCell.fGtuTimeLength");
44
45    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
46    return kTRUE;
47}
48
49//_____________________________________________________________________________
50Bool_t TrackDetectorPlaneModule::PreProcess() {
51    return kTRUE;
52}
53
54//_____________________________________________________________________________
55Bool_t TrackDetectorPlaneModule::Process(RecoEvent *ev) {
56    fEv = ev;
57    const RecoModuleData *gcm = fEv->GetModuleData("GTUClustering");
58    // check if GTUClusteringModule is not called and try to call HoughTransformModule
59    if ( gcm == NULL ) {
60         gcm = fEv->GetModuleData("HoughTransform");
61    }
62    if ( gcm==NULL ) {
63        Msg(EsafMsg::Panic) << "GTUClusteringModule is not found." << MsgDispatch;
64    } else {
65        if ( gcm->GetObj("CluPixels") == NULL ) {
66            Msg(EsafMsg::Warning) << "No data in event" << MsgDispatch;     
67            return kFALSE;
68        }
69        fPointsId = *(vector<Int_t>*)gcm->GetObj("CluPixels");
70        if (fPointsId.size() == 0) {
71            fNumPoints = 0;
72            Msg(EsafMsg::Warning) << "No pixels selected: terminated." << MsgDispatch;                 
73            return kFALSE;
74        } else {
75            fNumPoints = fPointsId.size();
76            Msg(EsafMsg::Debug) << "Processing " << fNumPoints << " points" << MsgDispatch;
77        }
78    }
79
80    FillPlaneData();
81    UseVelocities();
82
83    //save data
84    MyData()->Add("NormX",fNorm.X());
85    MyData()->Add("NormY",fNorm.Y());
86    MyData()->Add("NormZ",fNorm.Z());
87    MyData()->Add("WAxisX",fWAxis.X());
88    MyData()->Add("WAxisY",fWAxis.Y());
89    MyData()->Add("WAxisZ",fWAxis.Z());
90    MyData()->Add("UAxisX",fUAxis.X());
91    MyData()->Add("UAxisY",fUAxis.Y());
92    MyData()->Add("UAxisZ",fUAxis.Z());
93    MyData()->Add("XAxisX",fXAxis.X());
94    MyData()->Add("XAxisY",fXAxis.Y());
95    MyData()->Add("XAxisZ",fXAxis.Z());
96    MyData()->Add("YAxisX",fYAxis.X());
97    MyData()->Add("YAxisY",fYAxis.Y());
98    MyData()->Add("YAxisZ",fYAxis.Z());
99    vector<TVector3> *spPoints = &fPlaneData.fSpPos;
100    MyData()->Add("SpPoints",spPoints);
101    MyData()->Add("NumHits",fNumHits);
102   
103    return kTRUE;
104}
105
106//_____________________________________________________________________________
107Bool_t TrackDetectorPlaneModule::PostProcess() {
108    TVector3 centerPos = fEv->GetHeader().GetTrueShowerMaxPos();
109
110    if ( fNumPoints !=0 ) fPlaneErrors.push_back(fTestPlane);
111    if ( fNumPoints !=0 ) fEASCenterX.push_back(centerPos.X());
112    if ( fNumPoints !=0 ) fEASCenterY.push_back(centerPos.Y());
113
114    return kTRUE;
115}
116//_____________________________________________________________________________
117Bool_t TrackDetectorPlaneModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
118     return kTRUE;
119}
120//_____________________________________________________________________________
121Bool_t TrackDetectorPlaneModule::Done() {
122    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
123
124    Int_t numEv = (Int_t)fPlaneErrors.size();
125    if (numEv < 2) return true;
126
127    Float_t *xc, *yc, *pErr, *cpErr;
128    xc = new Float_t[numEv];
129    yc = new Float_t[numEv];
130    pErr = new Float_t[numEv];
131    cpErr = new Float_t[numEv];
132   
133    for( Int_t i=0; i<numEv; i++ ) {
134        xc[i] = fEASCenterX[i];
135        yc[i] = fEASCenterY[i];
136        cpErr[i] = fPlaneErrors[i];
137        pErr[i] = TMath::ACos(fPlaneErrors[i])*TMath::RadToDeg();
138    }
139
140    delete [] xc;
141    delete [] yc;
142    delete [] pErr;
143    delete [] cpErr;
144    return kTRUE;
145}
146
147//_____________________________________________________________________________
148void TrackDetectorPlaneModule::UserMemoryClean() {
149    if (fPlaneData.fSpPos.size()==0) return;
150    vector<TVector3> *spPoints = (vector<TVector3>*)MyData()->GetObj("SpPoints");
151    (*spPoints).clear();
152    MyData()->RemoveObj("SpPoints");
153}
154
155//_____________________________________________________________________________
156void TrackDetectorPlaneModule::FillPlaneData() {
157    fPlaneData.fNumPoints = fNumPoints;
158
159    TH2F *uSphere =  new TH2F("USphere", "Points on the unitary sphere", 100, -0.5, 0.5, 100, -0.5, 0.5);
160    // retrieve data from the event
161    for(Int_t i=0; i<fNumPoints; i++) {
162        RecoPixelData *pix = fEv->GetRecoPixelData(fPointsId[i]);
163        //FIXME
164        Double_t xx = TMath::Sin(pix->GetTheta())*TMath::Cos(pix->GetPhi());
165        Double_t yy = TMath::Sin(pix->GetTheta())*TMath::Sin(pix->GetPhi());
166        Double_t zz = TMath::Cos(pix->GetTheta());
167        Double_t dxx = TMath::Abs(TMath::Cos(pix->GetTheta()))*pix->GetThetaSigma() +
168                      TMath::Abs(TMath::Sin(pix->GetPhi()))*pix->GetPhiSigma();
169        Double_t dyy = TMath::Abs(TMath::Cos(pix->GetTheta()))*pix->GetThetaSigma() +
170                      TMath::Abs(TMath::Cos(pix->GetPhi()))*pix->GetPhiSigma();
171        Double_t dzz = TMath::Abs(TMath::Sin(pix->GetTheta()))*pix->GetThetaSigma();
172
173        TVector3 pos(-xx,-yy,-zz); //FIXME
174        TVector3 err(dxx,dyy,dzz); //FIXME
175       
176        fPlaneData.fSpPos.push_back(pos); 
177        fPlaneData.fSpErr.push_back(err); //FIXME
178        fPlaneData.fTime.push_back((Double_t)fGTU*pix->GetGtu()*ns);
179        fPlaneData.fHits.push_back(pix->GetCounts());
180
181        fPlaneData.fCentroid += pos*pix->GetCounts(); //FIXME
182        uSphere->Fill(pos[0],pos[1]); //FIXME
183    }
184//    uSphere->Write();
185    delete uSphere;
186
187    fNumHits = 0;
188    for( Int_t i=0; i<fNumPoints; i++ )
189        fNumHits += fPlaneData.fHits[i];
190
191    fPlaneData.fNumHits = fNumHits;
192    fPlaneData.fCentroid.SetMag(1.);
193}
194
195//____________________________________________________________________________
196void TrackDetectorPlaneModule::UseVelocities() {
197// points for fitting velocities
198    vector<Double_t> x;
199    vector<Double_t> y;
200    vector<Double_t> t;
201    vector<Double_t> errx;
202    vector<Double_t> erry;
203    TH2F *tProj =  new TH2F("TProj", "Proj on z axis tg Plane", 100, -0.5, 0.5, 100, -0.5, 0.5);
204
205    for( Int_t i=0; i<fNumPoints; i++ ) {
206        RecoPixelData *pix = fEv->GetRecoPixelData(fPointsId[i]);
207        TVector3 dummy = fPlaneData.fSpPos[i];
208        Double_t dummyZ = dummy.Z();
209        dummy.SetZ(dummyZ);
210        dummy *= 1/dummy.CosTheta();
211        tProj->Fill(dummy.X(),dummy.Y());
212        Double_t dx = TMath::Sqrt(
213                      TMath::Power(TMath::Sin(dummy.Phi())*pix->GetPhiSigma(),2.)
214                     +TMath::Power(TMath::Power(1/TMath::Cos(dummy.Theta()),2)*pix->GetThetaSigma(),2.)
215                      );
216        Double_t dy = TMath::Sqrt(
217                      TMath::Power(TMath::Cos(dummy.Phi())*pix->GetPhiSigma(),2.)
218                     +TMath::Power(TMath::Power(1/TMath::Cos(dummy.Theta()),2)*pix->GetThetaSigma(),2.)
219                     );
220        for( Int_t j=0; j<fPlaneData.fHits[i]; j++ ) {
221            x.push_back(dummy.X());
222            y.push_back(dummy.Y());
223            t.push_back(fPlaneData.fTime[i]);
224            errx.push_back(dx);
225            erry.push_back(dy);
226        }
227    }
228//    tProj->Write();
229    delete tProj;
230    // Median Fit
231    // fit in x-t
232    MedianFit *xMf = new MedianFit(fNumHits, t, x);
233    Double_t sX = xMf->GetSlope();
234    Double_t oX = xMf->GetOffset();
235//    Double_t aX = xMf->GetAbsoluteDeviation();
236
237    // fit in x-t
238    MedianFit *yMf = new MedianFit(fNumHits, t, y);
239    Double_t sY = yMf->GetSlope();
240    Double_t oY = yMf->GetOffset();
241    // normal to TDP
242    //fNorm.SetXYZ(vY,-vX,wY*vX-wX*vY); //with least squares
243    fNorm.SetXYZ(sY,-sX,oY*sX-oX*sY); // with median
244    fNorm.SetMag(1.);
245
246    if ( fNorm.Z() < 0 ) fNorm *= -1;
247   
248    // calculate TDP axis
249    // x axis is the centroid direction
250    fXAxis = fPlaneData.fCentroid.Unit();
251    //
252    fYAxis = (fNorm.Cross(fXAxis)).Unit(); 
253
254    fWAxis = (TVector3(0,0,1).Cross(fNorm)).Unit();
255    fUAxis = fNorm.Cross(fWAxis).Unit();
256
257    // save a TLine for cross check
258    TVector3 lVers = fYAxis;
259    TVector3 zVers;
260    zVers.SetXYZ(0.,0.,-1.);
261    TLine *line = new TLine(zVers.X()-lVers.X(),
262            zVers.Y()-lVers.Y(),
263            zVers.X()+lVers.X(),
264            zVers.Y()+lVers.Y());
265    delete line;
266    // MC truth
267    Double_t fTrueTheta = TMath::Pi()-fEv->GetHeader().GetTrueTheta();
268    Double_t fTruePhi = fEv->GetHeader().GetTruePhi();
269    fTrueDir.SetMagThetaPhi(1.,fTrueTheta,fTruePhi);
270    fTrueTheta *= TMath::RadToDeg();
271    fTruePhi *= TMath::RadToDeg();
272    fTestPlane = fTrueDir*fNorm;
273}
274
275//___________________________________________________
276void TDPData::Clear() {
277    // TDPData clear
278    fPos.clear();
279    fErr.clear();
280    fHits.clear();
281    fSpPos.clear();
282    fSpErr.clear();
283    fTime.clear();
284    fNumPoints = 0;
285    fNumHits = 0;
286    fCentroid.SetXYZ(0,0,0);
287}
Note: See TracBrowser for help on using the repository browser.