Class Reference for E1039 Core & Analysis Software
PatternDBUtil.cxx
Go to the documentation of this file.
1 
9 #include "PatternDBUtil.h"
10 
11 #include <TTree.h>
12 #include <TFile.h>
13 
14 #define LogInfo(message) std::cout << "DEBUG: " << __FILE__ << " " << __LINE__ << " " << __FUNCTION__ << " ::: " << message << std::endl
15 
16 #define _DEBUG_
17 
18 #define _D1_1_6_
19 
20 //#define _RESOLUTION1_ 2
21 //#define _RESOLUTION2_ 2
22 //#define _RESOLUTION3_ 2
23 
24 int PatternDBUtil::verbosity = 0;
25 bool PatternDBUtil::_loose_mode = false;
26 int PatternDBUtil::_RESOLUTION1_ = 2;
27 int PatternDBUtil::_RESOLUTION2_ = 2;
28 int PatternDBUtil::_RESOLUTION3_ = 2;
29 
30 std::map<unsigned int, unsigned int> PatternDBUtil::_detid_view = {
31  {3, 0},
32  {4, 1},
33  {1, 2},
34  {2, 3},
35  {5, 4},
36  {6, 5},
37 
38  {9, 0},
39  {10, 1},
40  {11, 2},
41  {12, 3},
42  {7, 4},
43  {8, 5},
44 
45  {16, 0},
46  {15, 1},
47  {17, 2},
48  {18, 3},
49  {13, 4},
50  {14, 5},
51 
52  {22, 0},
53  {21, 1},
54  {23, 2},
55  {24, 3},
56  {20, 4},
57  {19, 5},
58 
59  {28, 0},
60  {27, 1},
61  {30, 2},
62  {29, 3},
63  {26, 4},
64  {25, 5}
65 };
66 
67 int PatternDBUtil::BuildPatternDB(const std::string &fin, const std::string & fout, PatternDB& db) {
68 
69 #ifdef _DEBUG_
70  TFile *feval = TFile::Open("PatternDBUtilEval.root","recreate");
71  TTree *Teval = new TTree("T","PatternDBUtilEval");
72  int ntrack = 0;
73  float acc_prob_1 = 0;
74  float acc_prob_2 = 0;
75  float acc_prob_3 = 0;
76  float acc_prob_23 = 0;
77  float acc_prob_123 = 0;
78 
79  Teval->Branch("ntrack", &ntrack, "ntrack/I");
80  Teval->Branch("acc_prob_1", &acc_prob_1, "acc_prob_1/F");
81  Teval->Branch("acc_prob_2", &acc_prob_2, "acc_prob_2/F");
82  Teval->Branch("acc_prob_3", &acc_prob_3, "acc_prob_3/F");
83  Teval->Branch("acc_prob_23", &acc_prob_23, "acc_prob_23/F");
84  Teval->Branch("acc_prob_123", &acc_prob_123, "acc_prob_123/F");
85 
86  int nacc_1 = 0;
87  int nacc_2 = 0;
88  int nacc_3 = 0;
89  int nacc_23 = 0;
90  int nacc_123 = 0;
91  const int interval = 1000;
92 #endif
93 
94  if(verbosity >= 2) {
95  LogInfo("PatternDBUtil::BuildPatternDB from " << fin);
96  }
97 
98  TFile *f_in = TFile::Open(fin.c_str(), "read");
99  if(!f_in) {
100  LogInfo(fin << "not found!");
101  return -1;
102  }
103 
104  TTree *T = (TTree*) f_in->Get("T");
105  if(!T) {
106  LogInfo("TTree T not found in " << fin);
107  return -1;
108  }
109 
110  int n_particles = 0;
111  int gndc[1000];
112  int elmid [1000][55];
113 
114  T->SetBranchAddress("n_tracks", &n_particles);
115  T->SetBranchAddress("gelmid", &elmid);
116  T->SetBranchAddress("gndc", &gndc);
117 
118  //PatternDB db;
119 
120 #ifdef _DEBUG_
121  ntrack = 0;
122  nacc_1 = 0;
123  nacc_2 = 0;
124  nacc_3 = 0;
125  nacc_23 = 0;
126  nacc_123 = 0;
127 #endif
128 
129 // TrackletKey key = EncodeTrackletKey(PatternDB::DC1, 94, 0, 116, 0, 112, 0);
130 // db.St1.insert(key);
131 
132  //for(int ientry=0;ientry<0;++ientry) {
133  for(int ientry=0;ientry<T->GetEntries();++ientry) {
134  T->GetEntry(ientry);
135 
136  for(int ipar=0; ipar<n_particles; ++ipar) {
137 
138  if(!(gndc[ipar]>17)) continue;
139 
140 #ifdef _D1_1_6_
141  unsigned int D1U = elmid[ipar][1]; // 1
142  unsigned int D1Up = elmid[ipar][2]; // 2
143  unsigned int D1X = elmid[ipar][3]; // 3
144  unsigned int D1Xp = elmid[ipar][4]; // 4
145  unsigned int D1V = elmid[ipar][5]; // 5
146  unsigned int D1Vp = elmid[ipar][6]; // 6
147 #else
148  unsigned int D1V = elmid[ipar][7];
149  unsigned int D1Vp = elmid[ipar][8];
150  unsigned int D1X = elmid[ipar][9];
151  unsigned int D1Xp = elmid[ipar][10];
152  unsigned int D1U = elmid[ipar][11];
153  unsigned int D1Up = elmid[ipar][12];
154 #endif
155  unsigned int D2V = elmid[ipar][13];
156  unsigned int D2Vp = elmid[ipar][14];
157  unsigned int D2Xp = elmid[ipar][15];
158  unsigned int D2X = elmid[ipar][16];
159  unsigned int D2U = elmid[ipar][17];
160  unsigned int D2Up = elmid[ipar][18];
161 
162  unsigned int D3pVp = elmid[ipar][19];
163  unsigned int D3pV = elmid[ipar][20];
164  unsigned int D3pXp = elmid[ipar][21];
165  unsigned int D3pX = elmid[ipar][22];
166  unsigned int D3pUp = elmid[ipar][23];
167  unsigned int D3pU = elmid[ipar][24];
168 
169  unsigned int D3mVp = elmid[ipar][25];
170  unsigned int D3mV = elmid[ipar][26];
171  unsigned int D3mXp = elmid[ipar][27];
172  unsigned int D3mX = elmid[ipar][28];
173  unsigned int D3mUp = elmid[ipar][29];
174  unsigned int D3mU = elmid[ipar][30];
175 
176 
177 #ifdef _DEBUG_
178  auto size_1 = db.St1.size();
179  auto size_2 = db.St2.size();
180  auto size_3 = db.St3.size();
181  auto size_23 = db.St23.size();
182  auto size_123 = db.St123.size();
183 #endif
184 
185  // Single key
186  TrackletKey key1 = EncodeTrackletKey(PatternDB::DC1, D1X, D1Xp, D1U, D1Up, D1V, D1Vp);
187  TrackletKey key2 = EncodeTrackletKey(PatternDB::DC2, D2X, D2Xp, D2U, D2Up, D2V, D2Vp);
188  TrackletKey key3p = EncodeTrackletKey(PatternDB::DC3p, D3pX, D3pXp, D3pU, D3pUp, D3pV, D3pVp);
189  TrackletKey key3m = EncodeTrackletKey(PatternDB::DC3m, D3mX, D3mXp, D3mU, D3mUp, D3mV, D3mVp);
190 
191  if(key1 != PatternDB::ERR_KEY) db.St1.insert(key1);
192  if(key2 != PatternDB::ERR_KEY) db.St2.insert(key2);
193  if(key3p != PatternDB::ERR_KEY) db.St3.insert(key3p);
194  if(key3m != PatternDB::ERR_KEY) db.St3.insert(key3m);
195 
196  if(key2 != PatternDB::ERR_KEY and key3p != PatternDB::ERR_KEY) {
197  db.St23.insert(PartTrackKey(key2,key3p));
198  if(key1 != PatternDB::ERR_KEY) {
199  db.St123.insert(GlobTrackKey(key1,key2,key3p));
200  }
201  }
202 
203  if(key2 != PatternDB::ERR_KEY and key3m != PatternDB::ERR_KEY) {
204  db.St23.insert(PartTrackKey(key2,key3m));
205  if(key1 != PatternDB::ERR_KEY) {
206  db.St123.insert(GlobTrackKey(key1,key2,key3m));
207  }
208  }
209 
210 #ifdef _DEBUG_
211 
212  if(db.St1.size()>size_1) ++nacc_1;
213  if(db.St2.size()>size_2) ++nacc_2;
214  if(db.St3.size()>size_3) ++nacc_3;
215  if(db.St23.size()>size_23) ++nacc_23;
216  if(db.St123.size()>size_123) ++nacc_123;
217 
218  if(ntrack>0 and ntrack%interval==0) {
219  acc_prob_1 = 1.*nacc_1/interval;
220  acc_prob_2 = 1.*nacc_2/interval;
221  acc_prob_3 = 1.*nacc_3/interval;
222  acc_prob_23 = 1.*nacc_23/interval;
223  acc_prob_123 = 1.*nacc_123/interval;
224 
225  Teval->Fill();
226 
227  nacc_1 = 0;
228  nacc_2 = 0;
229  nacc_3 = 0;
230  nacc_23 = 0;
231  nacc_123 = 0;
232 
233  //LogInfo("ntrack: "<< ntrack);
234  //LogInfo("acc_prob_1: "<< acc_prob_1);
235  //LogInfo("acc_prob_2: "<< acc_prob_2);
236  //LogInfo("acc_prob_3: "<< acc_prob_3);
237  //LogInfo("acc_prob_23: "<< acc_prob_23);
238  //LogInfo("acc_prob_123: "<< acc_prob_123);
239  }
240  ++ntrack;
241 #endif
242  }
243  }
244 
245  if(verbosity >= 2) {
246  LogInfo("PatternDBUtil::BuildPatternDB from " << fin);
247  LogInfo("db_st23.size(): " << db.St23.size());
248  db.print();
249  }
250 
251  //TODO remove this debug code
252 // if(true) {
253 // LogInfo("");
254 // unsigned int e[3] = {78, 104, 93};
255 // //unsigned int e[3] = {77, 102, 91};
256 // TrackletKey tmp = PatternDBUtil::EncodeTrackletKey(PatternDB::DC1, e[0], 0, e[1], 0, e[2], 0);
257 // std::cout << tmp;
258 //
259 // if(db.St1.find(tmp)!=db.St1.end())
260 // std::cout << "Found!" << std::endl;
261 // else
262 // std::cout << "Not Found!" << std::endl;
263 //
264 // auto size = db.St1.size();
265 // db.St1.insert(tmp);
266 // if(db.St1.size()>size)
267 // std::cout << "Inserted!" << std::endl;
268 // else
269 // std::cout << "Not Inserted!" << std::endl;
270 //
271 // if(db.St1.find(tmp)!=db.St1.end())
272 // std::cout << "Found!" << std::endl;
273 // else
274 // std::cout << "Not Found!" << std::endl;
275 // }
276 
277 
278 
279  TFile *f_out = TFile::Open(fout.c_str(), "recreate");
280 
281  TTree *T_out_St1 = new TTree("St1","St1");
282  TrackletKey *b_St1 = new TrackletKey();
283  T_out_St1->Branch("key",&b_St1);
284  for(auto key : db.St1) {
285  b_St1 = &key;
286  T_out_St1->Fill();
287  }
288 
289  TTree *T_out_St2 = new TTree("St2","St2");
290  TrackletKey *b_St2 = new TrackletKey();
291  T_out_St2->Branch("key",&b_St2);
292  for(auto key : db.St2) {
293  b_St2 = &key;
294  T_out_St2->Fill();
295  }
296 
297  TTree *T_out_St3 = new TTree("St3","St3");
298  TrackletKey *b_St3 = new TrackletKey();
299  T_out_St3->Branch("key",&b_St3);
300  for(auto key : db.St3) {
301  b_St3 = &key;
302  T_out_St3->Fill();
303  }
304 
305  TTree *T_out_St23 = new TTree("St23","St23");
306  PartTrackKey *b_St23 = new PartTrackKey();
307  T_out_St23->Branch("key",&b_St23);
308  for(auto key : db.St23) {
309  b_St23 = &key;
310  T_out_St23->Fill();
311  }
312 
313  TTree *T_out_St123 = new TTree("St123","St123");
314  GlobTrackKey *b_St123 = new GlobTrackKey();
315  T_out_St123->Branch("key",&b_St123);
316  for(auto key : db.St123) {
317  b_St123 = &key;
318  T_out_St123->Fill();
319  }
320 
321  f_out->cd();
322  T_out_St1->Write();
323  T_out_St2->Write();
324  T_out_St3->Write();
325  T_out_St23->Write();
326  T_out_St123->Write();
327  f_out->Close();
328 
329 #ifdef _DEBUG_
330  feval->cd();
331  Teval->Write();
332  feval->Close();
333 #endif
334  return 0;
335 }
336 
337 PatternDB* PatternDBUtil::LoadPatternDB(const std::string& fin) {
338  if(verbosity >= 2) {
339  LogInfo("PatternDBUtil::BuildPatternDB from " << fin);
340  }
341 
342  TFile *f_in = TFile::Open(fin.c_str(), "read");
343  if(!f_in) {
344  LogInfo(fin << "not found!");
345  return nullptr;
346  }
347 
348  PatternDB* db = new PatternDB();
349 
350  TTree *St1 = (TTree*) f_in->Get("St1");
351  if(St1) {
352  TrackletKey *key = new TrackletKey();
353  St1->SetBranchAddress("key", &key);
354  for(int ientry=0;ientry<St1->GetEntries();++ientry) {
355  St1->GetEntry(ientry);
356  db->St1.insert(db->St1.end(), *key);
357  }
358  delete key;
359  }
360 
361  TTree *St2 = (TTree*) f_in->Get("St2");
362  if(St2) {
363  TrackletKey *key = new TrackletKey();
364  St2->SetBranchAddress("key", &key);
365  for(int ientry=0;ientry<St2->GetEntries();++ientry) {
366  St2->GetEntry(ientry);
367  db->St2.insert(db->St2.end(), *key);
368  }
369  delete key;
370  }
371 
372  TTree *St3 = (TTree*) f_in->Get("St3");
373  if(St3) {
374  TrackletKey *key = new TrackletKey();
375  St3->SetBranchAddress("key", &key);
376  for(int ientry=0;ientry<St3->GetEntries();++ientry) {
377  St3->GetEntry(ientry);
378  db->St3.insert(db->St3.end(), *key);
379  }
380  delete key;
381  }
382 
383  TTree *St23 = (TTree*) f_in->Get("St23");
384  if(St23) {
385  PartTrackKey *key = new PartTrackKey();
386  St23->SetBranchAddress("key", &key);
387  for(int ientry=0;ientry<St23->GetEntries();++ientry) {
388  St23->GetEntry(ientry);
389  db->St23.insert(db->St23.end(), PartTrackKey(*key));
390  }
391  delete key;
392  }
393 
394  TTree *St123 = (TTree*) f_in->Get("St123");
395  if(St123) {
396  GlobTrackKey *key = new GlobTrackKey();
397  St123->SetBranchAddress("key", &key);
398  for(int ientry=0;ientry<St123->GetEntries();++ientry) {
399  St123->GetEntry(ientry);
400  db->St123.insert(db->St123.end(), *key);
401  }
402  delete key;
403  }
404 
405  if(verbosity >= 2) {
406  LogInfo("PatternDBUtil::LoadPatternDB from " << fin);
407  if(db) db->print();
408  else LogInfo("PatternDB NOT found!!");
409  }
410 
411  return db;
412 }
413 
414 TrackletKey PatternDBUtil::EncodeTrackletKey(
415  PatternDB::STATION ST,
416  const unsigned int X, const unsigned int Xp,
417  const unsigned int U, const unsigned int Up,
418  const unsigned int V, const unsigned int Vp) {
419 
420  // TODO add range check
421  if(
422  ST > PatternDB::DC3m or ST < PatternDB::DC1 or
423  X > 255 or Xp > 255 or
424  U > 255 or Up > 255 or
425  V > 255 or Vp > 255
426  ) {
427  return PatternDB::ERR_KEY;
428  }
429 
430  int x = X;
431  int u = U;
432  int v = V;
433 
434  if(_loose_mode) {
435  x = X > 0 ? X : Xp;
436  u = U > 0 ? U : Up;
437  v = V > 0 ? V : Vp;
438  }
439 
440  if(x == 0 and u == 0 and v == 0) {
441  return PatternDB::ERR_KEY;
442  }
443 
444  if(ST == PatternDB::DC1)
445  return TrackletKey(ST, x/_RESOLUTION1_, u/_RESOLUTION1_, v/_RESOLUTION1_);
446  else if(ST == PatternDB::DC2)
447  return TrackletKey(ST, x/_RESOLUTION2_, u/_RESOLUTION2_, v/_RESOLUTION2_);
448  else
449  return TrackletKey(ST, x/_RESOLUTION3_, u/_RESOLUTION3_, v/_RESOLUTION3_);
450 }
451 
452 TrackletKey PatternDBUtil::GetTrackletKey(
453  const Tracklet tracklet,
454  const PatternDB::STATION station) {
455 
456  if(station > PatternDB::DC3m or station < PatternDB::DC1) return PatternDB::ERR_KEY;
457 
458  std::vector<unsigned int> elmids;
459  elmids.resize(6);
460 
461  for (auto ptr_hit = tracklet.hits.begin();
462  ptr_hit != tracklet.hits.end(); ++ptr_hit) {
463  auto hit = &ptr_hit->hit;
464  if (hit->index < 0) continue;
465  unsigned int det_id = hit->detectorID;
466  if(station==PatternDB::DC1 and !(det_id>=1 and det_id<=12)) continue;
467  if(station==PatternDB::DC2 and !(det_id>=13 and det_id<=18)) continue;
468  if(station==PatternDB::DC3p and !(det_id>=19 and det_id<=24)) continue;
469  if(station==PatternDB::DC3m and !(det_id>=25 and det_id<=30)) continue;
470 
471  elmids[_detid_view[det_id]] = hit->elementID;
472  }
473 
474  return EncodeTrackletKey(
475  station,
476  elmids[0],
477  elmids[1],
478  elmids[2],
479  elmids[3],
480  elmids[4],
481  elmids[5]
482  );
483 }
484 
485 TrackletKey PatternDBUtil::GetTrackletKey(
486  const std::vector< std::pair<unsigned int, unsigned int> >& det_elem_pairs,
487  const PatternDB::STATION& station) {
488  if(station > PatternDB::DC3m or station < PatternDB::DC1) return PatternDB::ERR_KEY;
489 
490  std::vector<unsigned int> elmids;
491  elmids.resize(6);
492 
493  for (auto pair : det_elem_pairs) {
494  unsigned int det_id = pair.first;
495  if(station==PatternDB::DC1 and !(det_id>=1 and det_id<=12)) continue;
496  if(station==PatternDB::DC2 and !(det_id>=13 and det_id<=18)) continue;
497  if(station==PatternDB::DC3p and !(det_id>=19 and det_id<=24)) continue;
498  if(station==PatternDB::DC3m and !(det_id>=25 and det_id<=30)) continue;
499 
500  elmids[_detid_view[det_id]] = pair.second;
501  }
502 
503  return EncodeTrackletKey(
504  station,
505  elmids[0],
506  elmids[1],
507  elmids[2],
508  elmids[3],
509  elmids[4],
510  elmids[5]
511  );
512 }
513 
514 
515 
516 
517 
LogInfo
#define LogInfo(message)
Definition: PatternDBUtil.cxx:14
PatternDBUtil::EncodeTrackletKey
static TrackletKey EncodeTrackletKey(PatternDB::STATION, const unsigned int X, const unsigned int Xp, const unsigned int U, const unsigned int Up, const unsigned int V, const unsigned int Vp)
Definition: PatternDBUtil.cxx:414
PatternDBUtil::BuildPatternDB
static int BuildPatternDB(const std::string &fin, const std::string &fout, PatternDB &db)
Definition: PatternDBUtil.cxx:67
PatternDBUtil::LoadPatternDB
static PatternDB * LoadPatternDB(const std::string &fin)
Definition: PatternDBUtil.cxx:337
PatternDBUtil::GetTrackletKey
static TrackletKey GetTrackletKey(const Tracklet tracklet, const PatternDB::STATION station)
Definition: PatternDBUtil.cxx:452
PatternDB
PatternDB interface objects.
Definition: PatternDB.h:117
PatternDB::St2
std::set< TrackletKey > St2
Definition: PatternDB.h:143
PatternDB::ERR_KEY
static const TrackletKey ERR_KEY
Definition: PatternDB.h:124
PatternDB::St123
std::set< GlobTrackKey > St123
Definition: PatternDB.h:146
PatternDB::print
void print()
Definition: PatternDB.cxx:157
PatternDB::St23
std::set< PartTrackKey > St23
Definition: PatternDB.h:145
PatternDB::St1
std::set< TrackletKey > St1
Definition: PatternDB.h:142
PatternDB::St3
std::set< TrackletKey > St3
Definition: PatternDB.h:144
Tracklet::hits
std::list< SignedHit > hits
Definition: FastTracklet.h:228