RegPixelMap.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file RegPixelMap.cxx
3 /// \brief RegPixelMap for RegCNN modifed from PixelMap.cxx
4 /// \author Ilsoo Seong - iseong@uci.edu
5 ///
6 /// \modified Wenjie Wu - wenjieww@uci.edu
7 ////////////////////////////////////////////////////////////////////////
8 
9 #include <cassert>
10 #include <iostream>
11 #include <ostream>
12 #include <iomanip>
14 
15 namespace cnn
16 {
17 
18  RegPixelMap::RegPixelMap(unsigned int nWire, unsigned int nWRes,
19  unsigned int nTdc, unsigned int nTRes, const RegCNNBoundary& bound, const bool& prongOnly):
20  fNWire(nWire),
21  fNWRes(nWRes),
22  fNTdc(nTdc),
23  fNTRes(nTRes),
24  fInPM(0),
25  fTPC(0),
26  fdist(100000),
27  fPE(nWire*nTdc),
28  fPEX(nWire*nTdc),
29  fPEY(nWire*nTdc),
30  fPEZ(nWire*nTdc),
31  fPur(nWire*nTdc),
32  fPurX(nWire*nTdc),
33  fPurY(nWire*nTdc),
34  fPurZ(nWire*nTdc),
35  fLab(nWire*nTdc),
36  fLabX(nWire*nTdc),
37  fLabY(nWire*nTdc),
38  fLabZ(nWire*nTdc),
39  fBound(bound),
40  fProngOnly(prongOnly),
41  fProngTagX(nWire*nTdc),
42  fProngTagY(nWire*nTdc),
43  fProngTagZ(nWire*nTdc)
44  {
45  std::cout<<"here :"<<fNWire<<", "<<fNTdc<<", "<<fNWRes<<", "<<fNTRes<<std::endl;
46  }
47 
49  {
50  unsigned int i = 0;
51 
52  for(const auto& pe:fPE){
53  input[i] = pe;
54  ++i;
55  }
56 
57  }
58 
59 
60  void RegPixelMap::Add(const int& wire, const int& tdc, const unsigned int& view, const double& pe, const unsigned int& tpc, int hit_prong_tag)
61  {
62  // keep these for now although we only use fPE
63  const HitType label = kEmptyHit;
64  const double purity=0.0;
65  if(fBound.IsWithin(wire, tdc, view)){
66  fInPM = 1; // any hit within the boundary
67  GetTPC(wire, tdc, view, tpc);
68  fPE[GlobalToIndex(wire, tdc, view)] += (float)pe;
69  fLab[GlobalToIndex(wire,tdc, view)] = label;
70  fPur[GlobalToIndex(wire,tdc, view)] = purity;
71  if(view==0){
72  fPEX[GlobalToIndexSingle(wire,tdc, view)] += (float)pe;
73  fProngTagX[GlobalToIndexSingle(wire, tdc, view)] = hit_prong_tag;
74  fLabX[GlobalToIndexSingle(wire,tdc, view)] = label;
75  fPurX[GlobalToIndexSingle(wire,tdc, view)] = purity;
76  // FIXIT
77  //if (pe!=0) {
78  // std::cout<<view<<" | ";
79  // std::cout<<wire<<", "<<tdc<<", "<<GlobalToIndexSingle(wire,tdc,view)<<", "<<pe<<std::endl;
80  //}
81  }
82  if(view==1){
83  fPEY[GlobalToIndexSingle(wire,tdc, view)] += (float)pe;
84  fProngTagY[GlobalToIndexSingle(wire, tdc, view)] = hit_prong_tag;
85  fLabY[GlobalToIndexSingle(wire,tdc, view)] = label;
86  fPurY[GlobalToIndexSingle(wire,tdc, view)] = purity;
87  }
88  if(view==2){
89  fPEZ[GlobalToIndexSingle(wire,tdc, view)] += (float)pe;
90  fProngTagZ[GlobalToIndexSingle(wire, tdc, view)] = hit_prong_tag;
91  fLabZ[GlobalToIndexSingle(wire,tdc, view)] = label;
92  fPurZ[GlobalToIndexSingle(wire,tdc, view)] = purity;
93  //std::cout << "Q = " << pe << " " << fPEZ[GlobalToIndexSingle(wire,tdc, view)] << std::endl;
94  }
95  }
96  }
97 
99  // fProngOnly=True means only the primary prong is selected to creat pixel maps
100  // and that prong needs to be either a muon or antimuon (FIXIT)?
101  // Caveat 1: the prong tag of that pixel is determined by the track id of the track
102  // associlated with the last hit that associated with that pixel. It means
103  // if this pixel-associated hits belong to different prongs, we could
104  // either add deposited energy from other prongs (if other prongs' hit is
105  // in the front), or throw the energy from primary prong away (if there is
106  // one or more hits from other prongs in the last). A better way could be
107  // determining the prong tag by looping the hits, instead of the pixels.
108  // That means we should create a new pixel map only with the spacepoints
109  // associated with the primary prong, instead of using the prong tag
110  // (little effect on the results, ignored for now)
111  if (fProngOnly) {
112  std::cout<<"Do Prong Only selection ......"<<std::endl;
113  for (unsigned int i_p= 0; i_p< fPE.size(); ++i_p) {
114  if (fProngTagX[i_p] != 0)
115  fPEX[i_p] = 0;
116  if (fProngTagY[i_p] != 0)
117  fPEY[i_p] = 0;
118  if (fProngTagZ[i_p] != 0)
119  fPEZ[i_p] = 0;
120  } // end of i_p
121  } // end of fProngOnly
122  }
123 
124  unsigned int RegPixelMap::GlobalToIndex(const int& wire,
125  const int& tdc,
126  const unsigned int& view)
127  {
128 
129  //int meanWire = (fBound.LastWire(view)+fBound.FirstWire(view))/2;
130  int meanWire = (fBound.LastWire(view)+fBound.FirstWire(view)+(int)fNWRes)/2;
131  int meanTDC = (fBound.LastTDC(view)+fBound.FirstTDC(view)+(int)fNTRes/2)/2;
132 
133  //unsigned int internalWire = (unsigned int)( (wire-meanWire) + fNWire/2 );
134  unsigned int internalWire = (unsigned int)(round((float)((unsigned int)(wire-meanWire)+fNWire*fNWRes/2)/(float)fNWRes) );
135  //unsigned int internalTdc = (unsigned int)( round((float)(tdc-meanTDC)/(float)fNTRes + fNTdc/2) );
136  unsigned int internalTdc = (unsigned int)( round((float)((unsigned int)(tdc-meanTDC)+fNTdc*fNTRes/2)/(float)fNTRes) );
137 
138  unsigned int index = internalWire * fNTdc + internalTdc % fNTdc;
139 
140  //if (internalTdc <1.1){
141  //std::cout << "====> " << meanWire <<" " << meanTDC << " " << internalWire << " " << internalTdc << " " << index << std::endl;
142  //std::cout << " => " << wire << " " << tdc << " " << wire-meanWire << " " << round((float)(tdc-meanTDC)/(float)fNTRes) << std::endl;
143  //}
144  assert(index < fPE.size());
145  return index;
146  }
147 
148  unsigned int RegPixelMap::LocalToIndex(const unsigned int& wire,
149  const unsigned int& tdc) const
150  {
151  unsigned int index = wire * fNTdc + tdc % fNTdc;
152 
153  assert(index < fPE.size());
154  return index;
155  }
156 
157  unsigned int RegPixelMap::GlobalToIndexSingle(const int& wire,
158  const int& tdc,
159  const unsigned int& view)
160 
161  {
162 
163  //int meanWire = (fBound.LastWire(view)+fBound.FirstWire(view))/2;
164  int meanWire = (fBound.LastWire(view)+fBound.FirstWire(view)+(int)fNWRes/2)/2;
165  int meanTDC = (fBound.LastTDC(view)+fBound.FirstTDC(view)+(int)fNTRes/2)/2;
166 
167  //unsigned int internalWire = (unsigned int)( (wire-meanWire) + fNWire/2 );
168  //unsigned int internalWire = (unsigned int)( round((float)((unsigned int)(wire-meanWire)+fNWire*fNWRes/2)/(float)fNWRes) );
169  unsigned int internalWire = (unsigned int)(round((float)(wire - meanWire + fNWire*fNWRes/2)/(float)fNWRes));
170 
171  //unsigned int internalTdc = (unsigned int)( round( (float)(tdc-meanTDC)/(float)fNTRes + fNTdc/2) );
172  //unsigned int internalTdc = (unsigned int)( round((float)((unsigned int)(tdc-meanTDC)+fNTdc*fNTRes/2)/(float)fNTRes) );
173  unsigned int internalTdc = (unsigned int)(round((float)(tdc - meanTDC + fNTdc*fNTRes/2)/(float)fNTRes));
174 
175  unsigned int index = internalWire * fNTdc + internalTdc % fNTdc;
176 
177  assert(index < fPEX.size());
178 
179  return index;
180  }
181 
182  void RegPixelMap::GetTPC(const int& wire,
183  const int& tdc,
184  const unsigned int& view,
185  const unsigned int& tpc)
186 
187  {
188 
189  //int meanWire = (fBound.LastWire(view)+fBound.FirstWire(view))/2;
190  int meanWire = (fBound.LastWire(view)+fBound.FirstWire(view)+(int)fNWRes)/2;
191  int meanTDC = (fBound.LastTDC(view)+fBound.FirstTDC(view)+(int)fNTRes/2)/2;
192  double dist = pow((wire-meanWire)*0.5,2)+pow((tdc-meanTDC)*0.08,2);
193  dist = sqrt(dist);
194  if (dist < fdist){
195  fdist = dist;
196  fTPC = tpc;
197  }
198  }
199 
200 
202  {
203 
204  // Start by doing even wires
205  for(unsigned int iTdc = 0; iTdc < fNTdc; ++iTdc)
206  {
207  for(unsigned int iWire = 0; iWire < fNWire; iWire += 2)
208  {
209  unsigned int index = LocalToIndex(iWire, iTdc);
210  if( fPE[index] > 0)
211  {
212  std::cout << "*";
213  }
214  else
215  {
216  std::cout << " ";
217  }
218 
219  }
220  std::cout << std::endl;
221  }
222  // Then do odd wires
223  for(unsigned int iTdc = 0; iTdc < fNTdc; ++iTdc)
224  {
225  for(unsigned int iWire = 1; iWire < fNWire; iWire += 2)
226  {
227  unsigned int index = LocalToIndex(iWire, iTdc);
228  if( fPE[index] > 0)
229  {
230  std::cout << "*";
231  }
232  else
233  {
234  std::cout << " ";
235  }
236 
237  }
238  std::cout << std::endl;
239  }
240 
241  }
242 
243  TH2F* RegPixelMap::ToTH2() const
244  {
245 
246  // Create a histogram, use twice as many tdcs to distinguish views
247  TH2F* hist = new TH2F("RegPixelMap", ";Wire;Tdc", fNWire, 0, fNWire,
248  fNTdc*3, 0, fNTdc*3);
249 
250  for(unsigned int iWire = 0; iWire < fNWire; ++iWire)
251  {
252  for(unsigned int iTdc = 0; iTdc < fNTdc; ++iTdc)
253  {
254  // Add 1 to in each bin to skip underflow
255  hist->SetBinContent(iWire+1, iTdc + fNTdc*(iWire%3) + 1,
256  fPE[LocalToIndex(iWire, iTdc)]);
257 
258  }
259  }
260  return hist;
261  }
262 
263  TH2F* RegPixelMap::ToLabTH2() const
264  {
265 
266  // Create a histogram, use twice as many tdcs to distinguish views
267  TH2F* hist = new TH2F("RegPixelMap", ";Wire;Tdc", fNWire, 0, fNWire,
268  fNTdc*3, 0, fNTdc*3);
269 
270  for(unsigned int iWire = 0; iWire < fNWire; ++iWire)
271  {
272  for(unsigned int iTdc = 0; iTdc < fNTdc; ++iTdc)
273  {
274  // Add 1 to in each bin to skip underflow
275  hist->SetBinContent(iWire+1, iTdc + fNTdc*(iWire%3) + 1,
276  (double)fLab[LocalToIndex(iWire, iTdc)]);
277 
278  }
279  }
280  return hist;
281  }
282 
283  TH2F* RegPixelMap::SingleViewToTH2(const unsigned int& view) const
284  {
285 
286  // Create a histogram
287  TH2F* hist = new TH2F("RegPixelMap", ";Wire;Tdc", fNWire, 0, fNWire,
288  fNTdc, 0, fNTdc);
289 
290  for(unsigned int iWire = 0; iWire < fNWire; ++iWire)
291  {
292  for(unsigned int iTdc = 0; iTdc < fNTdc; ++iTdc)
293  {
294  // Add 1 to in each bin to skip underflow
295  if(view==0){
296  hist->SetBinContent(iWire+1, iTdc + 1,
297  fPEX[LocalToIndex(iWire, iTdc)]);
298  // FIXIT
299  //if(fPEX[LocalToIndex(iWire, iTdc)]!=0) {
300  // std::cout<<iWire<<", "<<iTdc<<", "<<fPEX[LocalToIndex(iWire, iTdc)]/500<<std::endl;
301  //}
302  }
303  if(view==1){
304  hist->SetBinContent(iWire+1, iTdc + 1,
305  fPEY[LocalToIndex(iWire, iTdc)]);
306  }
307  if(view==2){
308  hist->SetBinContent(iWire+1, iTdc + 1,
309  fPEZ[LocalToIndex(iWire, iTdc)]);
310  }
311  }
312  }
313  return hist;
314  }
315 
316  std::ostream& operator<<(std::ostream& os, const RegPixelMap& m)
317  {
318  os << "RegPixelMap with " << m.NPixel() << " pixels, "
319  << m.NWire() << " wires"
320  << " by " << m.NTdc() << " tdcs" ;
321  return os;
322  }
323 }
int FirstWire(const unsigned int &view) const
unsigned int GlobalToIndex(const int &wire, const int &tdc, const unsigned int &view)
Take global wire, tdc (detector) and return index in fPE vector.
std::vector< HitType > fLab
Vector of Truth labels for pixels.
Definition: RegPixelMap.h:99
unsigned int GlobalToIndexSingle(const int &wire, const int &tdc, const unsigned int &view)
Take global wire, tdc (detector) and return index in fPE vector.
unsigned int fNTdc
Number of tdcs, width of pixel map.
Definition: RegPixelMap.h:86
RegPixelMap, basic input to CNN neural net.
Definition: RegPixelMap.h:22
std::vector< int > fProngTagY
Definition: RegPixelMap.h:107
std::vector< int > fProngTagZ
Definition: RegPixelMap.h:108
unsigned int fInPM
Definition: RegPixelMap.h:88
constexpr T pow(T x)
Definition: pow.h:72
unsigned int fNWire
Number of wires, length of pixel map.
Definition: RegPixelMap.h:84
std::vector< float > fPEZ
Vector of Y PE measurements for pixels.
Definition: RegPixelMap.h:94
unsigned int fNWRes
Definition: RegPixelMap.h:85
enum cnn::HType HitType
RegCNNBoundary fBound
Definition: RegPixelMap.h:104
void Add(const int &wire, const int &tdc, const unsigned int &view, const double &pe, const unsigned int &tpc, int hit_prong_tag)
Definition: RegPixelMap.cxx:60
std::vector< HitType > fLabY
Vector of Y Truth labels for pixels.
Definition: RegPixelMap.h:101
bool IsWithin(const int &wire, const int &tdc, const unsigned int &view)
unsigned int NPixel() const
Total number of pixels in map.
Definition: RegPixelMap.h:41
std::vector< double > fPurY
Vector of Y purity for pixels.
Definition: RegPixelMap.h:97
TH2F * ToLabTH2() const
static int input(void)
Definition: code.cpp:15695
TH2F * SingleViewToTH2(const unsigned int &view) const
std::vector< float > fPEX
Vector of X PE measurements for pixels.
Definition: RegPixelMap.h:92
std::vector< float > fPE
Vector of PE measurements for pixels.
Definition: RegPixelMap.h:91
std::vector< double > fPurX
Vector of X purity for pixels.
Definition: RegPixelMap.h:96
std::vector< float > fPEY
Vector of Y PE measurements for pixels.
Definition: RegPixelMap.h:93
unsigned int NTdc() const
Width in tdcs.
Definition: RegPixelMap.h:35
unsigned int fTPC
Definition: RegPixelMap.h:89
unsigned int NWire() const
Length in wires.
Definition: RegPixelMap.h:29
Defines an enumeration for cellhit classification.
std::vector< HitType > fLabX
Vector of X Truth labels for pixels.
Definition: RegPixelMap.h:100
int LastTDC(const unsigned int &view) const
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
std::vector< double > fPurZ
Vector of Y purity for pixels.
Definition: RegPixelMap.h:98
int LastWire(const unsigned int &view) const
std::ostream & operator<<(std::ostream &os, const RegPixelMap3DProducer &p)
std::vector< HitType > fLabZ
Vector of Y Truth labels for pixels.
Definition: RegPixelMap.h:102
unsigned int fNTRes
Definition: RegPixelMap.h:87
void GetTPC(const int &wire, const int &tdc, const unsigned int &view, const unsigned int &tpc)
RegPixelMap for RegCNN modified from PixelMap.h.
TH2F * ToTH2() const
Return the pixel map as a 2D histogram for visualization.
int FirstTDC(const unsigned int &view) const
void FillInputVector(float *input) const
Definition: RegPixelMap.cxx:48
std::vector< int > fProngTagX
Definition: RegPixelMap.h:106
QTextStream & endl(QTextStream &s)
unsigned int LocalToIndex(const unsigned int &wire, const unsigned int &tdc) const
Take local wire, tdc (within map) and return index in fPE vector.
std::vector< double > fPur
Vector of purity for pixels.
Definition: RegPixelMap.h:95