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 | |
---|
47 | using namespace TMath; |
---|
48 | using namespace EConst; |
---|
49 | using namespace sou; |
---|
50 | |
---|
51 | ClassImp(HoughModule) |
---|
52 | |
---|
53 | // FIXME memory leaks to be fixed asap. D.Naumov |
---|
54 | |
---|
55 | //_____________________________________________________________________________ |
---|
56 | HoughModule::HoughModule():RecoModule("HoughTransform") { |
---|
57 | // |
---|
58 | // ctor |
---|
59 | // |
---|
60 | } |
---|
61 | |
---|
62 | //_____________________________________________________________________________ |
---|
63 | HoughModule::~HoughModule() { |
---|
64 | // |
---|
65 | // dtor |
---|
66 | // |
---|
67 | } |
---|
68 | |
---|
69 | //_____________________________________________________________________________ |
---|
70 | Bool_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 | //_____________________________________________________________________________ |
---|
90 | Bool_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 | //_____________________________________________________________________________ |
---|
115 | Bool_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 | //_____________________________________________________________________________ |
---|
197 | Bool_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 | //_____________________________________________________________________________ |
---|
309 | Bool_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 | //_____________________________________________________________________________ |
---|
487 | Bool_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 | //_____________________________________________________________________________ |
---|
831 | Bool_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 | //_____________________________________________________________________________ |
---|
895 | Bool_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 | //____________________________________________________________________________ |
---|
910 | Bool_t HoughModule::SaveRootData(RecoRootEvent *fRecoRootEvent) { |
---|
911 | // |
---|
912 | // Save reco root data |
---|
913 | // |
---|
914 | |
---|
915 | return kTRUE; |
---|
916 | } |
---|
917 | |
---|
918 | //____________________________________________________________________________ |
---|
919 | Bool_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 | //_____________________________________________________________________________ |
---|
936 | void 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 | //_____________________________________________________________________________ |
---|
975 | void 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 | //_____________________________________________________________________________ |
---|
996 | void 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 | //_____________________________________________________________________________ |
---|
1026 | EusoCluster* 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 | |
---|