source: JEM-EUSO/esaf_cc_at_lal/packages/reconstruction/modules/base/clustering/src/GTUClusteringModule.cc @ 114

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

actual version of ESAF at CCin2p3
File size: 14.6 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: GTUClusteringModule.cc 2929 2011-06-16 20:50:44Z mabl $
3// R. Pesce created Mar, 11 2004
4
5#include "GTUClusteringModule.hh"
6#include "RecoEvent.hh"
7#include "RecoPixelData.hh"
8#include "EusoCluster.hh"
9#include "RecoRootEvent.hh"
10
11
12#include <TDirectory.h>
13#include <map>
14#include <vector>
15#include <TMath.h>
16#include <TH3F.h>
17
18ClassImp(GTUClusteringModule)
19//_____________________________________________________________________________
20// ctor
21GTUClusteringModule::GTUClusteringModule() : RecoModule("GTUClustering"){
22}
23
24//_____________________________________________________________________________
25// dtor
26GTUClusteringModule::~GTUClusteringModule() {
27}
28
29//_____________________________________________________________________________
30// init
31Bool_t GTUClusteringModule::Init() {
32    fEv = (RecoEvent*)0;
33    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
34    Msg(EsafMsg::Debug) << "with GTUClusteringModule.NaturalValues = " << Conf()->GetNum("GTUClusteringModule.NaturalValues")<< MsgDispatch;
35
36    fNumHitsMinimum = (Int_t)Conf()->GetNum("GTUClusteringModule.NumHitsMinimum");
37    if ( fNumHitsMinimum <= 0 ) {
38        Msg(EsafMsg::Panic) << "You must specify a positive number in GTUClusteringModule.NumHitsMinimum" << MsgDispatch;
39    }
40    if ((Int_t)Conf()->GetNum("GTUClusteringModule.NaturalValues")==1) {
41        fSignificanceLevel = (Int_t)Conf()->GetNum("GTUClusteringModule.SignificanceLevel");
42        if ( fSignificanceLevel <= 0 ) {
43            Msg(EsafMsg::Panic) << "You must specify a positive number in GTUClusteringModule.SignificanceLevel" << MsgDispatch;
44        }
45        fNumPointsMinimum = 0;
46        fUseNatural = kTRUE;
47    } else if ((Int_t)Conf()->GetNum("GTUClusteringModule.NaturalValues")==0) {
48        fNumPointsMinimum = (Int_t)Conf()->GetNum("GTUClusteringModule.NumPointsMinimum");
49        if ( fNumPointsMinimum <= 0 ) {
50            Msg(EsafMsg::Panic) << "You must specify a positive number in GTUClusteringModule.NumPointsMinimum" << MsgDispatch;
51         }
52        fSignificanceLevel = 0;
53        fUseNatural = kFALSE;
54        fThreshold = (Double_t)Conf()->GetNum("GTUClusteringModule.DistanceThreshold");
55        if ( fThreshold <= 0 ) {
56            Msg(EsafMsg::Panic) << "You must specify a positive number in GTUClusteringModule.DistanceThreshold" << MsgDispatch;
57        }
58    } else {
59        Msg(EsafMsg::Panic) << "You must specify 0 or 1 in GTUClusteringModule.NaturalValues" << MsgDispatch;
60    }
61   
62    return kTRUE;
63}
64//_____________________________________________________________________________
65// pre-process
66Bool_t GTUClusteringModule::PreProcess() {
67    fSpClusters.clear();
68    fPixels.clear();
69    fCluPixels.clear();
70    fId.clear();
71    fFlag1.clear();
72    fFlag2.clear();
73    fNumPoints = 0;
74    if ( fUseNatural==kTRUE ) {
75        fNumPointsMinimum = 0;
76        fThreshold = 0;
77        fDensity = 0;
78    }
79       
80    return kTRUE;
81}
82
83//_____________________________________________________________________________
84// process
85Bool_t GTUClusteringModule::Process( RecoEvent* ev ) {
86    fEv = ev;
87    const RecoModuleData *camd = fEv->GetModuleData("ClusterAnalysis");
88    if ( camd == NULL ) {
89        Msg(EsafMsg::Panic) << "ClusterAnalysisModule is not found." << MsgDispatch;
90    } else {
91        if ( camd->GetObj("Clusters") == NULL ) {
92            Msg(EsafMsg::Panic) << "No clusters to process." << MsgDispatch;
93            return kFALSE;
94        }
95        fSpClusters = *(vector<EusoCluster*>*) camd->GetObj("Clusters");
96        Msg(EsafMsg::Debug) << "ClusterAnalysisModule data found." << MsgDispatch;
97    }
98    Int_t hitsmin = camd->GetInt("NumHitsMinimum");
99    // search pixels in angular range with num hits minimum
100    if (hitsmin==1) {
101        for( Int_t i=0; i<(Int_t)fSpClusters.size(); i++ ) {
102            for( Int_t j=0; j<fSpClusters[i]->GetNumPoints(); j++ ) {
103                fPixels.push_back(fSpClusters[i]->GetPixelId(j));
104            }
105        }
106    } else {
107        for( Int_t i=0; i<(Int_t)fSpClusters.size(); i++ ) {
108            Double_t thmin = fSpClusters[i]->GetThetaMin();
109            Double_t thmax = fSpClusters[i]->GetThetaMax();
110            Double_t phimin = fSpClusters[i]->GetPhiMin();
111            Double_t phimax = fSpClusters[i]->GetPhiMax();
112            Double_t tmin = fSpClusters[i]->GetGtuMin();
113            Double_t tmax = fSpClusters[i]->GetGtuMax();
114            for( Int_t j=0; j<fEv->GetHeader().GetNumActiveFee(); j++ ) {
115                RecoPixelData *pix = fEv->GetRecoPixelData(j);
116                if ( pix->GetTheta() >= thmin && pix->GetTheta() <= thmax &&
117                    pix->GetPhi() >= phimin && pix->GetPhi() <= phimax &&
118                    pix->GetGtu() >= tmin && pix->GetGtu() <= tmax &&
119                    pix->GetCounts() >= fNumHitsMinimum )
120                    fPixels.push_back(j);
121            }
122
123        }
124    }
125    fNumPoints = fPixels.size();
126
127    Msg(EsafMsg::Debug) << "Pixels found = " << fNumPoints << MsgDispatch;
128
129    if ( fUseNatural == kTRUE ) {
130        CalculateDensity();
131        CalculateThreshold();
132        CalculateNumClustersUniform();
133        CalculateNumPointsMinimum();
134        //cout << "##### " << fDensity << " " << fThreshold <<" " << fNumClustersUniform
135             //<< " " << fNumPointsMinimum << endl;
136    }
137    Prepare();
138    SearchClusters();
139
140    Msg(EsafMsg::Debug) << "found " << fTmClusters.size() << " clusters." << MsgDispatch;
141    if ( fTmClusters.size() == 0 ) return kFALSE;
142
143    vector<Int_t> medTime;
144    vector<Int_t> numPC;
145    for( size_t i(0); i<fTmClusters.size(); i++ ) {
146        Int_t sum = 0;
147        for( size_t j(0); j<(fTmClusters[i]).size(); j++ ) {
148            sum += fEv->GetRecoPixelData((fTmClusters[i])[j])->GetGtu();
149        }
150        medTime.push_back(sum/(fTmClusters[i]).size());
151        numPC.push_back((fTmClusters[i]).size());
152    }
153    Int_t cluMostPop(0);
154    for( size_t i(0); i<numPC.size(); i++ ) {
155        if ( numPC[i] > numPC[cluMostPop] ) cluMostPop = i;
156    }
157    Int_t timeMP = medTime[cluMostPop];
158
159    vector<Int_t> selTClu;
160    selTClu.push_back(cluMostPop);
161
162    for( size_t i(0); i<medTime.size(); i++ ) {
163        if( TMath::Abs(timeMP-medTime[i]) < (size_t)2*fThreshold && i!=(size_t)cluMostPop ) {
164            selTClu.push_back(i);
165        }
166    }
167    Msg(EsafMsg::Debug) << "Selected " << selTClu.size() << " clusters." << MsgDispatch;       
168    vector<Int_t> tclu;
169    for( size_t i=0; i<selTClu.size(); i++ ) {
170        tclu = fTmClusters[selTClu[i]];
171        for(Int_t j=0; j<(Int_t)tclu.size(); j++) {
172            fCluPixels.push_back(tclu[j]);
173        }
174    }
175
176    return kTRUE;
177}
178
179//_____________________________________________________________________________
180// post-process
181Bool_t GTUClusteringModule::PostProcess() {
182    if ( fPixels.size() == 0 ) return kTRUE;
183    if ((Int_t)Conf()->GetNum("GTUClusteringModule.DoHistogram") == 1) 
184        Histo();
185
186    vector< vector<Int_t> > *pclu = &fTmClusters;
187    vector<Int_t> *pix = &fCluPixels;
188    MyData()->Add("TmClusters", pclu);
189    MyData()->Add("CluPixels", pix);
190   
191    fEv = (RecoEvent*)0;
192   
193    return kTRUE;
194}
195
196//_____________________________________________________________________________
197Bool_t GTUClusteringModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
198     return kTRUE;
199}
200
201//_____________________________________________________________________________
202// done
203Bool_t GTUClusteringModule::Done() {
204    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
205    return kTRUE;
206}
207
208//_____________________________________________________________________________
209// memory clean
210void GTUClusteringModule::UserMemoryClean() {
211    if ( fPixels.size() == 0 ) return;
212
213    if ( fTmClusters.size() != 0 ) {
214        vector< vector<Int_t> > *pclu = (vector< vector<Int_t> >*)MyData()->GetObj("TmClusters");
215        if ( pclu != 0 && (*pclu).size() != 0 ) {
216            (*pclu).clear();
217            MyData()->RemoveObj("TmClusters");
218        }
219    }   
220   
221    if ( fCluPixels.size() !=0 ) {
222        vector<Int_t> *pix = (vector<Int_t>*)MyData()->GetObj("CluPixels");
223        if ( pix != NULL && (*pix).size() != 0 ) {
224            (*pix).clear();
225            MyData()->RemoveObj("CluPixels");
226        }
227    }
228}
229
230// density: is simply the num of pixel minimum in a gtu
231void GTUClusteringModule::CalculateDensity() {
232    Int_t tmin = 0;
233    Int_t tmax = 0;
234    for( Int_t i=0; i<(Int_t)fPixels.size(); i++ ) {
235        Int_t gtuid = fEv->GetRecoPixelData(fPixels[i])->GetGtu();
236        if ( gtuid < tmin ) tmin = gtuid;
237        if ( gtuid > tmax ) tmax = gtuid;
238    }
239    fDensity = fPixels.size()/(tmax-tmin+1);
240}
241
242// gtu distance
243Int_t GTUClusteringModule::GTUDistance(Int_t n1, Int_t n2) {
244    return TMath::Abs( fEv->GetRecoPixelData(fPixels[n1])->GetGtu() - 
245                       fEv->GetRecoPixelData(fPixels[n2])->GetGtu() );
246}
247
248// natural threshold
249void GTUClusteringModule::CalculateThreshold() {
250    fThreshold = TMath::Log(2)/fDensity;
251}
252
253// calculate number of clusters expected in a distribution uniform
254void GTUClusteringModule::CalculateNumClustersUniform() {
255    fNumClustersUniform = (Int_t)( 1 + ( fNumPoints - 1 ) *
256        TMath::Exp( - fDensity * fThreshold ) );
257}
258
259// calculate the minimum number of points in a significant cluster
260void GTUClusteringModule::CalculateNumPointsMinimum() {
261    fNumPointsMinimum = (Int_t) ( fNumPoints / fNumClustersUniform
262        + fSignificanceLevel * TMath::Sqrt( (Double_t)( fNumPoints / fNumClustersUniform ) ) );
263}
264
265// preparation for clustering
266void GTUClusteringModule::Prepare() {
267    for( Int_t i=0; i<fNumPoints; i++ ) {
268        fId.push_back(0);
269        fFlag1.push_back(kFALSE);
270        fFlag2.push_back(kFALSE);
271    }
272    fNumClusters = 0;
273    fNumNodes = 0;
274    fNumSigClusters = 0;
275    fBookmark = 0;
276    fTmClusters.clear();
277}
278
279// search for clusters
280void GTUClusteringModule::SearchClusters() {
281    while( fNumNodes < fNumPoints ) {
282        NewCluster();
283        FillCluster();
284        WriteCluster();
285    }
286}
287
288// init a new cluster
289void GTUClusteringModule::NewCluster() {
290    fNumPointsCluster = 1;
291    fNumClusters++;
292    fNumNodes++;
293    fFlag2[fBookmark] = kTRUE;
294    fId[0] = fBookmark;
295    fCounter = 0;
296}
297
298// fill a cluster
299void GTUClusteringModule::FillCluster() {
300    while ( fCounter < fNumPointsCluster ) {
301   
302        Int_t k = fId[fCounter];
303        if ( !fFlag1[k] ) {
304            fFlag1[k] = kTRUE;
305            for (Int_t i=0; i<fNumPoints; i++) {
306                if ( !fFlag2[i] ) {
307                    if ( GTUDistance(i, k) <= fThreshold ) {
308                        fFlag2[i] = kTRUE;
309                        fNumPointsCluster++;
310                        fId[fNumPointsCluster-1] = i;
311                        fNumNodes++;
312                    }
313                    fBookmark = i;
314                }
315            }
316        }
317        fCounter++;
318    }
319}
320
321// write a cluster
322void GTUClusteringModule::WriteCluster() {
323    // check cluster significance
324    if ( fNumPointsCluster < fNumPointsMinimum )
325        return;
326   
327    vector<Int_t> clu;
328    for( Int_t i=0; i<fNumPointsCluster; i++ ) {
329        clu.push_back(fPixels[fId[i]]);
330    }
331    fTmClusters.push_back(clu);
332}
333
334// histogram
335void GTUClusteringModule::Histo() {
336
337    gDirectory->cd("/");
338    string pathname = Form("timecluster%d", fEv->GetHeader().GetNum());
339    TDirectory *path = new TDirectory(pathname.c_str(),pathname.c_str());
340    path->cd();
341
342    Float_t thmin = fEv->GetRecoPixelData(0)->GetTheta();
343    Float_t thmax = fEv->GetRecoPixelData(0)->GetTheta();
344    Float_t phimin = fEv->GetRecoPixelData(0)->GetPhi();
345    Float_t phimax = fEv->GetRecoPixelData(0)->GetPhi();
346    Float_t tmin = fEv->GetRecoPixelData(0)->GetGtu();
347    Float_t tmax = fEv->GetRecoPixelData(0)->GetGtu();
348
349    for( Int_t i=1; i<(Int_t)fEv->GetHeader().GetNumActiveFee(); i++ ) {
350        if ( fEv->GetRecoPixelData(i)->GetTheta() < thmin )
351            thmin = fEv->GetRecoPixelData(i)->GetTheta();
352        if ( fEv->GetRecoPixelData(i)->GetTheta() > thmax )
353            thmax = fEv->GetRecoPixelData(i)->GetTheta();
354        if ( fEv->GetRecoPixelData(i)->GetPhi() < phimin )
355            phimin = fEv->GetRecoPixelData(i)->GetPhi();
356        if ( fEv->GetRecoPixelData(i)->GetPhi() > phimax )
357            phimax = fEv->GetRecoPixelData(i)->GetPhi();
358        if ( fEv->GetRecoPixelData(i)->GetGtu() < tmin )
359            tmin = fEv->GetRecoPixelData(i)->GetGtu();
360        if ( fEv->GetRecoPixelData(i)->GetGtu() > tmax )
361            tmax = fEv->GetRecoPixelData(i)->GetGtu();
362}
363
364    Int_t bin = (Int_t)Conf()->GetNum("GTUClusteringModule.HistogramBinning");;
365    TH3F *histo = new TH3F("Histo","Points",bin,phimin,phimax,bin,thmin,thmax,bin,tmin,tmax);
366   
367    for( Int_t i=0; i<(Int_t)fEv->GetHeader().GetNumActiveFee(); i++ ) {
368        histo->Fill( fEv->GetRecoPixelData(i)->GetPhi(),fEv->GetRecoPixelData(i)->GetTheta(),
369                     fEv->GetRecoPixelData(i)->GetGtu() );
370    }
371    histo->GetXaxis()->SetTitle("#varphi (rad)");
372    histo->GetYaxis()->SetTitle("#theta (rad)");
373    histo->GetZaxis()->SetTitle("GTU");
374
375    TH3F *tmhisto = new TH3F("TmHisto","Clustered Points",bin,phimin,phimax,bin,thmin,thmax,
376                              bin,tmin,tmax);
377    thmin = fEv->GetRecoPixelData(fCluPixels[0])->GetTheta();
378    thmax = fEv->GetRecoPixelData(fCluPixels[0])->GetTheta();
379    phimin = fEv->GetRecoPixelData(fCluPixels[0])->GetPhi();
380    phimax = fEv->GetRecoPixelData(fCluPixels[0])->GetPhi();
381    tmin = fEv->GetRecoPixelData(fCluPixels[0])->GetGtu();
382    tmax = fEv->GetRecoPixelData(fCluPixels[0])->GetGtu();
383   
384    for(Int_t i=0; i<(Int_t)fCluPixels.size(); i++) {
385        RecoPixelData *pix = fEv->GetRecoPixelData(fCluPixels[i]);
386        if (pix->GetTheta() < thmin ) thmin = pix->GetTheta();
387        if (pix->GetTheta() > thmax ) thmax = pix->GetTheta();
388        if (pix->GetPhi() < phimin ) phimin = pix->GetPhi();
389        if (pix->GetPhi() > phimax ) phimax = pix->GetPhi();
390        if (pix->GetGtu() < tmin ) tmin = pix->GetGtu();
391        if (pix->GetGtu() > tmax ) tmax = pix->GetGtu();
392        tmhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
393    }
394    tmhisto->GetXaxis()->SetTitle("#varphi (rad)");
395    tmhisto->GetYaxis()->SetTitle("#theta (rad)");
396    tmhisto->GetZaxis()->SetTitle("GTU");
397
398    TH3F *zoomtmhisto = new TH3F("ZTmHisto","Clustered points",bin,phimin,phimax,bin,thmin,thmax,
399                                  bin,tmin,tmax);
400    for(Int_t i=0; i<(Int_t)fCluPixels.size(); i++) {
401        RecoPixelData *pix = fEv->GetRecoPixelData(fCluPixels[i]);
402        zoomtmhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
403    }
404    zoomtmhisto->GetXaxis()->SetTitle("#varphi (rad)");
405    zoomtmhisto->GetYaxis()->SetTitle("#theta (rad)");
406    zoomtmhisto->GetZaxis()->SetTitle("GTU");
407   
408    histo->Write();
409    tmhisto->Write();
410    zoomtmhisto->Write();
411    delete histo;
412    delete tmhisto;
413    delete zoomtmhisto;
414}
Note: See TracBrowser for help on using the repository browser.