FairRoot/PandaRoot
PndTrkGemCombinatorial.cxx
Go to the documentation of this file.
1 //
3 // PndTrkGemCombinatorial
4 //
5 // Class for GEM combinatorial suppression
6 //
7 // authors: Lia Lavezzi - University of Torino (2014)
8 //
10 
11 #include "PndTrkGemCombinatorial.h"
12 
13 #include "PndTrkTools.h"
14 
15 // fairroot
16 #include "FairRootManager.h"
17 #include "FairRunAna.h"
18 #include "FairRuntimeDb.h"
19 // ROOT
20 #include "TClonesArray.h"
21 #include "TVector3.h"
22 #include "TMarker.h"
23 
24 #include "PndGemHit.h"
25 #include "PndGemMCPoint.h"
26 
27 #include <iostream>
28 
29 using namespace std;
30 
31 int nhits = 0, npoints = 0, nassocorrect = 0, nassowrong = 0, ndontusecorrect = 0, ndontusewrong = 0, nasso = 0, nunasso = 0;
32 
33 // ----- Default constructor -------------------------------------------
34 PndTrkGemCombinatorial::PndTrkGemCombinatorial(): fGemHitArray(NULL), fGemPointArray(NULL), fDisplayOn(kFALSE), fMCEval(kFALSE) {}
35 
36 PndTrkGemCombinatorial::PndTrkGemCombinatorial(TClonesArray *gemhitarray, int verbose) : fGemHitArray(gemhitarray), fGemPointArray(NULL), fDisplayOn(kFALSE), fMCEval(kFALSE), fVerbose(verbose) {
37 }
38 
39 // -------------------------------------------------------------------------
40 
41 // ----- Destructor ----------------------------------------------------
42 PndTrkGemCombinatorial::~PndTrkGemCombinatorial() {
43 
44  delete fGemHitArray;
45  delete fGemPointArray;
46  delete hxy;
47  delete hxy1;
48  delete hxy2;
49  delete hxy3;
50  delete display;
51 
52 }
53 // -------------------------------------------------------------------------
54 
55 
56 
57 // ----- Public method Init --------------------------------------------
58 void PndTrkGemCombinatorial::SwitchOnDisplay() {
59  fDisplayOn = kTRUE;
60  display = new TCanvas("display", "display", 0, 0, 800, 800); // CHECK
61  display->Divide(2, 2);
62  hxy = new TH2F("hxy", "xy plane", 130, -65, 65, 130, -65, 65);
63  hxy1 = new TH2F("hxy1", "xy plane 1", 76, -38, 38, 76, -38, 38);
64  hxy2 = new TH2F("hxy2", "xy plane 2", 94, -47, 47, 94, -47, 47);
65  hxy3 = new TH2F("hxy3", "xy plane 3", 130, -65, 65, 130, -65, 65);
66 }
67 
68 // -------------------------------------------------------------------------
69 std::map< int, bool > PndTrkGemCombinatorial::CombinatorialSuppression()
70 {
71  std::map< int, bool > hitTousability;
72 
73  if(fMCEval == kTRUE && fGemPointArray == NULL) {
74  cerr << "PndTrkGemCombinatorial::CombinatorialSuppression youwant evaluate performances but you did not set the mc TCA!" << endl;
75  return hitTousability;
76  }
77 
78  if(fVerbose) cout << "---------------------- " << endl;
79  double distlimit = 999;
80  // fill maps ===================================
81  // hitid --> position (0...5)
82  // position --> hitid
83  std::map< int, int > hitidTopos;
84  std::multimap< int, int > posTohitid;
85  std::multimap< int, int > posTohitid1;
86  std::multimap< int, int > posTohitid2;
87 
88  if(fDisplayOn) {
89  DrawGeometry();
90  if(fMCEval) DrawMCPoints();
91  }
92 
93  nhits += fGemHitArray->GetEntriesFast();
94  if(fMCEval) npoints += fGemPointArray->GetEntriesFast();
95 
96  for(int ihit = 0; ihit < fGemHitArray->GetEntriesFast(); ihit++) {
97  PndGemHit *hit = (PndGemHit*) fGemHitArray->At(ihit);
98 
99  if(fMCEval) {
100  if(hit->GetRefIndex() == -1) nunasso++;
101  else nasso++;
102  }
103 
104  // statid = 1, 2, 3
105  // sensid = 1, 2
106  // posid = [sensid + 2 * (statid - 1)] - 1
107  // --> position 0, 1, 2, 3, 4, 5
108 
109  int statid = hit->GetStationNr();
110  int sensid = hit->GetSensorNr();
111  int posid = (sensid + 2 * (statid - 1)) - 1;
112 
113  // if(fVerbose) cout << ihit << " hit is on position " << posid << endl;
114  hitidTopos[ihit] = posid;
115  posTohitid.insert(std::pair<int, int>(posid, ihit));
116 
117  // two separate lists for first and second sensor on each station
118  if(sensid == 1) posTohitid1.insert(std::pair<int, int>(posid, ihit));
119  else posTohitid2.insert(std::pair<int, int>(posid, ihit));
120 
121 
122  if(fDisplayOn) {
123  TMarker *mrk = new TMarker(hit->GetPosition().X(), hit->GetPosition().Y(), 4);
124  display->cd(statid);
125  mrk->Draw("SAME");
126 
127  display->cd(4);
128  mrk->Draw("SAME");
129  display->Update();
130  display->Modified();
131 
132 
133  }
134 
135  }
136 
137  // std::multimap< int, int >::iterator it;
138  // for (it = posTohitid.begin(); it != posTohitid.end(); ++it) {
139  // if(fVerbose) cout << "position " << (*it).first << " has hit => " << (*it).second << endl;
140  // }
141  // ==============================================
142 
143 
144  // loop on "first" planes on each station 0, 2, 4
145  // and for each of them loop on the corresponding +1
146  std::multimap< int, int >::iterator it1;
147  std::multimap< int, int >::iterator it2;
148 
149  std::multimap< int, std::pair< int, double > > hit1Topairhit2_distance;
150  std::multimap< int, std::pair< int, double > > hit2Topairhit1_distance;
151  std::vector< int > alone1;
152  std::vector< int > alone2;
153  PndGemHit *hit1 = NULL, *hit2 = NULL;
154  for (it1 = posTohitid1.begin(); it1 != posTohitid1.end(); ++it1) {
155  int posid1 = (*it1).first;
156  int posid2 = posid1 + 1;
157  int hitid1 = (*it1).second;
158  hit1 = (PndGemHit*) fGemHitArray->At(hitid1);
159  TVector2 pos1 = hit1->GetPosition().XYvector();
160 
161  std::pair < std::multimap < int, int >::iterator, std::multimap< int, int >::iterator > ret;
162  ret = posTohitid2.equal_range(posid2);
163 
164  for (it2 = ret.first; it2 != ret.second; ++it2) {
165  // if(fVerbose) cout << "--- accoppiamenti1 " << (*it1).first << " " << (*it1).second << endl;
166  // if(fVerbose) cout << "+++ accoppiamenti2 " << (*it2).first << " " << (*it2).second << endl;
167  int hitid2 = (*it2).second;
168  if(hitid1 == hitid2) continue;
169  hit2 = (PndGemHit*) fGemHitArray->At(hitid2);
170  TVector2 pos2 = hit2->GetPosition().XYvector();
171 
172  // now we have hit1 and hit2 and their distance in xy
173  double distancexy = (pos2 - pos1).Mod();
174  std::pair< int, double > hit2_distance(hitid2, distancexy);
175  hit1Topairhit2_distance.insert(std::pair< int, std::pair< int, double > > (hitid1, hit2_distance));
176  std::pair< int, double > hit1_distance(hitid1, distancexy);
177  hit2Topairhit1_distance.insert(std::pair< int, std::pair< int, double > > (hitid2, hit1_distance));
178  }
179 
180  if(hit1Topairhit2_distance.count(hitid1) == 0) alone1.push_back(hitid1);
181  if(fVerbose) cout << "---> pos 1 hit " << hitid1 << " has counts " << hit1Topairhit2_distance.count(hitid1) << endl;
182  }
183 
184  for(size_t ihit = 0; ihit < alone1.size(); ihit++) {
185  hit1 = (PndGemHit*) fGemHitArray->At(alone1[ihit]);
186  if(fVerbose) cout << "**********alone1 " << alone1[ihit] << " " << hit1->GetRefIndex() << endl;
187  }
188 
189  for (it2 = posTohitid2.begin(); it2 != posTohitid2.end(); ++it2) {
190  int hitid2 = (*it2).second;
191  if(hit2Topairhit1_distance.count(hitid2) == 0) alone2.push_back(hitid2);
192  if(fVerbose) cout << "---> pos 2 hit " << hitid2 << " has counts " << hit2Topairhit1_distance.count(hitid2) << endl;
193  }
194 
195  for(size_t ihit = 0; ihit < alone2.size(); ihit++) {
196  hit2 = (PndGemHit*) fGemHitArray->At(alone2[ihit]);
197  if(fVerbose) cout << "**********alone2 " << alone2[ihit] << " " << hit2->GetRefIndex() << endl;
198  }
199 
200 
201  // create a chosenmap of chosen couples hit1 - hit2
202  std::multimap< int, std::pair< int, double > >::iterator iter;
203  int tmphitid1 = -1;
204  int tmphitid2 = -1;
205 
206  int tmpdigi1[2] = {-1, -1};
207  int tmpdigi2[2] = {-1, -1};
208 
209  double tmpdistance = 1000;
210  std::map< int, std::pair< int, double > > chosenmap;
211  std::multimap< int, int > digiTohitid1;
212  std::multimap< int, int > digiTohitid2;
213  for (iter = hit1Topairhit2_distance.begin(); iter != hit1Topairhit2_distance.end(); ++iter) {
214 
215  int hitid1 = (*iter).first;
216  int hitid2 = ((*iter).second).first;
217  double distance = ((*iter).second).second;
218  if(fVerbose) cout << "hit1: " << hitid1 << " associated to hit2: " << hitid2 << " with dist " << distance << endl;
219 
220  if(tmphitid1 == -1) tmphitid1 = hitid1;
221 
222  if(hitid1 == tmphitid1) { // comparison
223  if(distance < tmpdistance) {
224  tmpdistance = distance;
225  tmphitid2 = hitid2;
226  hit1 = (PndGemHit*) fGemHitArray->At(tmphitid1);
227  hit2 = (PndGemHit*) fGemHitArray->At(tmphitid2);
228  tmpdigi1[0] = hit1->GetDigiNr(0);
229  tmpdigi1[1] = hit1->GetDigiNr(1);
230  tmpdigi2[0] = hit2->GetDigiNr(0);
231  tmpdigi2[1] = hit2->GetDigiNr(1);
232  }
233  }
234  else { // new one
235  if(tmpdistance < distlimit) { // 1 cm is the limit
236  if(fVerbose) cout << "chosen " << tmphitid1 << " " << tmphitid2 << " " << tmpdistance << endl;
237  chosenmap[tmphitid1] = std::pair< int, double >(tmphitid2, tmpdistance);
238  chosenmap[tmphitid2] = std::pair< int, double >(tmphitid1, tmpdistance);
239 // if(fVerbose) cout << "digi1 " << tmpdigi1[0] << " " << tmpdigi1[1] << endl;
240 // if(fVerbose) cout << "digi2 " << tmpdigi2[0] << " " << tmpdigi2[1] << endl;
241 
242  digiTohitid1.insert(std::pair< int, int >(tmpdigi1[0], tmphitid1));
243  digiTohitid1.insert(std::pair< int, int >(tmpdigi1[1], tmphitid1));
244  digiTohitid2.insert(std::pair< int, int >(tmpdigi2[0], tmphitid2));
245  digiTohitid2.insert(std::pair< int, int >(tmpdigi2[1], tmphitid2));
246  }
247  tmpdistance = distance;
248  tmphitid2 = hitid2;
249  tmphitid1 = hitid1;
250  hit1 = (PndGemHit*) fGemHitArray->At(tmphitid1);
251  hit2 = (PndGemHit*) fGemHitArray->At(tmphitid2);
252  tmpdigi1[0] = hit1->GetDigiNr(0);
253  tmpdigi1[1] = hit1->GetDigiNr(1);
254  tmpdigi2[0] = hit2->GetDigiNr(0);
255  tmpdigi2[1] = hit2->GetDigiNr(1);
256  }
257  }
258 
259 
260 
261  // last one
262  if(tmpdistance < distlimit) { // 1 cm is the limit
263  if(fVerbose) cout << "chosen " << tmphitid1 << " " << tmphitid2 << " " << tmpdistance << endl;
264  chosenmap[tmphitid1] = std::pair< int, double >(tmphitid2, tmpdistance);
265  chosenmap[tmphitid2] = std::pair< int, double >(tmphitid1, tmpdistance);
266  // if(fVerbose) cout << "digi1 " << tmpdigi1[0] << " " << tmpdigi1[1] << endl;
267 // if(fVerbose) cout << "digi2 " << tmpdigi2[0] << " " << tmpdigi2[1] << endl;
268 
269  digiTohitid1.insert(std::pair< int, int >(tmpdigi1[0], tmphitid1));
270  digiTohitid1.insert(std::pair< int, int >(tmpdigi1[1], tmphitid1));
271  digiTohitid2.insert(std::pair< int, int >(tmpdigi2[0], tmphitid2));
272  digiTohitid2.insert(std::pair< int, int >(tmpdigi2[1], tmphitid2));
273  }
274 
275  std::map< int, int >::iterator ditr;
276  int tmpdigiid = -1;
277 
278  std::map< int, std::pair< int, double > > chosenmap2;
279  for (ditr = digiTohitid1.begin(); ditr != digiTohitid1.end(); ++ditr) {
280  int digiid = (*ditr).first;
281  if(digiid == tmpdigiid) continue;
282 
283  tmpdigiid = digiid;
284  int ntimes = digiTohitid1.count(digiid);
285  if(ntimes > 1) {
286  std::pair < std::multimap< int, int >::iterator, std::multimap< int, int >::iterator > ret2;
287  ret2 = digiTohitid1.equal_range(digiid);
288  std::multimap< int, int >::iterator itr3;
289 
290  int tmphit1 = -1, tmphit2 = -1 ;
291  double tmpdist = 1000;
292  for (itr3 = ret2.first; itr3 != ret2.second; ++itr3) {
293  int hitid1 = (*itr3).second;
294  int hitid2 = chosenmap[hitid1].first;
295  double distance = chosenmap[hitid1].second;
296  if(distance < tmpdist) {
297  tmphit1 = hitid1;
298  tmphit2 = hitid2;
299  tmpdist = distance;
300  }
301 // if(fVerbose) cout << "DIGI " << digiid << " ntimes " << ntimes << " hit1 " << hitid1 << " hit2 " << hitid2 << " distance " << distance << endl;
302  }
303  chosenmap2[tmphit1] = std::pair< int, double >(tmphit2, tmpdist);
304  chosenmap2[tmphit2] = std::pair< int, double >(tmphit1, tmpdist);
305  }
306  else {
307  int hitid1 = (*ditr).second;
308  int hitid2 = chosenmap[hitid1].first;
309  double distance = chosenmap[hitid1].second;
310  chosenmap2[hitid1] = std::pair< int, double >(hitid2, distance);
311  chosenmap2[hitid2] = std::pair< int, double >(hitid1, distance);
312 
313  }
314  }
315 
316 
317  std::map< int, std::pair < int, double > >::iterator citr;
318  for (citr = chosenmap2.begin(); citr != chosenmap2.end(); ++citr) {
319  if(fVerbose) cout << (*citr).first << " finally associated to " << ((*citr).second).first << " with dist " << ((*citr).second).second << endl;
320 
321  if(fMCEval) {
322  PndGemHit *hitA = (PndGemHit*) fGemHitArray->At((*citr).first);
323  // if(fVerbose) cout << "digi " << hit->GetDigiNr(0) << " " << hit->GetDigiNr(1) << endl;
324 
325  if(hitA->GetRefIndex() == -1) nassowrong++;
326  else nassocorrect++;
327 
328 
329 
330  }
331  }
332 
333 
334 
335 
336  std::vector< int > dontuse;
337  for(int ihit = 0; ihit < fGemHitArray->GetEntriesFast(); ihit++) {
338  hitTousability[ihit] = true;
339  if(chosenmap2.count(ihit) > 0) continue;
340  dontuse.push_back(ihit);
341  hitTousability[ihit] = false;
342  }
343 
344 
345  for(size_t ihit = 0; ihit < dontuse.size(); ihit++) {
346  PndGemHit *hit = (PndGemHit*) fGemHitArray->At(dontuse[ihit]);
347  if(fVerbose) cout << "dontuse " << dontuse[ihit] << " " << hit->GetRefIndex() << endl;
348 
349  if(fMCEval) {
350  if(hit->GetRefIndex() == -1) ndontusecorrect++;
351  else ndontusewrong++;
352  }
353 
354  if(fDisplayOn) {
355  TMarker *mrk = new TMarker(hit->GetPosition().X(), hit->GetPosition().Y(), 4);
356  display->cd(hit->GetStationNr());
357  mrk->SetMarkerColor(kRed);
358  mrk->Draw("SAME");
359 
360  display->cd(4);
361  mrk->Draw("SAME");
362  display->Update();
363  display->Modified();
364  }
365 
366 
367 
368  }
369 
370  if(fDisplayOn) {
371  char goOnChar;
372  cin >> goOnChar;
373  }
374  if(fMCEval) {
375  if(fVerbose > 2) {
376  cout << "TOTAL HITS " << nhits << " TOTAL MC POINTS " << npoints << endl;
377  cout << "ASSO " << nasso << " UN-ASSO " << nunasso << endl;
378  cout << "dontuse WRONG " << ndontusewrong << " RIGHT " << ndontusecorrect << endl;
379  cout << "asso WRONG " << nassowrong << " RIGHT " << nassocorrect << endl;
380  cout << endl;
381  }
382 
383  // if(fVerbose)
384  {
385  cout << "TOTAL number of hits = " << nhits << " of which " << 100. * nasso/nhits << "% TRUE, " << 100. * nunasso/nhits << "% FAKE" << endl;
386  cout << endl;
387  cout << "correctly assiged over the ones which should be " << 100. * nassocorrect/nasso << "%, of the totally assigned " << 100. * nassocorrect/(nassowrong + nassocorrect) << "%" << endl;
388  cout << "wrongly assiged over the totally assigned " << 100. * nassowrong/(nassowrong + nassocorrect) << "%" << endl;
389  cout << endl;
390  cout << "correctly NOT assiged over the ones which should not " << 100. * ndontusecorrect/nunasso << "%, of the totally un-assigned " << 100. * ndontusecorrect/(ndontusecorrect + ndontusewrong) << "%" << endl;
391  cout << "wrongly NOT assiged over the totally un-assigned " << 100. * ndontusewrong/(ndontusecorrect + ndontusewrong) << "%" << endl;
392  }
393  }
394 
395 
396  return hitTousability;
397 }
398 
399 void PndTrkGemCombinatorial::DrawGeometry() {
400 
401  if(hxy == NULL) hxy = new TH2F("hxy", "xy plane", 110, -55, 55, 110, -55, 55);
402  else hxy->Reset();
403  if(hxy1 == NULL) hxy1 = new TH2F("hxy1", "xy plane 1", 110, -55, 55, 110, -55, 55);
404  else hxy1->Reset();
405  if(hxy2 == NULL) hxy2 = new TH2F("hxy2", "xy plane 2", 110, -55, 55, 110, -55, 55);
406  else hxy2->Reset();
407  if(hxy3 == NULL) hxy3 = new TH2F("hxy3", "xy plane 3", 110, -55, 55, 110, -55, 55);
408  else hxy3->Reset();
409  display->cd(1);
410  hxy1->SetStats(kFALSE);
411  hxy1->Draw();
412  display->cd(2);
413  hxy2->SetStats(kFALSE);
414  hxy2->Draw();
415  display->cd(3);
416  hxy3->SetStats(kFALSE);
417  hxy3->Draw();
418  display->cd(4);
419  hxy->SetStats(kFALSE);
420  hxy->Draw();
421 
422  display->Update();
423  display->Modified();
424 
425 }
426 
427 void PndTrkGemCombinatorial::DrawMCPoints() {
428  for(int ipnt = 0; ipnt < fGemPointArray->GetEntriesFast(); ipnt++) {
429  PndGemMCPoint *point = (PndGemMCPoint*) fGemPointArray->At(ipnt);
430  int sensid = point->GetSensorId();
431  int statid = 0;
432 
433  // CHECK these are not fine
434  if(sensid == 273 || sensid == 289) statid = 1;
435  else if(sensid == 529 || sensid == 545) statid = 2;
436  if(sensid == 785 || sensid == 801) statid = 3;
437 
438  TVector3 pos;
439  point->Position(pos);
440  TMarker *mrk = new TMarker(pos.X(), pos.Y(), 3);
441  display->cd(statid);
442  mrk->SetMarkerColor(kBlue);
443  mrk->Draw("SAME");
444 
445  display->cd(4);
446  mrk->Draw("SAME");
447  display->Update();
448  display->Modified();
449  }
450 }
451 
452 
453 
454 ClassImp(PndTrkGemCombinatorial)
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
fVerbose
int fVerbose
Definition: poormantracks.C:24
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
PndGemHit::GetSensorNr
Int_t GetSensorNr() const
Definition: PndGemHit.h:83
PndTrkGemCombinatorial::CombinatorialSuppression
std::map< int, bool > CombinatorialSuppression()
Definition: PndTrkGemCombinatorial.cxx:69
ndontusecorrect
int ndontusecorrect
Definition: PndTrkGemCombinatorial.cxx:31
PndGemMCPoint::GetSensorId
Int_t GetSensorId() const
Definition: PndGemMCPoint.h:90
nassocorrect
int nassocorrect
Definition: PndTrkGemCombinatorial.cxx:31
point
TClonesArray * point
Definition: anaLmdDigi.C:29
PndGemHit::GetStationNr
Int_t GetStationNr() const
Definition: PndGemHit.h:81
ndontusewrong
int ndontusewrong
Definition: PndTrkGemCombinatorial.cxx:31
hit
int hit(Int_t nEvents=0, TString inFile="sim.root", TString parFile="par.root", TString inDigi="digi.root", TString outFile="hit.root", Int_t timeBased=0)
Definition: hit.C:1
nassowrong
int nassowrong
Definition: PndTrkGemCombinatorial.cxx:31
PndGemHit::GetPosition
TVector3 GetPosition() const
Definition: PndGemHit.h:71
ClassImp
ClassImp(PndAnaContFact)
PndGemHit::GetDigiNr
Int_t GetDigiNr(Int_t iside) const
Definition: PndGemHit.h:77
Generated on Sun Aug 18 2019 03:01:52 for FairRoot/PandaRoot by   doxygen 1.8.8