source: JEM-EUSO/esaf_cc_at_lal/packages/reconstruction/modules/base/clustering/src/HoughModule.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: 40.2 KB
Line 
1
2// $Id: HoughModule.cc
3// Author: Elena Taddei   
4
5/*****************************************************************************
6 * ESAF: Euso Simulation and Analysis Framework                              *
7 *                                                                           *
8 *  Id: HoughModule                                                          *
9 *                                                                           *
10 *****************************************************************************/
11
12//_____________________________________________________________________________
13//
14// HoughModule
15//
16// <extensive class description>
17//
18//   Config file parameters
19//   ======================
20//
21//   <parameter name>: <parameter description>
22//   -Valid options: <available options>
23//
24
25#include "HoughModule.hh"
26#include "RecoEvent.hh"
27#include "RecoRootEvent.hh"
28#include "RecoGlobalData.hh"
29#include "RecoCellInfo.hh"
30#include "TDirectory.h"
31#include "TMath.h"
32#include "Config.hh"
33#include "RecoPixelData.hh"
34#include "EusoCluster.hh"
35#include "MedianFit.hh"
36#include "HoughFit.hh"
37#include "ERunParameters.hh"
38#include "EDetector.hh"
39#include "EConst.hh"
40
41#include <TH3F.h>
42#include <TH2F.h>
43#include <TH1F.h>
44#include <TGraph.h>
45#include <TCanvas.h>
46
47using namespace TMath;
48using namespace EConst;
49using namespace sou;
50
51ClassImp(HoughModule)
52
53// FIXME memory leaks to be fixed asap. D.Naumov
54
55//_____________________________________________________________________________
56HoughModule::HoughModule():RecoModule("HoughTransform") {
57    //
58    // ctor
59    //
60}
61
62//_____________________________________________________________________________
63HoughModule::~HoughModule() {
64    //
65    // dtor
66    //
67}
68
69//_____________________________________________________________________________
70Bool_t HoughModule::Init() {
71    //
72    // init
73    //
74   
75   
76    fEv = (RecoEvent*)NULL;
77   
78    fNumPointsMinimum = (Int_t)Conf()->GetNum("HoughModule.fNumPointsMinimum");
79    fSignalToNoise    = (Double_t)Conf()->GetNum("HoughModule.fSignalToNoise");
80    fDebugPlots       = Conf()->GetBool("HoughModule.fDebugPlots");
81    fProcedure        = Conf()->GetStr("HoughModule.fProcedure");
82    Msg(EsafMsg::Info) << "Initializing... Procedure: \""<< fProcedure.c_str() << "\"" << MsgDispatch;
83   
84    fListObj = new TList;
85
86    return kTRUE;
87}
88
89//_____________________________________________________________________________
90Bool_t HoughModule::PreProcess() {
91    //
92    // pre-process
93    //
94       
95    fTotalNumPoints = 0;
96    fNumPointsAll = 0;
97    fClustersWritten = 0;
98   
99    fPreviousPattReco.clear();
100    fHittedFinal.clear();
101    fHittedFinalMacrocell.clear();
102    fHittedSelHough.clear();
103    fHittedSel2Hough.clear();
104    fHittedAll.clear();
105    fIdMacrocellsSelected.clear();
106    fHittedHoughMerged.clear();
107       
108    fClusters.clear();
109       
110    return kTRUE;
111}
112
113// process
114//_____________________________________________________________________________
115Bool_t HoughModule::Process( RecoEvent* ev ) {
116       
117    fEv = ev;
118    fTotalNumPoints = fEv->GetHeader().GetNumActiveFee();
119    fIdEv = fEv->GetHeader().GetNum();
120    fGtuLength = fEv->GetHeader().GetGtuLength();
121   
122    // check if any previous pattern recognition
123    fIsAlreadyPattReco = kFALSE;
124    RecoGlobalData* gdPattReco = fEv->FindGlobalData("PatternRecognition");
125    if ( gdPattReco ) {
126        fIsAlreadyPattReco = kTRUE;
127        fPreviousPattReco = *(vector<Int_t>*) gdPattReco->GetObj("SelectedPixels");
128        fTotalNumPoints = fPreviousPattReco.size();
129    }   
130   
131    if ( fTotalNumPoints <= fNumPointsMinimum ) {       
132        Msg(EsafMsg::Info) << "Not enough points in event! Exit from this event." << MsgDispatch;
133        return kFALSE;
134    }else{
135        Msg(EsafMsg::Debug) << "fTotalNumPoints = " << fTotalNumPoints << MsgDispatch;
136    }
137   
138    fMainMacroCellId = fEv->GetMainMacroCell(1);
139    if ( fEv->GetRecoCellInfo(fMainMacroCellId) == NULL )
140        Msg(EsafMsg::Info) << "Main macrocell does not correspond to any macrocell hit." << MsgDispatch;
141   
142     
143   
144    Bool_t ProcessOnlySignal = Config::Get()->GetCF("Reco","RootInputModule")->GetBool("RootInputModule.ProcessOnlySignal");
145    Bool_t FilterScatteredCherenkov = Config::Get()->GetCF("Reco","RootInputModule")->GetBool("RootInputModule.FilterScatteredCherenkov");
146
147    if (ProcessOnlySignal)  {
148        for( Int_t i=0; i<fTotalNumPoints; i++ ) {
149            RecoPixelData *pix;
150            if (!fIsAlreadyPattReco) pix = (RecoPixelData*)fEv->GetRecoPixelData(i);
151            else pix = (RecoPixelData*)fEv->GetRecoPixelData(fPreviousPattReco[i]);
152            if (FilterScatteredCherenkov) {
153                 pix->SetCountsSignal(pix->GetCountsSignal() - pix->GetCountsCherScatt());
154                 pix->SetCounts(pix->GetCounts() - pix->GetCountsCherScatt());
155             }
156             Int_t signal_counts = pix->GetCountsSignal();
157             Int_t noise_counts  = pix->GetCountsNoise();
158             Double_t ratio (-1);
159             if      (noise_counts > 0   )  ratio = 1.e0*signal_counts/Sqrt(noise_counts);
160             else if ( signal_counts > 0 )  ratio = 1.e10;
161             if (ratio >= fSignalToNoise) fHittedFinal.push_back(i);
162        }
163        fFinalNumHitsMinimum = 1;
164        WriteSingleCluster();
165        return kTRUE;
166    }
167   
168    // create a new directory only if debug plots are enabled
169    if (fDebugPlots) {
170        gDirectory->cd("/");
171        string pathnameevent = Form("HoughModule%d", fIdEv);
172        TDirectory *path =  new TDirectory(pathnameevent.c_str(),pathnameevent.c_str());
173        path->cd();
174    }
175
176    Bool_t ret = kFALSE;
177    // check if clustering exist
178    if (fIsAlreadyPattReco) ret = ProcessHoughAll();
179    else if( fProcedure == "single" )       ret = ProcessHoughSingleMacrocell();
180    else if( fProcedure == "velocity" ) ret = ProcessVelocity();
181    else if( fProcedure == "all" )          ret = ProcessHoughAll();
182    else Msg(EsafMsg::Panic) << "Wrong procedure in HoughModule.fProcedure. Can be single, velocity or all." << MsgDispatch;
183   
184    if(!ret){
185        Msg(EsafMsg::Info) << "Problems ret=0.. " << MsgDispatch;
186        return kFALSE;
187    } else {
188        if(fHittedFinal.size()<10 || fHittedFinal.size()>100000)        return kFALSE;
189            WriteSingleCluster();
190        return kTRUE;
191    } 
192    return kFALSE; 
193   
194}
195
196//_____________________________________________________________________________
197Bool_t HoughModule::ProcessHoughSingleMacrocell(){
198    //
199    // Do Hough Transform for each macrocell hitted and trying to select
200    // which are hitted by real signal
201    //
202       
203    Msg(EsafMsg::Info) << " ...ProcessHoughSingleMacrocell()" << MsgDispatch;
204   
205    Int_t numcountsmin = (Int_t)Conf()->GetNum("HoughModule.fMinCountsHoughSingleMacrocell");
206       
207    RecoEventHeader header = fEv->GetHeader();
208    ERunParameters *runpars = header.GetRunPars();
209       
210    map< Int_t,vector<Int_t> > mapIdfHittedSingleCell;
211    for(Int_t i=0;i<200;i++){
212        RecoCellInfo *info = fEv->GetRecoCellInfo(i);
213        if( info != NULL){
214                vector<Int_t> empty;
215                pair< Int_t,vector<Int_t> > p(i,empty);
216                mapIdfHittedSingleCell.insert(p);
217        }
218    }
219    Msg(EsafMsg::Info) << "Number of macrocells hitted = " << mapIdfHittedSingleCell.size() << MsgDispatch;
220
221    map< Int_t,vector<Int_t> >::iterator it;
222    for( Int_t i=0; i<fTotalNumPoints; i++ ) {
223        Int_t threshold(0);
224        if ( numcountsmin >=0 ) threshold = numcountsmin;
225        else {
226            Int_t uid = fEv->GetRecoPixelData(i)->GetPixelId();
227            Double_t rate = fEv->GetHeader().GetRunPars()->GetNightGlowRateByUId(uid)*fGtuLength/microsecond; 
228            threshold = (Int_t) Ceil(rate + Abs(numcountsmin)*Sqrt(rate));
229        }
230            if ( fEv->GetRecoPixelData(i)->GetCounts() >= threshold ){
231                Int_t cellid=fEv->GetRecoPixelData(i)->GetMacroCellId();
232                it = mapIdfHittedSingleCell.find(cellid);
233                it->second.push_back(i);
234            }
235    }
236   
237    if (numcountsmin >= 0) {
238        Msg(EsafMsg::Info) << "Selected pixels with counts > " << numcountsmin << MsgDispatch;
239    } else {
240        Msg(EsafMsg::Info) << "Selected pixels with counts over " << -numcountsmin << " sigma from background mean" << MsgDispatch;
241    }
242    for( it = mapIdfHittedSingleCell.begin(); it != mapIdfHittedSingleCell.end(); it++ ) {
243        Msg(EsafMsg::Info)<<"\n\t\t---->macrocellid = "<<it->first<<"\t with "<<(it->second).size()<<" points "<<MsgDispatch;
244        if((it->second).size()<10) continue;
245 
246        Int_t idmacrocell=it->first;
247        vector<Int_t> fHittedmc=it->second;
248        vector<Double_t> xmc;
249        vector<Double_t> ymc;
250        vector<Double_t> tmc;
251            vector<Int_t> cmc;
252               
253        TH1F *h1 = 0;
254            if (fDebugPlots ) fListObj->Add(h1 = new TH1F("h1","h1",10,0,10));
255        for(Int_t i=0;i<(Int_t)fHittedmc.size();i++){
256            if (fDebugPlots) h1->Fill(fEv->GetRecoPixelData(fHittedmc[i])->GetCounts());
257            Int_t pixid=fEv->GetRecoPixelData(fHittedmc[i])->GetPixelId();
258                    TVector3 dummy2 = runpars->PixelCenter(pixid);
259            xmc.push_back(dummy2.x()/m);//m
260            ymc.push_back(dummy2.y()/m);//m
261            Int_t tmcgtu = fEv->GetRecoPixelData(fHittedmc[i])->GetGtu();
262            tmc.push_back(tmcgtu*fGtuLength /microsecond*0.01); // ????? (RP)
263            cmc.push_back(fEv->GetRecoPixelData(fHittedmc[i])->GetCounts());
264        }
265       
266        if (fDebugPlots) {
267            string pathname = Form("/HoughModule%d/", fIdEv);
268            gDirectory->cd(pathname.c_str());
269                string pathnames = Form("SingleHoughMacrocell%d", idmacrocell);
270            TDirectory *path = new TDirectory(pathnames.c_str(),pathnames.c_str());
271            path->cd();
272            h1->Draw();
273                h1->Write(); 
274        }
275            SafeDelete(h1);
276       
277       
278        Bool_t possiblesignal = HoughTransform(xmc,ymc,tmc,cmc);
279       
280        if ( possiblesignal ) {
281            Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably HAS signal!!!" << MsgDispatch;
282            fIdMacrocellsSelected.push_back(idmacrocell);
283                fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
284        } else {
285            Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably DOESN'T HAVE a signal,";
286                if ( idmacrocell==fMainMacroCellId ){
287                        Msg(EsafMsg::Info) << " but is added because is the main!" << MsgDispatch;
288                fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
289                    fIdMacrocellsSelected.push_back(idmacrocell);
290                } else {
291                        Msg(EsafMsg::Info) << " discarded!" << MsgDispatch;
292                }
293        }
294       
295        xmc.clear();
296        ymc.clear();
297        tmc.clear();
298        cmc.clear();
299    }
300   
301    if (fIdMacrocellsSelected.size()<1 || fHittedFinalMacrocell.size()<10)      return kFALSE;
302    Msg(EsafMsg::Info) << "Re-process " << fIdMacrocellsSelected.size() << " selected macrocells by ProcessHoughSingleMacrocell()" << MsgDispatch; 
303    Bool_t ret = ProcessHoughAll();
304
305    return ret;   
306}
307
308//_____________________________________________________________________________
309Bool_t HoughModule::ProcessVelocity(){
310    //
311    // Do histograms of velocities for each macrocell hitted and trying to select
312    // which are hitted by the real signal
313    //
314       
315    Msg(EsafMsg::Info) << " ...ProcessVelocity()" << MsgDispatch;
316   
317    Int_t numcountsmin = (Int_t)Conf()->GetNum("HoughModule.fMinCountsVelocity");
318       
319    RecoEventHeader header = fEv->GetHeader();
320    ERunParameters *runpars = header.GetRunPars();
321       
322    map< Int_t,vector<Int_t> > mapIdfHittedSingleCell;
323    for(Int_t i=0;i<200;i++) {
324        RecoCellInfo *info = fEv->GetRecoCellInfo(i);
325        if( info != NULL) {
326                vector<Int_t> empty;
327                pair< Int_t,vector<Int_t> > p(i,empty);
328                mapIdfHittedSingleCell.insert(p);
329        }
330    }
331    Msg(EsafMsg::Info) << "Number of macrocells hitted = " << mapIdfHittedSingleCell.size() << MsgDispatch;
332   
333    map< Int_t,vector<Int_t> >::iterator it;
334    for( Int_t i=0; i<fTotalNumPoints; i++ ) {
335        Int_t threshold(0);
336        if ( numcountsmin >=0 ) threshold = numcountsmin;
337        else {
338            Int_t uid = fEv->GetRecoPixelData(i)->GetPixelId();
339            Double_t rate = fEv->GetHeader().GetRunPars()->GetNightGlowRateByUId(uid)*fGtuLength/microsecond; 
340            threshold = (Int_t) Ceil(rate + Abs(numcountsmin)*Sqrt(rate));
341        }
342        if( fEv->GetRecoPixelData(i)->GetCounts() >= threshold  ){
343                Int_t cellid=fEv->GetRecoPixelData(i)->GetMacroCellId();
344            it = mapIdfHittedSingleCell.find(cellid);
345                it->second.push_back(i);
346            }
347    }
348   
349    Int_t binv = 20;
350    if (numcountsmin >= 0) {
351        Msg(EsafMsg::Info) << "Selected pixels with counts > " << numcountsmin << MsgDispatch;
352    } else {
353        Msg(EsafMsg::Info) << "Selectet pixels with counts over " << -numcountsmin << " sigma from background mean" << MsgDispatch;
354    }
355    for( it = mapIdfHittedSingleCell.begin(); it != mapIdfHittedSingleCell.end(); it++ ) {
356        Int_t entriesTH2(0);           
357        Msg(EsafMsg::Info) << "macrocellid = " << it->first << "\t with " << (it->second).size() << " points" << MsgDispatch;
358       
359        if((it->second).size()<10) continue;
360       
361        TH2F *hsxt(0), *hsyt(0), *s(0), *sxy(0);
362        TH1F *sx(0), *sy(0);
363        if (fDebugPlots) {
364            fListObj->Add(hsxt = new TH2F("hsxt","hsxt;t(mus);x(mm)",1000,0,200,1000,-1000,1000));
365                fListObj->Add(hsyt = new TH2F("hsyt","hsyt;t(mus);y(mm)",1000,0,200,1000,-1000,1000));
366        }
367        fListObj->Add(s = new TH2F("s","s;velocita;coeffang",200,0,2,200,-4,4)); 
368        fListObj->Add(sx = new TH1F("sx","sx;velocita",binv,-2,2)); 
369        fListObj->Add(sy = new TH1F("sy","sy;velocita",binv,-2,2)); 
370        fListObj->Add(sxy = new TH2F("sxy","velocita",binv,-2,2,binv,-2,2)); 
371        Int_t idmacrocell = it->first;
372        vector<Int_t> fHittedmc = it->second;
373       
374        for( Int_t i=0;i<(Int_t)fHittedmc.size();i++ ) {
375            Int_t pixid = fEv->GetRecoPixelData(fHittedmc[i])->GetPixelId();
376                TVector3 dummy2 = runpars->PixelCenter(pixid);//mm
377            Int_t tmcgtu = fEv->GetRecoPixelData(fHittedmc[i])->GetGtu();
378            Int_t counts = fEv->GetRecoPixelData(fHittedmc[i])->GetCounts();
379            if (fDebugPlots) {
380                hsxt->Fill(tmcgtu*fGtuLength/microsecond,dummy2.x(),counts);
381                    hsyt->Fill(tmcgtu*fGtuLength/microsecond,dummy2.y(),counts);
382            }
383            Int_t pointstocouple = fHittedmc.size();
384            for( Int_t j=i+1;j<pointstocouple;j++ ) {
385                Int_t pixtocouple = fEv->GetRecoPixelData(fHittedmc[j])->GetPixelId();
386                        TVector3 dummycouple = runpars->PixelCenter(pixtocouple);
387                        Int_t tmccouple = fEv->GetRecoPixelData(fHittedmc[j])->GetGtu();
388                    if (Abs(dummy2.x()-dummycouple.x())>0 && Abs(dummy2.y()-dummycouple.y())>0 && tmccouple-tmcgtu>0 
389                        && tmccouple-tmcgtu<20 && Abs(dummy2.x()-dummycouple.x())<100 && Abs(dummy2.y()-dummycouple.y())<100 ) {
390                            Int_t countproduct = (Int_t)Sqrt(1.e0*counts*fEv->GetRecoPixelData(fHittedmc[j])->GetCounts());
391                            Double_t m = (dummy2.y()-dummycouple.y())/(dummy2.x()-dummycouple.x());
392                        Double_t v = Sqrt((dummy2.y()-dummycouple.y())*(dummy2.y()-dummycouple.y())+(dummy2.x()-dummycouple.x())
393                            *(dummy2.x()-dummycouple.x()))/((tmccouple-tmcgtu)*fGtuLength/microsecond);// mm/us
394                            Double_t vx = (dummy2.x()-dummycouple.x())/((tmcgtu-tmccouple)*fGtuLength/microsecond);
395                            Double_t vy = (dummy2.y()-dummycouple.y())/((tmcgtu-tmccouple)*fGtuLength/microsecond);
396                            entriesTH2++;
397                    s->Fill(v,m,countproduct);
398                                sx->Fill(vx,countproduct);
399                        sy->Fill(vy,countproduct);
400                        sxy->Fill(vx,vy,countproduct);
401                        }
402            }
403        }   
404   
405        Int_t binmaxv = sxy->GetMaximumBin();
406        Int_t nxv = (sxy->GetXaxis())->GetNbins()+2;
407        Int_t binxmaxv = binmaxv % nxv;
408        Int_t binymaxv = binmaxv / nxv;
409        Float_t vxmax[1], vymax[1];
410        vxmax[0] = (sxy->GetXaxis())->GetBinCenter(binxmaxv);
411        vymax[0] = (sxy->GetYaxis())->GetBinCenter(binymaxv);
412        TGraph *gvmax=new TGraph(1,vxmax,vymax);
413        gvmax->SetMarkerColor(2); gvmax->SetMarkerStyle(29); gvmax->SetMarkerSize(1.6);
414        //Float_t contmaxv=sxy->GetBinContent(binxmaxv,binymaxv);
415        Stat_t countx(0),county(0);
416        for(Int_t bin=1;bin<binv;bin++){
417            countx += sx->GetBinContent(bin);
418                county += sy->GetBinContent(bin);
419        }
420        Float_t countmaxx = sx->GetMaximum();
421        Float_t countmaxy = sy->GetMaximum();
422        Float_t countmeanx = (Float_t)((Float_t)countx/binv);
423        Float_t countmeany = (Float_t)((Float_t)county/binv);
424        Float_t ms = 20.;
425        /*cout<<"countmeanx="<<countmeanx<<"\tcontmaxx="<<countmaxx<<"\tcontmaxx-countmeanx="<<countmaxx-countmeanx<<"\tsig="<<ms*TMath::Sqrt(countmeanx)<<endl;
426        cout<<"countmeany="<<countmeany<<"\tcontmaxy="<<countmaxy<<"\tcontmaxy-countmeany="<<countmaxy-countmeany<<"\tsig="<<ms*TMath::Sqrt(countmeany)<<endl;*/
427        string a = "a"; string b = "b";
428        string pathname = Form("/HoughModule%d/", fIdEv);
429        if (fDebugPlots) gDirectory->cd(pathname.c_str());
430        string pathnames;
431        if( ((countmaxx-countmeanx) > ms*Sqrt(countmeanx) || (countmaxy-countmeany) > ms*Sqrt(countmeany)) && sx->GetEntries()>10 ) {
432            fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
433            fIdMacrocellsSelected.push_back(idmacrocell);
434            Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably HAS signal!!!" << MsgDispatch;
435            if (idmacrocell==fMainMacroCellId) {
436                pathnames = Form("VelocityMacrocellmain%d", idmacrocell);
437            } else {
438                pathnames = Form("VelocityMacrocell%d", idmacrocell);
439            }
440        } else{
441            Msg(EsafMsg::Info) << "Macrocell " << idmacrocell << " probably DOESN'T HAVE a signal,";
442            if (idmacrocell==fMainMacroCellId) {
443                pathnames = Form("VelocityMacrocellmain%d", idmacrocell);
444                Msg(EsafMsg::Info) <<" but is added because is the main!"<< MsgDispatch;
445                fHittedFinalMacrocell.insert(fHittedFinalMacrocell.end(),fHittedmc.begin(),fHittedmc.end());
446                fIdMacrocellsSelected.push_back(idmacrocell);
447            } else {
448                pathnames = Form("VelocityMacrocelldiscarded%d", idmacrocell);
449                Msg(EsafMsg::Info) <<" discarded!"<< MsgDispatch;
450            }
451        }
452               
453        s->SetTitle(pathnames.c_str());
454        if (fDebugPlots) {
455            TDirectory *path = new TDirectory(pathnames.c_str(),pathnames.c_str());
456            path->cd();
457            TCanvas *ac;
458            fListObj->Add(ac = new TCanvas(a.c_str(),a.c_str(),1)); ac->Divide(2,2);
459            ac->cd(1);s->Draw();    ac->cd(2);s->Draw("LEGO");
460            ac->cd(3);sx->Draw();       ac->cd(4);sy->Draw(); 
461            ac->Write();
462            SafeDelete(ac);
463            TCanvas *bc;
464            fListObj->Add(bc = new TCanvas(b.c_str(),b.c_str(),1)); bc->Divide(2,2);
465            bc->cd(1);hsxt->Draw(); 
466            bc->cd(2);hsyt->Draw(); 
467            bc->cd(3);sxy->Draw();gvmax->Draw("same"); 
468            bc->cd(4);sxy->Draw("LEGO"); 
469            bc->Write();
470            SafeDelete(bc);
471        }
472        SafeDelete(s);  SafeDelete(sx); SafeDelete(sy);
473        SafeDelete(hsxt);       SafeDelete(hsyt);
474        SafeDelete(sxy); SafeDelete(gvmax);
475    }
476   
477    if (fIdMacrocellsSelected.size()<1 || fHittedFinalMacrocell.size()<10)      return kFALSE;
478       
479    Msg(EsafMsg::Info) << "Re-process " << fIdMacrocellsSelected.size() << " selected macrocells by ProcessVelocity()" << MsgDispatch; 
480    Bool_t ret = ProcessHoughAll();
481   
482    return ret;
483
484}
485
486//_____________________________________________________________________________
487Bool_t HoughModule::ProcessHoughAll(){
488    //
489    // Doing Hough Transform to all macrocells selected by previous
490    // methods or by euso trigger
491    //
492       
493    Msg(EsafMsg::Info) << " ...ProcessHoughAll()" << MsgDispatch;
494       
495    if (fDebugPlots) {
496        string pathname = Form("/HoughModule%d/", fIdEv);
497        gDirectory->cd(pathname.c_str());
498        string pathnametotal = Form("HoughAll");
499        TDirectory *pathtotal = new TDirectory(pathnametotal.c_str(),pathnametotal.c_str());
500        pathtotal->cd();
501    }
502       
503    RecoEventHeader header = fEv->GetHeader();
504    ERunParameters *runpars = header.GetRunPars();
505       
506        fNumHitsMinimum = (Int_t)Conf()->GetNum("HoughModule.fMinCountsHoughAll");
507    if ( fIsAlreadyPattReco && fNumHitsMinimum < 0 ){
508        Msg(EsafMsg::Info) << "Using relative threshold at " << Abs(fNumHitsMinimum) << " sigma from the mean background value" << MsgDispatch; 
509    }else if (!fIsAlreadyPattReco ) {
510        Msg(EsafMsg::Info) << "Let's try beginning with MinCountsHoughAll = " << fNumHitsMinimum << MsgDispatch;
511    }
512   
513    Int_t regr = 0;
514   
515    Int_t gtumin = 10000;
516    Int_t gtumax = -1;
517    Int_t deltagtu = 0;
518    Int_t mean_cmin = 0; // mean number of hits minimum in pixels (differs from fNumHitsMinimum using relative threshold)
519   
520    fNumPointsAll = 0;
521    fHittedAll.clear();
522    if(!fIsAlreadyPattReco && !(fHittedFinalMacrocell.size())){
523            Msg(EsafMsg::Info) << "processing all macrocells!" << MsgDispatch;
524        mean_cmin = 0;
525            for( Int_t i=0; i<fTotalNumPoints; i++ ) {
526                if (fEv->GetRecoPixelData(i)->GetGtu()<gtumin) gtumin=fEv->GetRecoPixelData(i)->GetGtu();
527            if (fEv->GetRecoPixelData(i)->GetGtu()>gtumax) gtumax=fEv->GetRecoPixelData(i)->GetGtu();
528            Int_t countsmin = 0;
529            // compute relative threshold if requested
530            if ( fNumHitsMinimum >= 0 ) {
531                countsmin = fNumHitsMinimum;
532                mean_cmin = fNumHitsMinimum;
533            } else {
534                Int_t uid = fEv->GetRecoPixelData(i)->GetPixelId();
535                Double_t rate = fEv->GetHeader().GetRunPars()->GetNightGlowRateByUId(uid) * fGtuLength/microsecond;
536                countsmin = (Int_t)Ceil( rate + Abs(fNumHitsMinimum)*Sqrt(rate) );
537            }
538            if (fEv->GetRecoPixelData(i)->GetCounts() >= countsmin) {
539                fHittedAll.push_back(i);
540                if ( fNumHitsMinimum < 0 ) mean_cmin += countsmin;
541            }
542        }
543        if ( !fHittedAll.size() ) {
544            Msg(EsafMsg::Warning) << "No pixels over threshold. Exit from this event." << MsgDispatch;
545            fHittedAll.clear();
546            fMainMacroCellId = 0;
547            fNumPointsAll = 0;
548            return kFALSE;
549        }
550        if ( fNumHitsMinimum < 0 ) mean_cmin /= (Int_t)fHittedAll.size();
551    } else if (!fIsAlreadyPattReco){
552        Msg(EsafMsg::Info)<<"processing fHittedFinalMacrocell.size()="<<fHittedFinalMacrocell.size()<<MsgDispatch;
553        mean_cmin = 0;
554        Int_t add = 0;
555        for( Int_t j=0; j<(Int_t)(fHittedFinalMacrocell.size()); j++ ){
556            Int_t i = fHittedFinalMacrocell[j];
557            if (fEv->GetRecoPixelData(i)->GetGtu()<gtumin) gtumin=fEv->GetRecoPixelData(i)->GetGtu();
558            if (fEv->GetRecoPixelData(i)->GetGtu()>gtumax) gtumax=fEv->GetRecoPixelData(i)->GetGtu();
559            Int_t countsmin = 0;
560            // compute relative threshold if requested
561            if ( fNumHitsMinimum>=0 ) {
562                countsmin = fNumHitsMinimum;
563                mean_cmin = fNumHitsMinimum;
564            } else {
565                Int_t uid = fEv->GetRecoPixelData(i)->GetPixelId();
566                Double_t rate = fEv->GetHeader().GetRunPars()->GetNightGlowRateByUId(uid) * fGtuLength/microsecond;
567                countsmin = (Int_t)Ceil( rate + Abs(fNumHitsMinimum)*Sqrt(rate) );
568                //mean_cmin += countsmin;
569            }
570            if (fEv->GetRecoPixelData(i)->GetCounts() >= countsmin) {
571                fHittedAll.push_back(i);
572                add++;
573                if ( fNumHitsMinimum < 0 ) mean_cmin += countsmin;
574            }
575        }
576        if ( fNumHitsMinimum < 0 && add ) mean_cmin /= add;
577    } else if (fIsAlreadyPattReco) {
578        fHittedAll = fPreviousPattReco;
579        fNumHitsMinimum = fEv->FindGlobalData("PatternRecognition")->GetInt("NumPointsMinimum");
580        mean_cmin = fEv->FindGlobalData("PatternRecognition")->GetInt("MeanThreshold");
581        for (Int_t i(0); i<(Int_t)fHittedAll.size(); i++) {
582            Int_t j = fHittedAll[i];
583            if (fEv->GetRecoPixelData(j)->GetGtu()<gtumin) gtumin=fEv->GetRecoPixelData(j)->GetGtu();
584            if (fEv->GetRecoPixelData(j)->GetGtu()>gtumax) gtumax=fEv->GetRecoPixelData(j)->GetGtu();
585        }
586    }
587       
588   
589    fNumPointsAll = fHittedAll.size();
590       
591    fTmin = gtumin*fGtuLength/microsecond*0.01; //???? (RP)
592    fTmax = gtumax*fGtuLength/microsecond*0.01; //???? (RP)
593    deltagtu = gtumax - gtumin;
594   
595    Msg(EsafMsg::Info) << "Num. of pixels over threshold= " << fNumPointsAll << " over fTotalNumPoints="  << fTotalNumPoints << MsgDispatch;
596    if (fNumPointsAll<10 || deltagtu<5) {
597        if( fNumPointsAll<10 )   Msg(EsafMsg::Warning) << "Number of points "<< fNumPointsAll << " is less than 10, Exit from this event." << MsgDispatch;
598        if ( deltagtu<5 )       Msg(EsafMsg::Warning) << " deltagtu = "<< deltagtu << " is less than 5, Exit from this event." << MsgDispatch;
599        fHittedAll.clear();
600        fMainMacroCellId = 0;
601        fNumPointsAll = 0;
602        return kFALSE;
603    }
604   
605    vector<Double_t> x;
606    vector<Double_t> y;
607    vector<Double_t> t;
608    vector<Int_t> c;
609       
610    TH2F *single2(0), *single2xt(0), *single2yt(0);
611    Int_t bin = 2000;
612    if (fDebugPlots) {
613        fListObj->Add(single2 = new TH2F("single2","single2;x(m);y(m)",bin,-1,1,bin,-1,1));
614        fListObj->Add(single2xt = new TH2F("single2xt","single2xt;t;x(m)",bin,fTmin,fTmax,bin,-1,1));
615            fListObj->Add(single2yt = new TH2F("single2yt","single2yt;t;y(m)",bin,fTmin,fTmax,bin,-1,1));
616    }
617       
618        for( Int_t i=0; i < fNumPointsAll; i++ ) {
619        Int_t pixid=fEv->GetRecoPixelData(fHittedAll[i])->GetPixelId();
620                TVector3 dummy2 = runpars->PixelCenter(pixid); // mm
621        x.push_back(dummy2.x()/m);
622        y.push_back(dummy2.y()/m);
623        t.push_back((Double_t)fEv->GetRecoPixelData(fHittedAll[i])->GetGtu()*fGtuLength/microsecond*0.01); // ??? (RP)
624        c.push_back(fEv->GetRecoPixelData(fHittedAll[i])->GetCounts());
625        if (fDebugPlots) {
626            single2->Fill(x[i],y[i],c[i]);
627            single2xt->Fill(t[i],x[i],c[i]);
628            single2yt->Fill(t[i],y[i],c[i]);
629        }
630    }
631   
632    if (fDebugPlots) {
633        TCanvas *HMinit;
634        fListObj->Add(HMinit = new TCanvas("HMinit","HMinit",1));
635        HMinit->Divide(3,1);
636        HMinit->cd(1);  single2->SetMarkerStyle(6);  single2->Draw();
637        HMinit->cd(2);  single2xt->SetMarkerStyle(6);  single2xt->Draw();
638        HMinit->cd(3);  single2yt->SetMarkerStyle(6);  single2yt->Draw();
639        HMinit->Write();   
640        SafeDelete(HMinit);
641    }
642    SafeDelete(single2xt);
643    SafeDelete(single2yt);
644    SafeDelete(single2);
645   
646    HoughTransform(x,y,t,c); 
647   
648    x.clear();  y.clear();  t.clear();  c.clear();
649
650    if( (mean_cmin - regr) > 3) regr++;
651   
652    vector<Double_t> xsel;
653    vector<Double_t> ysel;
654    vector<Double_t> tsel;
655    vector<Int_t> csel;
656   
657    TH2F *single2sel(0), *single2xtsel(0), *single2ytsel(0);
658    if (fDebugPlots) {
659        fListObj->Add(single2sel = new TH2F("single2sel","single2sel;x(m);y(m)",bin,fXmin,fXmax,bin,fYmin,fYmax));
660        fListObj->Add(single2xtsel = new TH2F("single2xtsel","single2xtsel;t;x(m)",bin,fTmin,fTmax,bin,fXmin,fXmax));
661        fListObj->Add(single2ytsel = new TH2F("single2ytsel","single2ytsel;t;y(m)",bin,fTmin,fTmax,bin,fYmin,fYmax));
662    }
663    for( Int_t i=0; i<fTotalNumPoints; i++ ) {
664        Int_t pixid=fEv->GetRecoPixelData(i)->GetPixelId();
665        TVector3 dummy2 = runpars->PixelCenter(pixid);
666        Double_t x = dummy2.x()/m; // meters
667        Double_t y = dummy2.y()/m; // meters
668        Double_t t=(Double_t)fEv->GetRecoPixelData(i)->GetGtu()*fGtuLength/microsecond*0.01; // ???
669        Int_t c=fEv->GetRecoPixelData(i)->GetCounts();
670        if ( x>=fXmin && x<=fXmax && y>=fYmin && y<=fYmax && c>=(mean_cmin-regr) ) {
671            if (fDebugPlots) {
672                single2sel->Fill(x,y,c);
673                single2xtsel->Fill(t,x,c);
674                single2ytsel->Fill(t,y,c);
675            }
676            xsel.push_back(x);
677            ysel.push_back(y);
678            tsel.push_back(t);
679            csel.push_back(c);
680            fHittedSelHough.push_back(i);
681        }
682    }
683   
684    Msg(EsafMsg::Info) << "first selected zone: pixels with at least " << mean_cmin-regr << " hits = " << fHittedSelHough.size() << MsgDispatch;
685    if(fHittedSelHough.size()<10 && !fIsAlreadyPattReco){
686        Msg(EsafMsg::Info) << "too few! exit from this event" <<MsgDispatch;
687        return kFALSE;
688    } else if (fHittedSelHough.size()<10 && fIsAlreadyPattReco) {
689        fHittedFinal = fHittedFinal;
690        fFinalNumHitsMinimum = mean_cmin - regr;
691        return kTRUE;
692    } 
693   
694    if (fDebugPlots) {
695        TCanvas *HMsel;
696        fListObj->Add(HMsel = new TCanvas("HMsel","HMsel",1));
697        HMsel->Divide(3,1);
698        HMsel->cd(1);  single2sel->SetMarkerStyle(6);  single2sel->Draw();
699        HMsel->cd(2);  single2xtsel->SetMarkerStyle(6);  single2xtsel->Draw();
700        HMsel->cd(3);  single2ytsel->SetMarkerStyle(6);  single2ytsel->Draw();
701        HMsel->Write();
702        SafeDelete(HMsel);
703    }
704    SafeDelete(single2xtsel);
705    SafeDelete(single2ytsel);
706    SafeDelete(single2sel);
707   
708    HoughTransform(xsel,ysel,tsel,csel); 
709   
710    xsel.clear();  ysel.clear();  tsel.clear();  csel.clear();
711    if(mean_cmin-regr > 2)      regr++;
712    vector<Double_t> xselsel;
713    vector<Double_t> yselsel;
714    vector<Double_t> tselsel;
715    vector<Int_t> cselsel;
716   
717    TH2F *single2selsel(0), *single2xtselsel(0), *single2ytselsel(0);
718    if (fDebugPlots) { 
719        fListObj->Add(single2selsel = new TH2F("single2selsel","single2selsel;x(m);y(m)",bin,fXmin,fXmax,bin,fYmin,fYmax));
720        fListObj->Add(single2xtselsel = new TH2F("single2xtselsel","single2xtselsel;t;x(m)",bin,fTmin,fTmax,bin,fXmin,fXmax));
721        fListObj->Add(single2ytselsel = new TH2F("single2ytselsel","single2ytselsel;t;y(m)",bin,fTmin,fTmax,bin,fYmin,fYmax));
722    }
723    for( Int_t i=0; i<fTotalNumPoints; i++ ) {
724        Int_t pixid=fEv->GetRecoPixelData(i)->GetPixelId();
725            TVector3 dummy2 = runpars->PixelCenter(pixid);
726        Double_t x = dummy2.x()/m; //meters
727        Double_t y = dummy2.y()/m; //meters
728        Double_t t = (Double_t)fEv->GetRecoPixelData(i)->GetGtu()*fGtuLength/microsecond*0.01;
729        Int_t c=fEv->GetRecoPixelData(i)->GetCounts();
730        if ( x>=fXmin && x<=fXmax && y>=fYmin && y<=fYmax && c>=(mean_cmin-regr) ) {   
731            if (fDebugPlots) {
732                single2selsel->Fill(x,y,c);
733                single2xtselsel->Fill(t,x,c);
734                single2ytselsel->Fill(t,y,c);
735            }
736                xselsel.push_back(x);
737                yselsel.push_back(y);
738                tselsel.push_back(t);
739                cselsel.push_back(c);
740                fHittedSel2Hough.push_back(i);
741        }
742    }
743    Msg(EsafMsg::Info) << "second selected zone: pixels with at least " << mean_cmin-regr << " hits = " << fHittedSel2Hough.size()<< MsgDispatch;
744    if(fHittedSel2Hough.size()<10 && !fIsAlreadyPattReco){
745        Msg(EsafMsg::Info) << "too few! exit from this event" <<MsgDispatch;
746        return kFALSE;
747    } else if (fHittedSel2Hough.size()<10 && fIsAlreadyPattReco) {
748        fHittedFinal = fHittedSelHough;
749        fFinalNumHitsMinimum = mean_cmin - regr;
750        return kTRUE;
751    } 
752   
753    if (fDebugPlots) {
754        TCanvas *HMselsel;
755        fListObj->Add(HMselsel = new TCanvas("HMselsel","HMselsel",1));
756        HMselsel->Divide(3,1);
757        HMselsel->cd(1);  single2selsel->SetMarkerStyle(6);  single2selsel->Draw();
758        HMselsel->cd(2);  single2xtselsel->SetMarkerStyle(6);  single2xtselsel->Draw();
759        HMselsel->cd(3);  single2ytselsel->SetMarkerStyle(6);  single2ytselsel->Draw();
760        HMselsel->Write();
761        SafeDelete(HMselsel);
762    }
763    SafeDelete(single2xtselsel);
764    SafeDelete(single2ytselsel);
765    SafeDelete(single2selsel);
766   
767   
768    HoughTransform(xselsel,yselsel,tselsel,cselsel,fHittedSel2Hough); 
769    fHittedFinal = fHittedHoughMerged;
770    xselsel.clear();  yselsel.clear();  tselsel.clear();  cselsel.clear();
771   
772    //now not used
773    /*vector<Double_t> xend;
774    vector<Double_t> yend;
775    vector<Double_t> tend;
776    vector<Int_t> cend;*/
777   
778    TH2F *single2end(0), *single2xtend(0), *single2ytend(0);
779    if (fDebugPlots) {
780        fListObj->Add(single2end = new TH2F("single2end","single2end;x(m);y(m)",bin,fXmin,fXmax,bin,fYmin,fYmax));
781        fListObj->Add(single2xtend = new TH2F("single2xtend","single2xtend;t;x(m)",bin,fTmin,fTmax,bin,fXmin,fXmax));
782        fListObj->Add(single2ytend = new TH2F("single2ytend","single2ytend;t;y(m)",bin,fTmin,fTmax,bin,fYmin,fYmax));
783        for( Int_t i=0; i<(Int_t)fHittedHoughMerged.size(); i++ ) {
784        Int_t pixid=fEv->GetRecoPixelData(fHittedHoughMerged[i])->GetPixelId();
785                TVector3 dummy2 = runpars->PixelCenter(pixid);
786        Double_t x = dummy2.x()/m; // meters
787        Double_t y = dummy2.y()/m; // meters
788        Double_t t = (Double_t)fEv->GetRecoPixelData(fHittedHoughMerged[i])->GetGtu()*fGtuLength/microsecond*0.01; // ???
789        Int_t c = fEv->GetRecoPixelData(fHittedHoughMerged[i])->GetCounts();
790        single2end->Fill(x,y,c);
791        single2xtend->Fill(t,x,c);
792        single2ytend->Fill(t,y,c);
793        /*xend.push_back(x);
794        yend.push_back(y);
795        tend.push_back(t);
796        cend.push_back(c);*/
797    }
798    }
799    Msg(EsafMsg::Info) << "final selected zone: pixels with at least " << mean_cmin-regr << " hits = " << fHittedHoughMerged.size()<< MsgDispatch;
800    if(fHittedHoughMerged.size()<10 && !fIsAlreadyPattReco){
801        Msg(EsafMsg::Info) << "too few! exit from this event" <<MsgDispatch;
802        return kFALSE;
803    } else if (fHittedHoughMerged.size()<10 && fIsAlreadyPattReco) {
804        fHittedFinal = fHittedSel2Hough;
805        fFinalNumHitsMinimum = mean_cmin - regr;
806        return kTRUE;
807    } 
808   
809    if (fDebugPlots) {
810        TCanvas *HMfinal;
811        fListObj->Add(HMfinal = new TCanvas("HMfinal","HMfinal",1));
812        HMfinal->Divide(3,1);
813        HMfinal->cd(1);  single2end->SetMarkerStyle(6);  single2end->Draw();
814        HMfinal->cd(2);  single2xtend->SetMarkerStyle(6);  single2xtend->Draw();
815        HMfinal->cd(3);  single2ytend->SetMarkerStyle(6);  single2ytend->Draw();
816        HMfinal->Write();
817        SafeDelete(HMfinal);
818    }
819    SafeDelete(single2xtend);
820    SafeDelete(single2ytend);
821    SafeDelete(single2end);
822   
823    //xend.clear(); yend.clear(); tend.clear(); cend.clear();
824    fFinalNumHitsMinimum = mean_cmin - regr;
825    return kTRUE;
826   
827}
828
829
830//_____________________________________________________________________________
831Bool_t HoughModule::HoughTransform(vector<Double_t> x,vector<Double_t> y,vector<Double_t> t,vector<Int_t> c,vector<Int_t> id){
832    //
833    // Implements the Hough Transform
834    // Return true if there is a possible signal
835    //
836
837    if( x.size()<10 )   return 0;
838       
839    Int_t optionselection = 1;
840       
841    Float_t ex = 0.004;
842    Float_t ey = 0.004;
843    Float_t et = fGtuLength/microsecond*0.01; // ???
844       
845    HoughFit *xtFit = new HoughFit(t, x, c, optionselection,et,ex,id);
846    //Double_t vx = xtFit->GetSlope();//Double_t wx = xtFit->GetOffset();
847    fXmin = xtFit->GetMinY();
848    fXmax = xtFit->GetMaxY();
849    Double_t tmin1 = xtFit->GetMinX();
850    Double_t tmax1 = xtFit->GetMaxX();
851    fAbsDevX = xtFit->GetAbsDev();
852    Bool_t possiblesignalxt = xtFit->IsSignal();
853    fCmaxToCmeanX = xtFit->GetMaxToMedium();
854    vector<Int_t> idselx = xtFit->GetIdSel();
855    delete xtFit;
856       
857    HoughFit *ytFit = new HoughFit(t, y, c, optionselection,et,ey,id);
858    //Double_t vy = ytFit->GetSlope();//Double_t wy = ytFit->GetOffset();
859    fYmin = ytFit->GetMinY();
860    fYmax = ytFit->GetMaxY();
861    Double_t tmin2 = ytFit->GetMinX();
862    Double_t tmax2 = ytFit->GetMaxX();
863    fAbsDevY = ytFit->GetAbsDev();
864    Bool_t possiblesignalyt = ytFit->IsSignal();
865    fCmaxToCmeanY = ytFit->GetMaxToMedium();
866    vector<Int_t> idsely = ytFit->GetIdSel();
867    delete ytFit;
868       
869    if(id.size()>0) {
870        Int_t sumsel=idsely.size()+idselx.size();
871        //cout<<"id.size()="<<id.size()<<"\tidselx.size()="<<idselx.size()<<"\tidsely.size()="<<idsely.size()<<endl;
872    vector<Int_t> idselectedmerged(sumsel);
873        merge(idselx.begin(), idselx.end(),idsely.begin(), idsely.end(),idselectedmerged.begin());
874        vector<Int_t>::iterator diversi=unique(idselectedmerged.begin(),idselectedmerged.end());
875        idselectedmerged.erase(diversi,idselectedmerged.end());
876        //cout<<"idselectedmerged.size()="<<idselectedmerged.size()<<endl;
877        fHittedHoughMerged = idselectedmerged;
878    }
879       
880    //Msg(EsafMsg::Info) <<"possiblesignalxt="<<possiblesignalxt<<"\tpossiblesignalyt="<<possiblesignalyt<<MsgDispatch;
881    //Msg(EsafMsg::Info) <<"fCmaxToCmeanX="<<fCmaxToCmeanX<<"\tcmaxtocmedioy="<<fCmaxToCmeanY<<MsgDispatch;
882    Bool_t possible = kFALSE;
883    if( possiblesignalxt || possiblesignalyt ) possible = kTRUE;
884       
885    fTmin = (tmin1<tmin2) ? tmin1 : tmin2;
886    fTmax = (tmax1<tmax2) ? tmax2 : tmax1;
887       
888    x.clear();  y.clear();  t.clear();  c.clear();
889   
890    return possible;
891   
892}
893
894//_____________________________________________________________________________
895Bool_t HoughModule::PostProcess() {
896    //
897    // post-process
898    //
899   
900    if ( !fTotalNumPoints ) return kTRUE;
901   
902    TmpMemoryClean();
903   
904    fEv = (RecoEvent*)0;
905    return kTRUE;
906       
907}
908
909//____________________________________________________________________________
910Bool_t HoughModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
911    //
912    // Save reco root data
913    //
914   
915    return kTRUE;
916}
917
918//____________________________________________________________________________
919Bool_t HoughModule::Done() {
920    //
921    // done
922    //
923   
924    fHittedAll.clear();
925    fHittedFinal.clear();
926    fHittedSelHough.clear();
927    fHittedSel2Hough.clear();
928    fIdMacrocellsSelected.clear();
929   
930    fClusters.clear();
931    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
932    return kTRUE;
933}
934
935//_____________________________________________________________________________
936void HoughModule::UserMemoryClean() {
937    //
938    // memory clean
939    //
940
941    Msg(EsafMsg::Debug) << "UserMemoryClean()" << MsgDispatch;
942    //fListObj->Delete();
943    //fListObj->Clear();
944       
945//     Int_t fHittedFinalsize = fHittedFinal.size();
946
947//     if(fHittedFinalsize<10 || !fTotalNumPoints || fHittedFinalsize>100000)           return;
948       
949    vector<EusoCluster*> *pclu = (vector<EusoCluster*>*)MyData()->GetObj("Clusters");
950    if (pclu) {
951        EusoCluster *dummy(NULL);
952        for(size_t i(0); i<pclu->size(); i++) {
953            dummy = (*pclu)[i];
954            delete dummy;
955        }
956        pclu->clear();
957        MyData()->RemoveObj("Clusters");
958        vector<EusoCluster*> dummyV;
959        dummyV.swap(*pclu);
960
961    }
962   
963    vector<Int_t> *pix = (vector<Int_t>*)MyData()->GetObj("CluPixels");
964    if (pix) {
965        pix->clear();
966        MyData()->RemoveObj("CluPixels");
967        vector<Int_t> dummyV;
968        dummyV.swap(*pix);
969    }
970   
971}
972
973
974//_____________________________________________________________________________
975void HoughModule::TmpMemoryClean() {
976    //
977    // Clean all the temporary buffers
978    //
979   
980    vector<Int_t> dummy1;
981    dummy1.swap(fHittedAll);
982    vector<Int_t> dummy2;
983    dummy2.swap(fHittedSelHough);
984    vector<Int_t> dummy3;
985    dummy3.swap(fHittedSel2Hough);
986    vector<Int_t> dummy4;
987    dummy4.swap(fHittedFinalMacrocell);
988    vector<Int_t> dummy5;
989    dummy5.swap(fIdMacrocellsSelected);
990    vector<Int_t> dummy6;
991    dummy6.swap(fHittedHoughMerged);
992
993}
994
995//_____________________________________________________________________________
996void HoughModule::WriteSingleCluster() {
997    //
998    // close and write the current cluster if is significant
999    //
1000   
1001    if ( !fHittedFinal.size() ) return; 
1002    // write the cluster
1003    EusoCluster *cl = new EusoCluster();
1004    for( UInt_t i=0; i<fHittedFinal.size(); i++ ) {
1005        cl->AddPoint( fHittedFinal[i] );
1006    }
1007    fClusters.push_back(cl);
1008    Msg(EsafMsg::Info) << "SingleCluster saved with "<< fHittedFinal.size() << " points" << MsgDispatch;
1009   
1010    // add data to event
1011    vector<Int_t> *pix = &fHittedFinal;
1012    MyData()->Add("CluPixels", pix);
1013    vector<EusoCluster*> *pclu = &fClusters;
1014    MyData()->Add("Clusters", pclu);
1015    MyData()->Add("NumHitsMinimum",fFinalNumHitsMinimum);
1016
1017    // add global data to the event
1018    RecoGlobalData* data = fEv->FindCreateGlobalData("PatternRecognition");
1019    data->CleanMaps();
1020    data->SetIssuer(GetName());
1021    data->Add("SelectedPixels",&fHittedFinal);
1022
1023}
1024
1025//_____________________________________________________________________________
1026EusoCluster* HoughModule::GetCluster( Int_t i ) {
1027    //
1028    // get a specified cluster
1029    //
1030
1031    if ( i >= 0 && i < (Int_t)fClusters.size() )
1032        return fClusters[i];
1033    else
1034        return NULL;
1035}
1036
1037
1038
Note: See TracBrowser for help on using the repository browser.