source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/reconstruction/modules/base/clustering/src/ClusterAnalysisModule.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: 18.0 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: ClusterAnalysisModule.cc 2949 2011-06-26 18:39:34Z mabl $
3// R. Pesce created Mar,  4 2004
4
5#include "TMath.h"
6#include "RecoFramework.hh"
7#include "ClusterAnalysisModule.hh"
8#include "EusoCluster.hh"
9#include "AutomatedClustering.hh"
10#include "RecoEvent.hh"
11#include "RecoPixelData.hh"
12#include "RecoRootEvent.hh"
13
14#include <TH3F.h>
15#include <TH2F.h>
16#include <TDirectory.h>
17
18ClassImp(ClusterAnalysisModule)
19
20//_____________________________________________________________________________
21// ctor
22ClusterAnalysisModule::ClusterAnalysisModule() : RecoModule("ClusterAnalysis") {
23}
24
25//_____________________________________________________________________________
26// dtor
27ClusterAnalysisModule::~ClusterAnalysisModule() {
28}
29
30//_____________________________________________________________________________
31// get a specified cluster
32EusoCluster* ClusterAnalysisModule::GetCluster( Int_t i ) {
33    if ( i >= 0 && i < (Int_t)fClusters.size() )
34        return fClusters[i];
35    else
36        return NULL;
37}
38
39//_____________________________________________________________________________
40// init
41Bool_t ClusterAnalysisModule::Init() {
42    fEv = (RecoEvent*)0;
43    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
44    fNumThreshold = (Double_t)Conf()->GetNum("ClusterAnalysisModule.NumThreshold");
45    if ( fNumThreshold <= 0 ) {
46        Msg(EsafMsg::Panic) << "You must specify a positive number in ClusterAnalysisModule.NumThreshold" << MsgDispatch;
47    }
48
49    fRelNumPointsMin = (Double_t)Conf()->GetNum("ClusterAnalysisModule.RelNumPointsMin");
50    if ( fRelNumPointsMin <= 0 ) {
51        Msg(EsafMsg::Panic) << "You must specify a positive number in ClusterAnalysisModule.RelNumPointsMin" << MsgDispatch;
52    }
53   
54    fRelNumHits = (Double_t)Conf()->GetNum("ClusterAnalysisModule.RelNumHits");
55    if ( fRelNumHits <= 0 ) {
56        Msg(EsafMsg::Panic) << "You must specify a positive number in ClusterAnalysisModule.RelNumHits" << MsgDispatch;
57    }
58
59    gDoFit = (Bool_t)Conf()->GetNum("ClusterAnalysisModule.DoMedianFit");
60    return true;
61}
62
63//_____________________________________________________________________________
64// pre-process
65Bool_t ClusterAnalysisModule::PreProcess() {
66    fMainCluster = 0;
67    fThetaMean.clear();
68    fPhiMean.clear();
69    fClusters.clear();
70    fSelectedClusters.clear();
71    fId.clear();
72    fFlag1.clear();
73    fFlag2.clear();
74    return true;
75}
76
77//_____________________________________________________________________________
78// process
79Bool_t ClusterAnalysisModule::Process( RecoEvent *ev ) {
80    fEv = ev;
81    Int_t siglevel(0), hitsmin(0), pointsmin(0);
82    const RecoModuleData *bcmd = fEv->GetModuleData("BaseClustering");
83    if ( bcmd == NULL ) {
84        Msg(EsafMsg::Panic) << "BaseClusteringModule is not found" << MsgDispatch;
85    } else {
86        if ( bcmd->GetObj("Clusters") == NULL ) return kFALSE;
87        siglevel = bcmd->GetInt("SignificanceLevel");
88        hitsmin = bcmd->GetInt("NumHitsMinimum");
89        fThreshold = bcmd->GetDouble("DistanceThreshold");
90        pointsmin = bcmd->GetInt("NumPointsMinimum");
91        fClusters = *(vector<EusoCluster*>*) bcmd->GetObj("Clusters");
92        Msg(EsafMsg::Debug) << "BaseClusteringModule data found." << MsgDispatch;
93    }
94   
95    if ( fClusters.size() == 0 ) {
96        Msg(EsafMsg::Panic) << "There are no significant clusters. " << fSelectedClusters.size() << MsgDispatch;
97        return kFALSE;
98    } else if ( fClusters.size() == 1 ) {
99        Msg(EsafMsg::Debug) << "There is only one significant cluster." << MsgDispatch;
100        fSelectedClusters.push_back(0);
101        MyData()->Add("NumHitsMinimum", hitsmin);
102        return kTRUE;
103    } else {
104        Msg(EsafMsg::Debug)<< "There are " << fClusters.size() << " significant clusters." << MsgDispatch;
105    }
106   
107
108    if(Conf()->GetStr("ClusterAnalysisModule.ReClustering") == "yes") {
109        fClusters.clear();
110        pointsmin = (Int_t) (pointsmin*fRelNumPointsMin);
111        hitsmin = (Int_t) (hitsmin*fRelNumHits);
112        AutomatedClustering *ac = 
113            new AutomatedClustering( fEv, hitsmin, fThreshold*fNumThreshold, pointsmin, gDoFit );
114        fClusters = ac->GetClusters();
115        Msg(EsafMsg::Debug)<< "ReClustering: found " << fClusters.size() << " significant clusters." << MsgDispatch;
116        for( Int_t i=0; i<(Int_t)fClusters.size(); i++ ) {
117            Msg(EsafMsg::Debug) << "Cluster n. " << i << ": Number of points = " << fClusters[i]->GetNumPoints() << MsgDispatch;
118            /*if (gDoFit) {
119                 Msg(EsafMsg::Debug)<< " Slope = " << fClusters[i]->GetSlope() << " Offset = "
120                     << fClusters[i]->GetOffset() << " AbsoluteDeviation = "
121                     << fClusters[i]->GetAbsoluteDeviation() << MsgDispatch;
122            }*/
123        }
124    } else if(Conf()->GetStr("ClusterAnalysisModule.ReClustering") == "no") {
125        Msg(EsafMsg::Debug) << "No ReClustering." << MsgDispatch;
126    } else {
127        Msg(EsafMsg::Panic) << "You must specify yes or no in ClusterAnalysisModule.RiClustering" << MsgDispatch; 
128    }
129
130    SearchMainCluster();
131   
132    for(Int_t i=0; i<(Int_t)fClusters.size(); i++) {
133        Double_t thmin = fClusters[i]->GetThetaMin();
134        Double_t thmax = fClusters[i]->GetThetaMax();
135        Double_t phimin = fClusters[i]->GetPhiMin();
136        Double_t phimax = fClusters[i]->GetPhiMax();
137        fThetaMean.push_back( (thmin + thmax) / 2 );
138        fPhiMean.push_back( (phimin + phimax) / 2 );
139        Double_t sum = AngularDistance(fThetaMean[i], fPhiMean[i], thmin, phimin);
140        sum += AngularDistance(fThetaMean[i], fPhiMean[i], thmin, phimax); 
141        sum += AngularDistance(fThetaMean[i], fPhiMean[i], thmax, phimin); 
142        sum += AngularDistance(fThetaMean[i], fPhiMean[i], thmax, phimax);
143        fDim.push_back(sum / 4.);
144    }
145
146    InitSuperClustering();
147    SearchSuperClusters();
148
149    if ( fSelectedClusters.size() == 0 ) {
150        fSelectedClusters.push_back(fMainCluster);
151        Msg(EsafMsg::Debug) << "Selected main cluster." << MsgDispatch;
152    } else {
153        Msg(EsafMsg::Debug) << "Selected " << fSelectedClusters.size() << " clusters" << MsgDispatch;
154        Int_t pointsum = 0;
155        Bool_t mainselect = kFALSE;
156        for(Int_t i=0; i<(Int_t)fSelectedClusters.size(); i++) {
157            pointsum += fClusters[fSelectedClusters[i]]->GetNumPoints();
158            if (fSelectedClusters[i] == fMainCluster) mainselect = kTRUE;
159        }
160        Msg(EsafMsg::Debug) << "Number of points clusterized = " << pointsum << MsgDispatch;
161        if (mainselect==kFALSE) {
162            Msg(EsafMsg::Warning) << "Main Cluster is not selected." << MsgDispatch;
163            Double_t thmean(0), phimean(0), dim(0);
164            for(Int_t i=0; i<(Int_t)fSelectedClusters.size(); i++) {
165                thmean += fThetaMean[fSelectedClusters[i]];
166                phimean += fPhiMean[fSelectedClusters[i]];
167                dim += fDim[fSelectedClusters[i]]+fThreshold;
168            }
169            Double_t dist = AngularDistance(thmean, phimean, fThetaMean[fMainCluster],
170                            fPhiMean[fMainCluster]);
171            dist -= (fDim[fMainCluster]+dim);
172            if (dist <= 2.*fThreshold) {
173                Msg(EsafMsg::Debug) << "Checking: Main Cluster selected" << MsgDispatch;
174                pointsum += fClusters[fMainCluster]->GetNumPoints();
175                Msg(EsafMsg::Debug) << "Total points clusterized = " << pointsum << MsgDispatch;
176                fSelectedClusters.push_back(fMainCluster);
177            } else {
178                if (fClusters[fMainCluster]->GetNumPoints() > 2*pointsum) {
179                    fSelectedClusters.clear();
180                    fSelectedClusters.push_back(fMainCluster);
181                    Msg(EsafMsg::Debug) << "Kept Main Cluster." << MsgDispatch;
182                } else if (pointsum > 2*fClusters[fMainCluster]->GetNumPoints()) {
183                    Msg(EsafMsg::Debug) << "Kept selected clusters." << MsgDispatch;
184                } else {
185                    fSelectedClusters.push_back(fMainCluster);
186                    Msg(EsafMsg::Debug) << "Doubtful case: kept selected clusters and main cluster" << MsgDispatch;
187                }
188            }
189        }
190    }
191    MyData()->Add("NumHitsMinimum", hitsmin);
192
193    return kTRUE;
194}
195
196//_____________________________________________________________________________
197// post-process
198Bool_t ClusterAnalysisModule::PostProcess() {
199    if (fSelectedClusters.size() == 0) return true;
200
201    if ((Bool_t)Conf()->GetNum("ClusterAnalysisModule.DoHistogram")) HistoClusters();
202   
203    // keep only selected clusters
204    vector<EusoCluster*> temp;
205    EusoCluster *dummy;
206    for(Int_t i=0; i<(Int_t)fSelectedClusters.size(); i++) {
207        dummy = new EusoCluster();
208        (*dummy) = *(fClusters[fSelectedClusters[i]]);
209        temp.push_back(dummy);
210    }
211    fClusters.clear();
212    fClusters = temp;
213    temp.clear();
214
215 
216    // add data to event
217    vector<EusoCluster*> *pclu = &fClusters;
218    MyData()->Add("Clusters", pclu);
219
220    fEv = (RecoEvent*)0;
221
222    return true;
223}
224
225Bool_t ClusterAnalysisModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
226     return kTRUE;
227}
228// done
229Bool_t ClusterAnalysisModule::Done() {
230   Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
231   return true;
232}
233
234// user memory clean
235void ClusterAnalysisModule::UserMemoryClean() {
236    if ( fSelectedClusters.size() == 0 ) return;
237   
238    if (fClusters.size() != 0) {
239        vector<EusoCluster*> *pclu = (vector<EusoCluster*>*)MyData()->GetObj("Clusters");
240        EusoCluster *dummy;
241        for(size_t i(0); i<pclu->size(); i++) {
242            dummy = (*pclu)[i];
243            delete dummy;
244        }
245        (*pclu).clear();
246    MyData()->RemoveObj("Clusters");
247    }
248}
249
250// get main cluster id
251void ClusterAnalysisModule::SearchMainCluster() {
252    fMainCluster = 0;
253    Int_t pointsmain = fClusters[0]->GetNumPoints();
254    for( Int_t i=0; i<(Int_t)fClusters.size(); i++ ) {
255        if ( fClusters[i]->GetNumPoints() > pointsmain ) {
256            fMainCluster = i;
257            pointsmain = fClusters[i]->GetNumPoints();
258        }
259    }
260    Msg(EsafMsg::Debug) << "Cluster most populous: n." << fMainCluster << " with " << pointsmain << " points" << MsgDispatch;
261}
262
263// angular distance between two points
264Double_t ClusterAnalysisModule::AngularDistance(Double_t th1, Double_t phi1, Double_t th2, Double_t phi2) {
265    return TMath::ACos( TMath::Sin(th1)*TMath::Sin(th2)*TMath::Cos(phi1)*TMath::Cos(phi2) +
266           TMath::Sin(th1)*TMath::Sin(th2)*TMath::Sin(phi1)*TMath::Sin(phi2) +
267           TMath::Cos(th1)*TMath::Cos(th2) );
268}
269
270// init superclustering process
271void ClusterAnalysisModule::InitSuperClustering() {
272    fNumSuperClusters = 0;
273    fNumNodes = 0;
274    fNumPointsCluster = 0;
275    fBookmark = 0;
276    for(Int_t i=0; i<(Int_t)fClusters.size(); i++) {
277        fId.push_back(0);
278        fFlag1.push_back(kFALSE);
279        fFlag2.push_back(kFALSE);
280    }
281}
282
283// search for superclusters
284void ClusterAnalysisModule::SearchSuperClusters() {
285    while( fNumNodes < (Int_t)fClusters.size() ) {
286        NewSuperCluster();
287        FillSuperCluster();
288        WriteSuperCluster();
289    }
290}
291
292// init a new supercluster
293void ClusterAnalysisModule::NewSuperCluster() {
294    fNumPointsCluster = 1;
295    fNumSuperClusters++;
296    fNumNodes++;
297    fFlag2[fBookmark] = kTRUE;
298    fId[0] = fBookmark;
299    fCounter = 0;
300}
301
302// fill a supercluster
303void ClusterAnalysisModule::FillSuperCluster() {
304    while ( fCounter < fNumPointsCluster ) {
305   
306        Int_t k = fId[fCounter];
307        if ( !fFlag1[k] ) {
308            fFlag1[k] = kTRUE;
309            for (Int_t i=0; i<(Int_t)fClusters.size(); i++) {
310                if ( !fFlag2[i] ) {
311                    Double_t dist = AngularDistance(fThetaMean[i], fPhiMean[i], fThetaMean[k], fPhiMean[k]);
312                    dist -= (fDim[i]+fDim[k]);
313                    if ( dist <= 2.*fThreshold ) {
314                        fFlag2[i] = kTRUE;
315                        fNumPointsCluster++;
316                        fId[fNumPointsCluster-1] = i;
317                        fNumNodes++;
318                    }
319                    fBookmark = i;
320                }
321            }
322        }
323        fCounter++;
324    }
325}
326
327// write a supercluster
328void ClusterAnalysisModule::WriteSuperCluster() {
329    // check cluster significance
330    if ( fNumPointsCluster < 2 )
331        return;
332   
333    for( Int_t i=0; i<fNumPointsCluster; i++ ) {
334        fSelectedClusters.push_back(fId[i]);
335    }
336}
337
338// histogram of clustering
339void ClusterAnalysisModule::HistoClusters() {
340
341    gDirectory->cd("/");
342    string pathname = Form("selectcluster%d", fEv->GetHeader().GetNum());
343    TDirectory *path = new TDirectory(pathname.c_str(),pathname.c_str());
344    path->cd();
345
346    Float_t thmin = fEv->GetRecoPixelData(0)->GetTheta();
347    Float_t thmax = fEv->GetRecoPixelData(0)->GetTheta();
348    Float_t phimin = fEv->GetRecoPixelData(0)->GetPhi();
349    Float_t phimax = fEv->GetRecoPixelData(0)->GetPhi();
350    Float_t tmin = fEv->GetRecoPixelData(0)->GetGtu();
351    Float_t tmax = fEv->GetRecoPixelData(0)->GetGtu();
352
353    for( Int_t i=1; i<(Int_t)fEv->GetHeader().GetNumActiveFee(); i++ ) {
354        if ( fEv->GetRecoPixelData(i)->GetTheta() < thmin )
355            thmin = fEv->GetRecoPixelData(i)->GetTheta();
356        if ( fEv->GetRecoPixelData(i)->GetTheta() > thmax )
357            thmax = fEv->GetRecoPixelData(i)->GetTheta();
358        if ( fEv->GetRecoPixelData(i)->GetPhi() < phimin )
359            phimin = fEv->GetRecoPixelData(i)->GetPhi();
360        if ( fEv->GetRecoPixelData(i)->GetPhi() > phimax )
361            phimax = fEv->GetRecoPixelData(i)->GetPhi();
362        if ( fEv->GetRecoPixelData(i)->GetGtu() < tmin )
363            tmin = fEv->GetRecoPixelData(i)->GetGtu();
364        if ( fEv->GetRecoPixelData(i)->GetGtu() > tmax )
365            tmax = fEv->GetRecoPixelData(i)->GetGtu();
366}
367
368    Int_t bin = (Int_t)Conf()->GetNum("ClusterAnalysisModule.HistogramBinning");
369    TH3F *selhisto = new TH3F("SelHisto","Selected Clustered Points",bin,phimin,phimax,bin,thmin,thmax,
370                              bin,tmin,tmax);
371    TH3F *reclhisto = new TH3F("RiClHisto","RiClusteredPoints",bin,phimin,phimax,bin,thmin,thmax,
372                               bin,tmin,tmax);
373    TH3F *histo = new TH3F("Histo","Points",bin,phimin,phimax,bin,thmin,thmax,bin,tmin,tmax);
374    TH2F *selhisto2 = new TH2F("SelHisto2","Selected Clustered Points",bin,phimin,phimax,bin,thmin,thmax);
375    TH2F *reclhisto2 = new TH2F("RiClHisto2","RiClusteredPoints",bin,phimin,phimax,bin,thmin,thmax);
376
377
378    for( Int_t i=0; i<(Int_t)fEv->GetHeader().GetNumActiveFee(); i++ ) {
379        histo->Fill( fEv->GetRecoPixelData(i)->GetPhi(),fEv->GetRecoPixelData(i)->GetTheta(),
380                     fEv->GetRecoPixelData(i)->GetGtu() );
381    }
382    histo->GetXaxis()->SetTitle("#varphi (rad)");
383    histo->GetYaxis()->SetTitle("#theta (rad)");
384    histo->GetZaxis()->SetTitle("GTU");
385
386    for(Int_t i=0; i<(Int_t)fClusters.size(); i++) {
387        EusoCluster *clu = fClusters[i];
388        for( Int_t j=0; j<clu->GetNumPoints(); j++) {
389            RecoPixelData *pix = fEv->GetRecoPixelData(clu->GetPixelId(j));
390            reclhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
391            reclhisto2->Fill( pix->GetPhi(), pix->GetTheta());
392        }
393    }
394    reclhisto->GetXaxis()->SetTitle("#varphi (rad)");
395    reclhisto->GetYaxis()->SetTitle("#theta (rad)");
396    reclhisto->GetZaxis()->SetTitle("GTU");
397    reclhisto2->GetXaxis()->SetTitle("#varphi (rad)");
398    reclhisto2->GetYaxis()->SetTitle("#theta (rad)");
399
400    thmin = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetTheta();
401    thmax = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetTheta();
402    phimin = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetPhi();
403    phimax = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetPhi();
404    tmin = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetGtu();
405    tmax = fEv->GetRecoPixelData(fClusters[fSelectedClusters[0]]->GetPixelId(0))->GetGtu();
406   
407    for(Int_t i=0; i<(Int_t)fSelectedClusters.size(); i++) {
408        EusoCluster *clu = fClusters[fSelectedClusters[i]];
409        for( Int_t j=0; j<clu->GetNumPoints(); j++) {
410            RecoPixelData *pix = fEv->GetRecoPixelData(clu->GetPixelId(j));
411            if (pix->GetTheta() < thmin ) thmin = pix->GetTheta();
412            if (pix->GetTheta() > thmax ) thmax = pix->GetTheta();
413            if (pix->GetPhi() < phimin ) phimin = pix->GetPhi();
414            if (pix->GetPhi() > phimax ) phimax = pix->GetPhi();
415            if (pix->GetGtu() < tmin ) tmin = pix->GetGtu();
416            if (pix->GetGtu() > tmax ) tmax = pix->GetGtu();
417            selhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
418            selhisto2->Fill( pix->GetPhi(), pix->GetTheta());
419        }
420    }
421    selhisto->GetXaxis()->SetTitle("#varphi (rad)");
422    selhisto->GetYaxis()->SetTitle("#theta (rad)");
423    selhisto->GetZaxis()->SetTitle("GTU");
424    selhisto->GetXaxis()->SetTitle("#varphi (rad)");
425    selhisto->GetYaxis()->SetTitle("#theta (rad)");
426
427
428    TH3F *zoomselhisto = new TH3F("ZSelHisto","Selected Clustered points",bin,phimin,phimax,bin,thmin,thmax,
429                                  bin,tmin,tmax);
430    TH2F *zoomselhisto2 = new TH2F("ZSelHisto2","Selected Clustered points",bin,phimin,phimax,bin,thmin,thmax);
431    for(Int_t i=0; i<(Int_t)fSelectedClusters.size(); i++) {
432        EusoCluster *clu = fClusters[fSelectedClusters[i]];
433        for( Int_t j=0; j<clu->GetNumPoints(); j++) {
434            RecoPixelData *pix = fEv->GetRecoPixelData(clu->GetPixelId(j));
435            zoomselhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
436            zoomselhisto2->Fill( pix->GetPhi(), pix->GetTheta() );
437        }
438    }
439    zoomselhisto->GetXaxis()->SetTitle("#varphi (rad)");
440    zoomselhisto->GetYaxis()->SetTitle("#theta (rad)");
441    zoomselhisto->GetZaxis()->SetTitle("GTU");
442    zoomselhisto2->GetXaxis()->SetTitle("#varphi (rad)");
443    zoomselhisto2->GetYaxis()->SetTitle("#theta (rad)");
444   
445    histo->Write();
446    selhisto->Write();
447    reclhisto->Write();
448    zoomselhisto->Write();
449    delete histo;
450    delete selhisto;
451    delete zoomselhisto;
452    delete reclhisto;
453    selhisto2->Write();
454    reclhisto2->Write();
455    zoomselhisto2->Write();
456    delete selhisto2;
457    delete zoomselhisto2;
458    delete reclhisto2;
459}
Note: See TracBrowser for help on using the repository browser.