source: JEM-EUSO/esaf_cc_at_lal/packages/reconstruction/modules/base/clustering/src/BaseClusteringModule.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: 32.2 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: BaseClusteringModule.cc 2949 2011-06-26 18:39:34Z mabl $
3// R. Pesce created Jan, 19 2004
4
5#include "BaseClusteringModule.hh"
6#include "RecoEvent.hh"
7#include "RecoGlobalData.hh"
8#include "RecoRootEvent.hh"
9#include "RecoCellInfo.hh"
10#include "TDirectory.h"
11#include "TMath.h"
12#include "Config.hh"
13#include "RecoPixelData.hh"
14#include "EusoCluster.hh"
15#include "ERunParameters.hh"
16#include "MedianFit.hh"
17#include "LeastSquaresFit.hh"
18#include "HoughFit.hh"
19
20#include <TH3F.h>
21#include <TH2F.h>
22#include <TCanvas.h>
23#include <TGraph.h>
24
25ClassImp(BaseClusteringModule)
26using namespace TMath;
27using namespace sou;
28
29//____________________________________________________________________________
30BaseClusteringModule::BaseClusteringModule() : 
31    RecoModule(string("BaseClustering")) {
32    //
33    // ctor
34    //
35}
36
37// dtor
38//____________________________________________________________________________
39BaseClusteringModule::~BaseClusteringModule() {
40}
41
42//____________________________________________________________________________
43Bool_t BaseClusteringModule::Init() {
44    //
45    // init
46    //
47    fEv = (RecoEvent*)0;
48    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
49    fAlpha = 0.55;
50    fNumHitsMinimum = Conf()->GetNum("BaseClusteringModule.fNumHitsMinimum");
51    fSignificanceLevel = (Int_t)Conf()->GetNum("BaseClusteringModule.fSignificanceLevel");
52    if ( fSignificanceLevel <= 0 )
53        Msg(EsafMsg::Panic) << "You must specify a positive number in BaseClusteringModule.fSignificanceLevel" << MsgDispatch; 
54    fSignalToNoise    = (Double_t)Conf()->GetNum("BaseClusteringModule.fSignalToNoise");   
55    fFOVThreshold = (Double_t)Conf()->GetNum("BaseClusteringModule.fFOVThreshold");
56    fNSigmaToSelect = (Double_t)Conf()->GetNum("BaseClusteringModule.fNSigmaToSelect");
57    fGTUDiffThreshold = (Int_t)Conf()->GetNum("BaseClusteringModule.fGTUDiffThreshold");
58    return kTRUE;
59}
60
61//____________________________________________________________________________
62Bool_t BaseClusteringModule::PreProcess() {
63    //
64    // pre-process
65    //
66    fTotalNumPoints = 0;
67    fNumPoints = 0;
68    fNumHits = 0;
69    fNumClusters = 0;
70    fNumSigClusters = 0;
71    fNumNodes = 0;
72    fNumPointsCluster = 0;
73    fBookmark = 0;
74    fHitted.clear();
75    fId.clear();
76    fFlag1.clear();
77    fFlag2.clear();
78    fClusters.clear();
79    fSignalFraction = 0.;
80    fContamination = 1.;
81    fGtuMin = 100000;
82    fGtuMax = -100000;
83    fPixelGtuMap.clear();
84       
85    return kTRUE;
86}
87
88//____________________________________________________________________________
89Bool_t BaseClusteringModule::Process( RecoEvent* ev ) {
90    //
91    // process
92    //
93    fEv = ev;
94    fTotalNumPoints = fEv->GetHeader().GetNumActiveFee();
95    fGtuLength = fEv->GetHeader().GetGtuLength();
96   
97
98    if ( fTotalNumPoints <= 0 ) {       
99        Msg(EsafMsg::Info) << "No points in event! Exit from this event." << MsgDispatch;
100        return kFALSE;
101    }
102   
103    // get id of most populous macrocell (main macrocell)
104    fMainMacroCellId = fEv->GetMainMacroCell(1);
105    if ( fEv->GetRecoCellInfo(fMainMacroCellId) == NULL )
106        Msg(EsafMsg::Panic) << "Main macrocell does not correspond to any macrocell simulated in the event." << MsgDispatch;
107
108    // process only signal and filter scattered cherenkov
109    Bool_t ProcessOnlySignal = Config::Get()->GetCF("Reco","RootInputModule")->GetBool("RootInputModule.ProcessOnlySignal");
110    Bool_t FilterScatteredCherenkov = Config::Get()->GetCF("Reco","RootInputModule")->GetBool("RootInputModule.FilterScatteredCherenkov");
111
112    // check if any previous pattern recognition
113    Bool_t isAlreadyPattReco = kFALSE;
114    vector<Int_t> selpix; selpix.clear();
115    RecoGlobalData* gdPattReco = fEv->FindGlobalData("PatternRecognition");
116    if ( gdPattReco ) {
117        isAlreadyPattReco = kTRUE;
118        selpix = *(vector<Int_t>*) gdPattReco->GetObj("SelectedPixels");
119        fTotalNumPoints = selpix.size();
120    }
121       
122
123   
124    // search pixels with num hits > minimum
125    // and calculate the total number of signal hits
126    // in all the event and in the pixels over threshold
127    Int_t numsel(0);
128    fMeanThreshold = 0;
129    Int_t totalSignals(0);
130    Int_t ovThSignals(0);
131    for( Int_t i=0; i<fTotalNumPoints; i++ ) {
132        RecoPixelData *pix;
133        if (!isAlreadyPattReco) pix = (RecoPixelData*)fEv->GetRecoPixelData(i);
134        else pix = (RecoPixelData*)fEv->GetRecoPixelData(selpix[i]); 
135        if (FilterScatteredCherenkov) {
136            pix->SetCountsSignal(pix->GetCountsSignal() - pix->GetCountsCherScatt());
137            pix->SetCounts(pix->GetCounts() - pix->GetCountsCherScatt());
138        }
139       
140        totalSignals += pix->GetCountsSignal();
141       
142        Bool_t thereIsSignal = kTRUE;
143        if (ProcessOnlySignal)  {
144             Int_t signal_counts = pix->GetCountsSignal();
145             Int_t noise_counts  = pix->GetCountsNoise();
146             Double_t ratio (-1);
147             if      (noise_counts > 0   )  ratio = 1.e0*signal_counts/Sqrt(noise_counts);
148             else if ( signal_counts > 0 )  ratio = 1.e10;
149             if (ratio <= fSignalToNoise) thereIsSignal = kFALSE;
150        }
151        if (!thereIsSignal) continue;
152
153        Int_t threshold(0);
154        if ( fNumHitsMinimum >=0 ) threshold = (Int_t)fNumHitsMinimum;
155        else {
156            Int_t uid = pix->GetPixelId();
157            Double_t rate = fEv->GetHeader().GetRunPars()->GetNightGlowRateByUId(uid)*fGtuLength/microsecond;
158            threshold = (Int_t)Ceil( rate + Abs(fNumHitsMinimum)*Sqrt(rate));
159        }
160        numsel++; fMeanThreshold += threshold;
161        if ( pix->GetCounts() >= threshold ) {
162            fHitted.push_back(i);
163            fNumPoints++;
164            fNumHits += pix->GetCounts();
165            ovThSignals += pix->GetCountsSignal();
166        }
167    }
168    if (numsel) fMeanThreshold /= numsel;
169    else fMeanThreshold = 100000;
170   
171    Msg(EsafMsg::Debug) << "Num. of pixels with at least " << fNumHitsMinimum << " hits = " << fNumPoints << MsgDispatch;
172    if (fNumPoints<10) {
173        Msg(EsafMsg::Warning) << "Number of points "<< fNumPoints << " is less then 10 "<< " Exit from this event." << MsgDispatch;
174        fHitted.clear();
175        fMainMacroCellId = 0;
176        fNumPoints = 0;
177        return kFALSE;
178    }
179
180    // fill id and flag vectors
181    for( Int_t i=0; i<fNumPoints; i++ ) {
182        fId.push_back(0);
183        fFlag1.push_back(kFALSE);
184        fFlag2.push_back(kFALSE);
185    }
186   
187    Msg(EsafMsg::Debug) << "Mean Threshold: " << fMeanThreshold << MsgDispatch;
188   
189    //make data maps
190    MakeMaps();
191   
192    // clustering method
193    if (Conf()->GetStr("BaseClusteringModule.fThresholdType")=="natural") {
194        CalculateDensity();
195        CalculateNaturalThreshold();
196            Msg(EsafMsg::Debug) <<"Density: " << fDensity << "  Natural Distance Threshold: " << fThreshold << MsgDispatch;
197        CalculateNumClustersUniform();
198        CalculateNumPointsMinimum();
199        Msg(EsafMsg::Debug) << "Number of clusters in uniform distribution: " << fNumClustersUniform << MsgDispatch;
200        Msg(EsafMsg::Debug) << "Number of points minimum in a significant cluster: " << fNumPointsMinimum << MsgDispatch;
201    } else if (Conf()->GetStr("BaseClusteringModule.fThresholdType")=="fixed") {
202        fThreshold = Conf()->GetNum("BaseClusteringModule.fThreshold");
203            fNumPointsMinimum = (Int_t)Conf()->GetNum("BaseClusteringModule.fNumPointsMinimum");
204        Msg(EsafMsg::Debug) << "Distance Threshold: " << fThreshold << MsgDispatch;
205    } else if (Conf()->GetStr("BaseClusteringModule.fThresholdType")=="fov") {
206        fThreshold = fEv->GetRecoPixelData(fHitted[0])->GetThetaSigma() * fEv->GetRecoPixelData(fHitted[0])->GetPhiSigma();
207        fThreshold = Abs(fFOVThreshold) * Sqrt(fThreshold);
208        Msg(EsafMsg::Debug) << "Distance Threshold: " << fThreshold << MsgDispatch;
209        fNumPointsMinimum = (Int_t)Conf()->GetNum("BaseClusteringModule.fNumPointsMinimum"); 
210    } else {
211        Msg(EsafMsg::Panic) << "Wrong value in BaseClusteringModule.fThresholdType" << MsgDispatch;
212    }
213   
214    SearchClusters(); 
215   
216    //FIXME: just a try
217    /*while ( (fNumSigClusters == 0) && (fSignificanceLevel>1) ) {
218        fSignificanceLevel /= 2;
219        CalculateNumPointsMinimum();
220        SearchClusters();
221        Msg(EsafMsg::Debug) << "Significance level: " << fSignificanceLevel << MsgDispatch;
222        Msg(EsafMsg::Debug) << "Number of significant clusters: " << fNumSigClusters << " in " << fNumClusters << " clusters" << MsgDispatch;
223    }*/
224
225    //find the main cluster (cluster with maximum number of points)
226    Int_t mainclu_id = -1;
227    Int_t maxnumpoints = -1;
228    for(Int_t i(0); i<fNumSigClusters; i++) {
229        if ( fClusters[i]->GetNumPoints() > maxnumpoints ) {
230            maxnumpoints = fClusters[i]->GetNumPoints();
231            mainclu_id = i;
232        }
233    }
234
235    // calculate the background contamination in the main cluster
236    // fit the cluster and do the selection of points
237    if ( mainclu_id >=0 ) {
238        EusoCluster* mainclu = fClusters[mainclu_id];
239        FitCluster(mainclu);
240   
241        fSlopeMainTheta = mainclu->GetSlopeTheta();
242        fOffsetMainTheta = mainclu->GetOffsetTheta();
243        fAbsDevMainTheta = mainclu->GetAbsoluteDeviationTheta();
244        fSlopeMainPhi = mainclu->GetSlopePhi();
245        fOffsetMainPhi = mainclu->GetOffsetPhi();
246        fAbsDevMainPhi = mainclu->GetAbsoluteDeviationPhi();
247
248        // contamination in the cluster
249        Int_t nNoisePixCluster(0);
250        for( Int_t i(0); i<mainclu->GetNumPoints(); i++) {
251            RecoPixelData *pix = (RecoPixelData*)fEv->GetRecoPixelData(mainclu->GetPixelId(i));
252            if (!pix->GetCountsSignal()) nNoisePixCluster++;
253        }
254        fContaminationInCluster = (Double_t)nNoisePixCluster/mainclu->GetNumPoints();
255       
256        // Selection of points
257        for( Int_t i(0); i<fNumPoints; i++ ) {
258            Double_t theta = fEv->GetRecoPixelData(fHitted[i])->GetTheta(); 
259            Double_t phi   = fEv->GetRecoPixelData(fHitted[i])->GetPhi();
260            Double_t time   = fEv->GetRecoPixelData(fHitted[i])->GetGtu()*fGtuLength;
261            if (Conf()->GetStr("BaseClusteringModule.fFitMethod")=="spacetime") {
262                if ( ( Abs(theta -fSlopeMainTheta*time -fOffsetMainTheta) < Abs(fNSigmaToSelect)*fAbsDevMainTheta ) &&
263                 ( Abs(phi -fSlopeMainPhi*time -fOffsetMainPhi) < Abs(fNSigmaToSelect)*fAbsDevMainPhi ) ) 
264                    fSelectedPixels.push_back(fHitted[i]);
265            } else if (Conf()->GetStr("BaseClusteringModule.fFitMethod")=="space") {
266                if  ( Abs(theta -fSlopeMainTheta*phi -fOffsetMainTheta) < Abs(fNSigmaToSelect)*fAbsDevMainTheta ) 
267                    fSelectedPixels.push_back(fHitted[i]);
268            }
269
270        }
271    }
272        if ( fNumSigClusters != 0 ) {
273            for( Int_t i=0; i<fNumSigClusters; i++ ) {
274                MsgForm(EsafMsg::Debug, "Cluster n.%d with %d points",i,fClusters[i]->GetNumPoints());
275        }
276
277        // calculate the quality parameters
278        Int_t nPixNoise(0); // number of pure noise pixels in selection
279        Int_t selectedSignals(0); // number of signal counts in selection
280        for( UInt_t i(0); i<fSelectedPixels.size(); i++ ) {
281            RecoPixelData *pix = (RecoPixelData*)fEv->GetRecoPixelData(fSelectedPixels[i]);
282            if ( pix->GetCountsSignal() ) selectedSignals += pix->GetCountsSignal();
283            else nPixNoise++;
284        }
285        fSignalFraction = (Double_t)selectedSignals/totalSignals;
286        fSigFracOverThresh = (Double_t)selectedSignals/ovThSignals;
287        fContamination = (Double_t)nPixNoise/fSelectedPixels.size();
288        Msg(EsafMsg::Info) << "SignalFraction = " << fSignalFraction << " Contamination = " << fContamination << MsgDispatch;
289        Msg(EsafMsg::Info) << "SignalFraction over threshold = " << fSigFracOverThresh << MsgDispatch;
290        Msg(EsafMsg::Info) << "Contamination in main cluster = " << fContaminationInCluster << MsgDispatch;
291   
292        if ( Conf()->GetBool("BaseClusteringModule.fDoHistogram") ) {
293            HistoCluster();
294        }
295       
296        // add data to event
297        vector<EusoCluster*> *pclu = &fClusters;
298        MyData()->Add("Clusters", pclu);
299        vector<Int_t> *selpix = &fSelectedPixels;
300        MyData()->Add("SelectedPixels", selpix);
301        MyData()->Add("NumHitsMinimum",fNumHitsMinimum);
302        MyData()->Add("MeanThreshold",fMeanThreshold);
303        MyData()->Add("SignificanceLevel",fSignificanceLevel);
304        MyData()->Add("NumPointsMinimum",fNumPointsMinimum);
305        MyData()->Add("DistanceThreshold",fThreshold);
306        fSignificanceLevel = (Int_t)Conf()->GetNum("BaseClusteringModule.fSignificanceLevel"); //reset
307
308        return kTRUE;
309    }
310       
311    if ( Conf()->GetBool("BaseClusteringModule.fDoHistogram") ) {
312        HistoCluster();
313    }   
314
315    fSignificanceLevel = (Int_t)Conf()->GetNum("BaseClusteringModule.fSignificanceLevel"); //reset
316    return kFALSE;
317}
318
319//___________________________________________________________________________
320void BaseClusteringModule::MakeMaps() {
321    //
322    // Organize the event data to optimize the clustering process
323    //
324
325    // sort fHitted by gtu
326    /*Int_t n = fHitted.size();
327    for (Int_t j=n/2; 0<j; j/=2) {
328        for(Int_t i(0); i<n; i++) {
329            for(Int_t k=i-j; 0<=k; k-=j) {
330                if (fEv->GetRecoPixelData(fHitted[k+j])->GetGtu() < fEv->GetRecoPixelData(fHitted[k])->GetGtu() ) {
331                    Int_t dummy = fHitted[k];
332                    fHitted[k] = fHitted[k+j];
333                    fHitted[k+j] = dummy;
334                }
335            }
336        }
337    }*/
338   
339    // fill the gtu-pixel map
340    // for each gtu store the corresponding fHitted entries
341    // (fHitted is the vector of the ids of recopixeldata over threshold)
342    for(Int_t i(0); i<fNumPoints; i++) {
343        Int_t gtu = fEv->GetRecoPixelData(fHitted[i])->GetGtu();
344        if ( gtu < fGtuMin ) fGtuMin = gtu;
345        if ( gtu > fGtuMax ) fGtuMax = gtu;
346        fPixelGtuMap[gtu].push_back(i);
347    }
348   
349   
350}
351
352//___________________________________________________________________________
353Bool_t BaseClusteringModule::PostProcess() {
354    //
355    // post-process
356    //
357
358    if ( fTotalNumPoints == 0 ) return kTRUE;
359    if ( fNumPoints < 10 ) {
360        fEv = (RecoEvent*)0; 
361        return kTRUE;
362    }
363   
364    // Update data in PatternRecognition Global data
365    RecoGlobalData* data = fEv->FindCreateGlobalData("PatternRecognition");
366    data->CleanMaps();
367    data->SetIssuer(GetName());
368    data->Add("SelectedPixels",&fSelectedPixels);
369    data->Add("NumPointsMinimum",fNumPointsMinimum);
370    data->Add("MeanThreshold",fMeanThreshold);
371   
372       
373    fEv = (RecoEvent*)0;
374    return kTRUE;
375     
376}
377
378//___________________________________________________________________________
379Bool_t BaseClusteringModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
380    fRecoRootEvent->GetRecoClustering().SetNumPoints((Int_t)fSelectedPixels.size());
381    fRecoRootEvent->GetRecoClustering().SetSignalFraction(fSignalFraction);
382    fRecoRootEvent->GetRecoClustering().SetSigFracOverThresh(fSigFracOverThresh);
383    fRecoRootEvent->GetRecoClustering().SetContamination(fContamination);
384    fRecoRootEvent->GetRecoClustering().SetContaminationInCluster(fContaminationInCluster);
385    fRecoRootEvent->GetRecoClustering().SetAbsDevMainTheta(fAbsDevMainTheta);
386    fRecoRootEvent->GetRecoClustering().SetAbsDevMainPhi(fAbsDevMainPhi);
387   
388    return kTRUE;
389}
390
391//___________________________________________________________________________
392Bool_t BaseClusteringModule::Done() {
393    //
394    // done
395    //
396    fHitted.clear();
397    fId.clear();
398    fFlag1.clear();
399    fFlag2.clear();
400    fClusters.clear();
401    fPixelGtuMap.clear();
402    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
403    return kTRUE;
404}
405
406
407//___________________________________________________________________________
408void BaseClusteringModule::UserMemoryClean() {
409    //
410    // memory clean
411    //
412    if ( fTotalNumPoints == 0 || fNumPoints < 10 || fNumSigClusters == 0 ) return;
413
414    vector<EusoCluster*> *pclu = (vector<EusoCluster*>*)MyData()->GetObj("Clusters");
415    EusoCluster *dummy(NULL);
416    for(size_t i(0); i<pclu->size(); i++) {
417        dummy = (*pclu)[i];
418        delete dummy;
419    }
420    (*pclu).clear();
421    MyData()->RemoveObj("Clusters");
422
423    vector<Int_t> *selpix = (vector<Int_t>*)MyData()->GetObj("SelectedPixels");
424    if (selpix) {
425        selpix->clear();
426        MyData()->RemoveObj("SelectedPixels");
427        vector<Int_t> dummyV;
428        dummyV.swap(*selpix);
429    }
430
431}
432
433//___________________________________________________________________________
434void BaseClusteringModule::CalculateDensity() {
435    //
436    // calculate density of points
437    //
438    if ( Conf()->GetStr("BaseClusteringModule.fDensityMethod")=="mainmc" ) {
439        DensityMainMC();
440    } else if ( Conf()->GetStr("BaseClusteringModule.fDensityMethod")=="event" ) {
441        DensityEvent();
442    } else {
443        Msg(EsafMsg::Panic) <<"Wrong config value BaseClusteringModule.fDensityMethod" << MsgDispatch; 
444    }
445}
446
447//___________________________________________________________________________
448void BaseClusteringModule::DensityMainMC() {
449    //
450    // Calculate density from main macrocell
451    //
452    Double_t nph = (Double_t)fEv->GetNumPx(fMainMacroCellId, fMeanThreshold);
453    Double_t area = (Double_t)fEv->GetThetaRange(fMainMacroCellId, fMeanThreshold);
454    area *= (Double_t)fEv->GetPhiRange(fMainMacroCellId, fMeanThreshold);
455    fDensity = nph / area;
456    if (fDensity==0) DensityEvent();
457}
458
459//___________________________________________________________________________
460void BaseClusteringModule::DensityEvent() {
461    //
462    // Calculate density from complete event
463    //
464    Float_t thmin = fEv->GetRecoPixelData(fHitted[0])->GetTheta();
465    Float_t thmax = fEv->GetRecoPixelData(fHitted[0])->GetTheta();
466    Float_t phimin = fEv->GetRecoPixelData(fHitted[0])->GetPhi();
467    Float_t phimax = fEv->GetRecoPixelData(fHitted[0])->GetPhi();
468   
469    for(Int_t i=1; i<fNumPoints; i++) {
470        Int_t j = fHitted[i];
471        if ( fEv->GetRecoPixelData(j)->GetTheta() < thmin )
472            thmin = fEv->GetRecoPixelData(j)->GetTheta();
473        if ( fEv->GetRecoPixelData(j)->GetTheta() > thmax )
474            thmax = fEv->GetRecoPixelData(j)->GetTheta();
475        if ( fEv->GetRecoPixelData(j)->GetPhi() < phimin )
476            phimin = fEv->GetRecoPixelData(i)->GetPhi();
477        if ( fEv->GetRecoPixelData(j)->GetPhi() > phimax )
478            phimax = fEv->GetRecoPixelData(i)->GetPhi();
479    }
480    Double_t area = Abs(( fEv->GetEndGtu()-fEv->GetStartGtu() )*(thmax-thmin)*(phimax-phimin));
481    fDensity = (Double_t) fNumPoints / area;
482}
483
484//___________________________________________________________________________
485void BaseClusteringModule::CalculateNaturalThreshold() {
486    //
487    // calculate natural threshold
488    //
489    fThreshold = Sqrt(Log(2)/
490        (fDensity*fAlpha*Pi()));
491}
492
493//___________________________________________________________________________
494Double_t BaseClusteringModule::CalculateDistance(Int_t n1, Int_t n2) {
495    //
496    // calculate distance between two points
497    //
498    RecoPixelData *px1=fEv->GetRecoPixelData(n1);
499    RecoPixelData *px2=fEv->GetRecoPixelData(n2);
500    return ACos( Sin(px1->GetTheta())*Sin(px2->GetTheta())*
501           Cos(px1->GetPhi())*Cos(px2->GetPhi()) +
502           Sin(px1->GetTheta())*Sin(px2->GetTheta())*
503           Sin(px1->GetPhi())*Sin(px2->GetPhi()) +
504           Cos(px1->GetTheta())*Cos(px2->GetTheta()) );
505}
506
507//___________________________________________________________________________
508Int_t BaseClusteringModule::CalculateTimeDifference(Int_t n1, Int_t n2) {
509    //
510    // calculate the time difference between two points in GTU
511    //
512    RecoPixelData *px1=fEv->GetRecoPixelData(n1);
513    RecoPixelData *px2=fEv->GetRecoPixelData(n2);
514    return Abs(px1->GetGtu() - px2->GetGtu());
515}
516
517//___________________________________________________________________________
518void BaseClusteringModule::CalculateNumClustersUniform() {
519    //
520    // calculate number of clusters expected in a distribution uniform
521    //
522    fNumClustersUniform = (Int_t)( 1 + ( fNumPoints - 1 ) *
523        Exp( -fAlpha * fDensity * fThreshold * fThreshold * Pi() ) );
524}
525
526//___________________________________________________________________________
527void BaseClusteringModule::CalculateNumPointsMinimum() {
528    //
529    // calculate the minimum number of points in a significant cluster
530    //
531    fNumPointsMinimum = (Int_t) ( fNumPoints / fNumClustersUniform
532        + fSignificanceLevel * Sqrt( (Double_t)( fNumPoints / fNumClustersUniform ) ) );
533}
534
535//___________________________________________________________________________
536void BaseClusteringModule::SearchClusters() {
537    //
538    // search for clusters
539    //
540    while( fNumNodes < fNumPoints ) {
541        NewCluster();
542        FillCluster();
543        WriteCluster();
544    }
545}
546
547//___________________________________________________________________________
548void BaseClusteringModule::NewCluster() {
549    //
550    // init a new cluster
551    //
552    fNumPointsCluster = 1;
553    fNumClusters++;
554    fNumNodes++;
555    fFlag2[fBookmark] = kTRUE;
556    fId[0] = fBookmark;
557    fCounter = 0;
558}
559
560//___________________________________________________________________________
561void BaseClusteringModule::FillCluster() {
562    //
563    // fill a cluster
564    //
565    /*while ( fCounter < fNumPointsCluster ) {
566   
567        Int_t k = fId[fCounter];
568        if ( !fFlag1[k] ) {
569            fFlag1[k] = kTRUE;
570            for (Int_t i=0; i<fNumPoints; i++) {
571                if ( !fFlag2[i] ) {
572                    if ( (CalculateDistance(fHitted[i], fHitted[k]) <= fThreshold) &&
573                            (CalculateTimeDifference(fHitted[i], fHitted[k]) <= fGTUDiffThreshold) ) {   
574                        fFlag2[i] = kTRUE;
575                        fNumPointsCluster++;
576                        fId[fNumPointsCluster-1] = i;
577                        fNumNodes++;
578                    }
579                    fBookmark = i;
580                }
581            }
582        }
583        fCounter++;
584    }*/
585    while ( fCounter < fNumPointsCluster ) {
586   
587        Int_t k = fId[fCounter];
588        if ( !fFlag1[k] ) {
589            fFlag1[k] = kTRUE;
590            Int_t gtu = fEv->GetRecoPixelData(fHitted[k])->GetGtu();
591            for (Int_t j=-1; j<2; j++) {
592                if ( (gtu+j)<fGtuMin || (gtu+j)>fGtuMax ) continue;
593                for (Int_t i=0; i<(Int_t)(fPixelGtuMap[gtu+j]).size(); i++) {
594                    Int_t q = (fPixelGtuMap[gtu+j])[i];
595                    if ( !fFlag2[q] ) {
596                        if  (CalculateDistance(fHitted[q], fHitted[k]) <= fThreshold) {
597                            fFlag2[q] = kTRUE;
598                            fNumPointsCluster++;
599                            fId[fNumPointsCluster-1] = q;
600                            fNumNodes++;
601                        }
602                        //fBookmark = q;
603                    }
604                }
605            }
606        }
607        fCounter++;
608
609    }
610    fBookmark++;
611}
612
613//___________________________________________________________________________
614void BaseClusteringModule::WriteCluster() {
615    //
616    // close and write the current cluster if is significant
617    //
618
619    // check cluster significance
620    if ( fNumPointsCluster < fNumPointsMinimum )
621        return;
622
623    // write the cluster
624    fNumSigClusters++;
625    EusoCluster *cl = new EusoCluster();
626    Double_t thmin, thmax, phimin, phimax;
627    Int_t gtumin, gtumax;
628    thmin = fEv->GetRecoPixelData(fHitted[fId[0]])->GetTheta(); 
629    thmax = fEv->GetRecoPixelData(fHitted[fId[0]])->GetTheta();
630    phimin = fEv->GetRecoPixelData(fHitted[fId[0]])->GetPhi();
631    phimax = fEv->GetRecoPixelData(fHitted[fId[0]])->GetPhi();
632    gtumin = fEv->GetRecoPixelData(fHitted[fId[0]])->GetGtu();
633    gtumax = fEv->GetRecoPixelData(fHitted[fId[0]])->GetGtu();
634   
635    Double_t theta(0),phi(0);
636    Int_t gtu(0);
637    for( Int_t i=0; i<fNumPointsCluster; i++ ) {
638        theta = fEv->GetRecoPixelData(fHitted[fId[i]])->GetTheta();
639        phi = fEv->GetRecoPixelData(fHitted[fId[i]])->GetPhi(); 
640        gtu = fEv->GetRecoPixelData(fHitted[fId[i]])->GetGtu();
641        if ( theta < thmin ) thmin = theta;
642        if ( theta > thmax ) thmax = theta;
643        if ( phi < phimin ) phimin = phi;
644        if ( phi > phimax ) phimax = phi;
645        if ( gtu < gtumin ) gtumin = gtu;
646        if ( gtu > gtumax ) gtumax = gtu;
647        cl->AddPoint( fHitted[fId[i]] );
648    }
649    cl->SetThetaMin( thmin );
650    cl->SetThetaMax( thmax );
651    cl->SetPhiMin( phimin );
652    cl->SetPhiMax( phimax );
653    cl->SetGtuMin( gtumin );
654    cl->SetGtuMax( gtumax );
655    cl->SetThreshold( (Int_t)fNumHitsMinimum );
656    fClusters.push_back(cl);
657}
658
659
660//___________________________________________________________________________
661EusoCluster* BaseClusteringModule::GetCluster( Int_t i ) {
662    //
663    // get a specified cluster
664    //
665    if ( i >= 0 && i < (Int_t)fClusters.size() )
666        return fClusters[i];
667    else
668        return NULL;
669}
670
671//___________________________________________________________________________
672void BaseClusteringModule::FitCluster( EusoCluster *clu ) {
673    //
674    // fit the points of a cluster; use linear or median or hough fit
675    //
676
677    vector<Double_t> thetavector;
678    vector<Double_t> phivector;
679    vector<Double_t> timevector;
680    for( Int_t i=0; i<clu->GetNumPoints(); i++ ) {
681        RecoPixelData *pix = fEv->GetRecoPixelData( clu->GetPixelId(i) );
682        thetavector.push_back(pix->GetTheta());
683        phivector.push_back(pix->GetPhi());
684        timevector.push_back(pix->GetGtu()*fGtuLength);
685    }
686   
687    if ( Conf()->GetStr("BaseClusteringModule.fFitMethod")=="spacetime" ) {
688        MedianFit *mftheta = new MedianFit( clu->GetNumPoints(), timevector, thetavector );
689        clu->SetFitted( kTRUE );
690        clu->SetSlopeTheta( mftheta->GetSlope() );
691        clu->SetOffsetTheta( mftheta->GetOffset() );
692        clu->SetAbsoluteDeviationTheta( mftheta->GetAbsoluteDeviation() );
693        MedianFit *mfphi = new MedianFit( clu->GetNumPoints(), timevector, phivector );
694        clu->SetSlopePhi( mfphi->GetSlope() );
695        clu->SetOffsetPhi( mfphi->GetOffset() );
696        clu->SetAbsoluteDeviationPhi( mfphi->GetAbsoluteDeviation() );
697        //cout << mftheta->GetSlope() << " " << mftheta->GetOffset() << " " << mftheta->GetAbsoluteDeviation() << endl;
698        //cout << mfphi->GetSlope() << " " << mfphi->GetOffset() << " " << mfphi->GetAbsoluteDeviation() << endl;
699    } else if ( Conf()->GetStr("BaseClusteringModule.fFitMethod")=="space" ) {
700        MedianFit *mf = new MedianFit( clu->GetNumPoints(), phivector, thetavector );
701        clu->SetFitted( kTRUE );
702        clu->SetSlopeTheta( mf->GetSlope() );
703        clu->SetOffsetTheta( mf->GetOffset() );
704        clu->SetAbsoluteDeviationTheta( mf->GetAbsoluteDeviation() );
705    } else {
706        Msg(EsafMsg::Panic) << "Wrong config value BaseClusteringModule.fFitMethod" << MsgDispatch; 
707    }
708}
709
710//___________________________________________________________________________
711void BaseClusteringModule::HistoCluster() {
712    //
713    // make an histogram of the cluster
714    //
715
716    gDirectory->cd("/");
717    string pathname = Form("cluster%d", fEv->GetHeader().GetNum());
718    TDirectory *path = new TDirectory(pathname.c_str(),pathname.c_str());
719    path->cd();
720   
721    Float_t thmin = fEv->GetRecoPixelData(0)->GetTheta();
722    Float_t thmax = fEv->GetRecoPixelData(0)->GetTheta();
723    Float_t phimin = fEv->GetRecoPixelData(0)->GetPhi();
724    Float_t phimax = fEv->GetRecoPixelData(0)->GetPhi();
725    Float_t tmin = fEv->GetRecoPixelData(0)->GetGtu();
726    Float_t tmax = fEv->GetRecoPixelData(0)->GetGtu();
727   
728    for( Int_t i=1; i<fTotalNumPoints; i++ ) {
729        if ( fEv->GetRecoPixelData(i)->GetTheta() < thmin )
730            thmin = fEv->GetRecoPixelData(i)->GetTheta();
731        if ( fEv->GetRecoPixelData(i)->GetTheta() > thmax )
732            thmax = fEv->GetRecoPixelData(i)->GetTheta();
733        if ( fEv->GetRecoPixelData(i)->GetPhi() < phimin )
734            phimin = fEv->GetRecoPixelData(i)->GetPhi();
735        if ( fEv->GetRecoPixelData(i)->GetPhi() > phimax )
736            phimax = fEv->GetRecoPixelData(i)->GetPhi();
737        if ( fEv->GetRecoPixelData(i)->GetGtu() < tmin )
738            tmin = fEv->GetRecoPixelData(i)->GetGtu();
739        if ( fEv->GetRecoPixelData(i)->GetGtu() > tmax )
740            tmax = fEv->GetRecoPixelData(i)->GetGtu();
741    }
742       
743    Int_t bin = (Int_t)Conf()->GetNum("BaseClusteringModule.fHistogramBinning");
744    TH3F *clhisto = new TH3F("ClHisto","Clustered Points",bin,phimin,phimax,bin,thmin,thmax,bin,tmin,tmax);
745    TH3F *hithisto = new TH3F("HitHisto","Points with minimum hits",bin,phimin,phimax,bin,thmin,thmax,
746                              bin,tmin,tmax);
747    TH3F *histo = new TH3F("Histo","Points",bin,phimin,phimax,bin,thmin,thmax,bin,tmin,tmax);
748    TH2F *clhisto2 = new TH2F("ClHisto2","Clustered Points",bin,phimin,phimax,bin,thmin,thmax);
749    TH2F *hithisto2 = new TH2F("HitHisto2","Points with minimum hits",bin,phimin,phimax,bin,thmin,thmax);
750    TH2F *histo2 = new TH2F("Histo2","Points",bin,phimin,phimax,bin,thmin,thmax);
751    TH3F *selhisto = new TH3F("SelHisto","Selected Points",bin,phimin,phimax,bin,thmin,thmax,bin,tmin,tmax);
752    TH2F *selhisto2 = new TH2F("SelHisto2","Selected Points",bin,phimin,phimax,bin,thmin,thmax);
753 
754 
755    for( Int_t i=0; i<fTotalNumPoints; i++ ) {
756        histo->Fill( fEv->GetRecoPixelData(i)->GetPhi(),fEv->GetRecoPixelData(i)->GetTheta(), 
757                     fEv->GetRecoPixelData(i)->GetGtu() );
758        histo2->Fill( fEv->GetRecoPixelData(i)->GetPhi(),fEv->GetRecoPixelData(i)->GetTheta());
759    }
760    histo->GetXaxis()->SetTitle("#varphi (rad)");
761    histo->GetYaxis()->SetTitle("#theta (rad)");
762    histo->GetZaxis()->SetTitle("GTU");
763    histo2->GetXaxis()->SetTitle("#varphi (rad)");
764    histo2->GetYaxis()->SetTitle("#theta (rad)");
765   
766    for( Int_t i=0; i<fNumPoints; i++ ) {
767        hithisto->Fill(fEv->GetRecoPixelData(fHitted[i])->GetPhi(),
768            fEv->GetRecoPixelData(fHitted[i])->GetTheta(), fEv->GetRecoPixelData(fHitted[i])->GetGtu() );
769        hithisto2->Fill(fEv->GetRecoPixelData(fHitted[i])->GetPhi(),
770            fEv->GetRecoPixelData(fHitted[i])->GetTheta());
771    }
772    hithisto->GetXaxis()->SetTitle("#varphi (rad)");
773    hithisto->GetYaxis()->SetTitle("#theta (rad)");
774    hithisto->GetZaxis()->SetTitle("GTU");
775    hithisto2->GetXaxis()->SetTitle("#varphi (rad)");
776    hithisto2->GetYaxis()->SetTitle("#theta (rad)");
777   
778    if (!fClusters.size()) {
779        hithisto->Write();
780        histo->Write();
781        delete hithisto;
782        delete histo;
783        hithisto2->Write();
784        histo2->Write();
785        delete hithisto2;
786        delete histo2;
787        return;
788    }
789   
790    for(Int_t i=0; i<(Int_t)fClusters.size(); i++) {
791        EusoCluster *clu = fClusters[i];
792        for( Int_t j=0; j<clu->GetNumPoints(); j++) {
793            RecoPixelData *pix = fEv->GetRecoPixelData(clu->GetPixelId(j));
794            clhisto->Fill( pix->GetPhi(), pix->GetTheta(), pix->GetGtu() );
795            clhisto2->Fill( pix->GetPhi(), pix->GetTheta());
796        }
797    }
798    clhisto->GetXaxis()->SetTitle("#varphi (rad)");
799    clhisto->GetYaxis()->SetTitle("#theta (rad)");
800    clhisto->GetZaxis()->SetTitle("GTU");
801    clhisto2->GetXaxis()->SetTitle("#varphi (rad)");
802    clhisto2->GetYaxis()->SetTitle("#theta (rad)");
803   
804    Int_t numsel = fSelectedPixels.size();
805    Double_t th[numsel]; Double_t phi[numsel];
806    Double_t thline[numsel];
807    for( Int_t i(0); i<numsel; i++ ) {
808        th[i] = fEv->GetRecoPixelData(fSelectedPixels[i])->GetTheta();
809        phi[i] = fEv->GetRecoPixelData(fSelectedPixels[i])->GetPhi();
810        selhisto2->Fill(phi[i], th[i]);
811        selhisto->Fill(phi[i], th[i], fEv->GetRecoPixelData(fSelectedPixels[i])->GetGtu());
812    }
813    selhisto2->GetXaxis()->SetTitle("#varphi (rad)");
814    selhisto2->GetYaxis()->SetTitle("#theta (rad)");
815    selhisto->GetXaxis()->SetTitle("#varphi (rad)");
816    selhisto->GetYaxis()->SetTitle("#theta (rad)");
817    selhisto->GetZaxis()->SetTitle("GTU");
818   
819    TGraph *selpoints = new TGraph(numsel,phi,th);
820    TGraph *mainline = new TGraph(numsel,phi,thline);
821    TCanvas *cv = new TCanvas(); cv->Divide(1,1); cv->cd(1);
822    selpoints->Draw("AP");
823    mainline->Draw("PLsame");
824    cv->Write(); delete cv;
825    delete selpoints;
826
827    TCanvas *cv1 = new TCanvas(); cv1->Divide(1,2); cv1->cd(1);
828    selhisto2->Draw();
829    cv1->cd(2); selhisto->Draw();
830    cv1->Write(); delete cv1;
831
832    clhisto->Write();
833    hithisto->Write();
834    histo->Write();
835    delete clhisto;
836    delete hithisto;
837    delete histo;
838    clhisto2->Write();
839    hithisto2->Write();
840    histo2->Write();
841    delete clhisto2;
842    delete hithisto2;
843    delete histo2;
844}
Note: See TracBrowser for help on using the repository browser.