Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FillBlobHist.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FillBlobHist.C
1 
2 // FillBlobHist.C //
4 
5 // This makes a BH (Blob Bistogram) that can show the raw hits.
6 
7 #include "FillBlobHist.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 #include "ABlob.h"
11 #include "AZigzag.h"
12 
13 #include "groot.h"
14 
15 #include <iostream>
16 #include <cmath>
17 #include <string>
18 
19 TH1D* BH_Q = 0; // How many blobs of a certain charge
20 TH1D* BH_Size = 0; // How many blobs made of a certain # of zig-zag pads
21 TH2D* BH_QvsSize = 0; // Charge of a blob vs. how many pads it's made of
22 
23 TH1D* BH_PhiTic = 0; // Just the phi location of each pad
24 TH2D* BH_PhivsR = 0; // Phi vs. Radius
25 
26 TH1D* BH_Tmax = 0; // This is the time of the zig with the most charge
27 TH1D* BH_Tdiff = 0;
28 
29 TH1D* BH_5Residual = 0; // Minor cheating based on arc-length for a quick preliminary resolution
30 TH2D* BH_5Residual_vsPhi = 0;
31  TH1D* BH_5Residual_PC = 0;
32  TH2D* BH_5Residual_vsPhi_PC = 0;
33  TH1D* BH_5Efficiency = 0; // Checks if a blob even exists on the 5th radii if 4 & 6 are hit
34 
35 
36 // Histograming ALL Radii through an array
37 char name[100];
38 TH1D* BH_Residual_[16]; // Initialize rows 0 and 15 (in case they're useful later)
39 TH2D* BH_Residual_vsPhi_[16];
40  TH1D* BH_Residual_PC_[16]; // PC = PhiCut (efficiency/track attempt before Hough)
41  TH2D* BH_Residual_vsPhi_PC_[16];
42  TH1D* BH_Efficiency_[16];
43 
44 
45 using namespace std;
46 
47 void FillBlobHist()
48 {
49  groot* Tree = groot::instance();
50 
51  // phi=90 and phi=91 are instrumented (connected to hybrids)
52  double phi90 = Tree->ZigzagMap2[0][90]->PCenter();
53  double phi91 = Tree->ZigzagMap2[0][91]->PCenter();
54 
55  double phi0 = 91.0*phi90 - 90.0*phi91;
56  double phi127 = phi0 + 127.0*(phi91-phi90);
57 
58  double BIN_FACTOR = 10.0;
59  double dPhi_BIN = (phi91-phi90)/BIN_FACTOR;
60 
61  double min = phi0 - dPhi_BIN/2.0;
62  double max = phi127 + dPhi_BIN/2.0;
63  int N_BIN = 127*BIN_FACTOR;
64 
65  if (!BH_Q)
66  {
67  BH_Q = new TH1D("BH_Q", "BH_Q", 256, -0.5, 2047.5);
68  BH_Size = new TH1D("BH_Size", "BH_Size", 21, -0.5, 20.5);
69  BH_QvsSize = new TH2D("BH_QvsSize", "BH_QvsSize", 256, -0.5, 2047.5, 21, -0.5, 20.5);
70  BH_PhiTic = new TH1D("BH_PhiTic", "BH_PhiTic", 2*N_BIN, min, max); // "normalized" tics
71  BH_PhivsR = new TH2D("BH_PhivsR", "BH_PhivsR", N_BIN, min, max, 16, -0.5, 15.5);
72  BH_Tmax = new TH1D("BH_Tmax", "BH_Tmax", 300, -5.0, 30.0); // "normalized" tics
73  BH_Tdiff = new TH1D("BH_Tdiff", "BH_Tdiff", 20000, -30.0, 30.0); // "normalized" tic difference
74 
75  // Resolution-study Histograms
76  BH_5Residual = new TH1D("BH_5Residual", "BH_5Residual", 200, -1000.0, 1000.0); //arc length in microns
77  BH_5Residual_vsPhi = new TH2D("BH_5Residual_vsPhi", "BH_5Residual_vsPhi", (int)(0.2/dPhi_BIN), 1.50, 1.66, 200, -1300.0, 1300.0);
78  BH_5Residual_PC = new TH1D("BH_5Residual_PC", "BH_5Residual_PC", 200, -1000.0, 1000.0); //arc length in microns
79  BH_5Residual_vsPhi_PC = new TH2D("BH_5Residual_vsPhi_PC", "BH_5Residual_vsPhi_PC", (int)(0.2/dPhi_BIN), 1.53, 1.63, 200, -1000.0, 1000.0);
80  BH_5Efficiency = new TH1D("BH_5Efficiency", "BH_5Efficiency", 2, 0, 2);
81 
82 
83  // Histograming ALL Radii through an array
84  for (int i=1; i<15; i++) // Can't do radii 0 or 15 since we don't have radii -1 or 16 for comparison
85  {
86  sprintf(name, "BH_Residual_%d", i);
87  BH_Residual_[i] = new TH1D(name, name, 200,-1000.0, 1000.0);
88  sprintf(name, "BH_Residual_vsPhi_%d", i);
89  BH_Residual_vsPhi_[i] = new TH2D(name, name, (int)(0.2/dPhi_BIN), 1.50, 1.66, 200, -1300.0, 1300.0);
90  sprintf(name, "BH_Residual_PC_%d", i);
91  BH_Residual_PC_[i] = new TH1D(name, name, 200,-1000.0, 1000.0);
92  sprintf(name, "BH_Residual_vsPhi_PC_%d", i);
93  BH_Residual_vsPhi_PC_[i] = new TH2D(name, name, (int)(0.2/dPhi_BIN), 1.50, 1.66, 200, -1300.0, 1300.0);
94  sprintf(name, "BH_Efficiency_%d", i);
95  BH_Efficiency_[i] = new TH1D(name, name, 2, 0, 2);
96  }
97 
98 
99  // Used for Condor
100  Tree->theHistograms.push_back(BH_Q);
101  Tree->theHistograms.push_back(BH_Size);
102  Tree->theHistograms.push_back(BH_QvsSize);
103  Tree->theHistograms.push_back(BH_PhiTic );
104  Tree->theHistograms.push_back(BH_PhivsR);
105  Tree->theHistograms.push_back(BH_Tmax);
106  Tree->theHistograms.push_back(BH_Tdiff);
107  Tree->theHistograms.push_back(BH_5Residual);
108  Tree->theHistograms.push_back(BH_5Residual_PC);
109  Tree->theHistograms.push_back(BH_5Residual_vsPhi);
110  Tree->theHistograms.push_back(BH_5Residual_vsPhi_PC);
111  Tree->theHistograms.push_back(BH_5Efficiency);
112  for (int i=1; i<15; i++)
113  {
114  Tree->theHistograms.push_back(BH_Residual_[i]);
115  Tree->theHistograms.push_back(BH_Residual_vsPhi_[i]);
116  Tree->theHistograms.push_back(BH_Residual_PC_[i]);
117  Tree->theHistograms.push_back(BH_Residual_vsPhi_PC_[i]);
118  Tree->theHistograms.push_back(BH_Efficiency_[i]);
119  }
120  } // End of creating histograms
121 
122 
123  if ( BH_PhiTic->GetMaximum() != 1) //Gets the positions of the pads. Hist hight 0 to give better BH_->Scale(int) control
124  {
125  for (int i=0; i<Nphi; i++)
126  {
127  if (Tree->ZigzagMap2[0][i] ) { BH_PhiTic->Fill(Tree->ZigzagMap2[0][i]->PCenter()); }
128  }
129  }
130 
131  // All radii are in mm. Only convert to microns in resolution since all x,y values are also in mm
132  double RadiusValue[16];
133  for (int k=0; k<16; k++) { RadiusValue[k] = Tree->theZigzags[k]->RCenter(); }
134 
135 
136  for (int i=0; i<Nr; i++)
137  {
138  for (int j=0; j< Tree->theBlobs[i].size(); j++)
139  {
140  double Q = Tree->theBlobs[i][j]->Q();
141  double Phi = Tree->theBlobs[i][j]->CentroidP();
142  int Radius = i;
143 
144  int numZigs = Tree->theBlobs[i][j]->numZigs();
145 
146  if ( Tree->theBlobs[i].size()==1 ) // Only 1 blob per radius
147  {
148  BH_Q->Fill(Q);
149  BH_Size->Fill(numZigs);
150  BH_QvsSize->Fill(Q, numZigs);
151 
152  BH_PhivsR->Fill(Phi, Radius);
153 
154  double TMAX = Tree->theBlobs[i][0]->maxT();
155  BH_Tmax->Fill(TMAX);
156  for (int k=0; k< Tree->theBlobs[i][0]->manyZigs.size(); k++)
157  {
158  BH_Tdiff->Fill( Tree->theBlobs[i][0]->manyZigs[k]->T() - TMAX );
159  }
160  }
161  }
162  } // End of priliminary histograms
163 
164 
165 
166  // Original full Resolution studies (over all angles)
167  if ( Tree->theBlobs[3].size()==1 && Tree->theBlobs[4].size()==1 && Tree->theBlobs[5].size()==1 )
168  {
169  // Calculation to handle polar coordinates correctly
170  double X4 = Tree->theBlobs[3][0]->CentroidX(), Y4 = Tree->theBlobs[3][0]->CentroidY();
171  double X6 = Tree->theBlobs[5][0]->CentroidX(), Y6 = Tree->theBlobs[5][0]->CentroidY();
172 
173  double x5a, y5a, x5b, y5b; //line must intersect the circle in 2 points
174  int found = line_circle(X4, Y4, X6, Y6, RadiusValue[4], x5a, y5a, x5b, y5b); //takes 9 arguments
175 
176  if (found == 2) // Nonesense unless both are found
177  {
178  double X5p = x5a, Y5p = y5a;
179  if (y5b > 0)
180  {
181  X5p = x5b; Y5p = y5b;
182  }
183 
184  double Phi5p = atan2(Y5p, X5p); // predicted/expected angle, should be around 1.58 [rad]
185  double Phi5a = Tree->theBlobs[4][0]->CentroidP(); // actual angle
186  double Residual5 = Phi5a - Phi5p; // difference in angles
187 
188  BH_5Residual->Fill(RadiusValue[4]*Residual5*1000.0); //angle * radius * conversion = arc length [microns]
189  BH_5Residual_vsPhi->Fill(Phi5p, RadiusValue[4]*Residual5*1000.0); //Fill (predicted/expected angle, arc length[microns]);
190  } // End of filling actual histographs
191  } // End of making sure only one blob per radii
192 
193 
194  // New efficiency histograms with a Phi cut
195  if ( Tree->theBlobs[3].size()==1 && Tree->theBlobs[5].size()==1 )
196  {
197  if ( (Tree->theBlobs[3][0]->CentroidP() >= 1.53) && (Tree->theBlobs[3][0]->CentroidP() <= 1.63) &&
198  (Tree->theBlobs[5][0]->CentroidP() >= 1.53) && (Tree->theBlobs[5][0]->CentroidP() <= 1.63) )
199  {
200  if ( Tree->theBlobs[4].size()==1 )
201  {
202  BH_5Efficiency->Fill(1.5);
203 
204  // Calculation to handle polar coordinates correctly
205  double X4 = Tree->theBlobs[3][0]->CentroidX(), Y4 = Tree->theBlobs[3][0]->CentroidY();
206  double X6 = Tree->theBlobs[5][0]->CentroidX(), Y6 = Tree->theBlobs[5][0]->CentroidY();
207 
208  double x5a, y5a, x5b, y5b; //line must intersect the circle in 2 points
209  int found = line_circle(X4, Y4, X6, Y6, RadiusValue[4], x5a, y5a, x5b, y5b); //takes 9 arguments
210 
211  if (found == 2) // Nonesense unless both are found
212  {
213  double X5p = x5a, Y5p = y5a;
214  if (y5b > 0)
215  {
216  X5p = x5b; Y5p = y5b;
217  }
218 
219  double Phi5p = atan2(Y5p, X5p); // predicted/expected angle, should be around 1.58 [rad]
220  double Phi5a = Tree->theBlobs[4][0]->CentroidP(); // actual angle
221  double Residual5 = Phi5a - Phi5p; // difference in angles
222 
223  BH_5Residual_PC->Fill(RadiusValue[4]*Residual5*1000.0); //angle * radius * conversion = arc length [microns]
224  BH_5Residual_vsPhi_PC->Fill(Phi5p, RadiusValue[4]*Residual5*1000.0); //Fill (predicted/expected angle, arc length[microns]);
225  } // End of filling actual histographs
226  } // End of one blob on 5
227  else { BH_5Efficiency->Fill(0.5); }
228  } // End of one blob on 4 & 6 within angle cut
229  } // End of one blob per adjacent radii
230 
231 
232 
233 
234 
235 
236 
237 
238 //---------------------------------------------------------------------------------------------------------------------------------------
239 //------------------------------------------------|ALL RESOLUTIONS|----------------------------------------------------------------------
240 //---------------------------------------------------------------------------------------------------------------------------------------
241 
242  for (int i=1; i<15; i++)
243  {
244  // Original full Resolution studies (over all angles)
245  if ( Tree->theBlobs[i-1].size()==1 && Tree->theBlobs[i].size()==1 && Tree->theBlobs[i+1].size()==1 )
246  {
247  // Calculation to handle polar coordinates correctly
248  double Xpre = Tree->theBlobs[i-1][0]->CentroidX(), Ypre = Tree->theBlobs[i-1][0]->CentroidY();
249  double Xpost = Tree->theBlobs[i+1][0]->CentroidX(), Ypost = Tree->theBlobs[i+1][0]->CentroidY();
250 
251  double xia, yia, xib, yib; //line must intersect the circle in 2 points
252  int found = line_circle(Xpre, Ypre, Xpost, Ypost, RadiusValue[i], xia, yia, xib, yib); //takes 9 arguments
253 
254  if (found == 2) // Nonesense unless both are found
255  {
256  double Xip = xia, Yip = yia;
257  if (yib > 0) { Xip = xib; Yip = yib; }
258 
259  double Phi_ip = atan2(Yip, Xip); // predicted/expected angle, should be around 1.58 [rad]
260  double Phi_ia = Tree->theBlobs[i][0]->CentroidP(); // actual angle
261  double Residuali = Phi_ia - Phi_ip; // difference in angles
262 
263  BH_Residual_[i]->Fill(RadiusValue[i]*Residuali*1000.0); //angle * radius * conversion = arc length [microns]
264  BH_Residual_vsPhi_[i]->Fill(Phi_ip, RadiusValue[i]*Residuali*1000.0); //Fill (predicted/expected angle, arc length[microns]);
265  } // End of filling actual histographs
266  } // End of making sure only one blob per radii
267 
268 
269  // New efficiency histograms with a Phi cut
270  if ( Tree->theBlobs[i-1].size()==1 && Tree->theBlobs[i+1].size()==1 )
271  {
272  double PC_LOW = 1.53;
273  double PC_HIGH = 1.63;
274  if ( (Tree->theBlobs[i-1][0]->CentroidP() >= PC_LOW) && (Tree->theBlobs[i-1][0]->CentroidP() <= PC_HIGH) &&
275  (Tree->theBlobs[i+1][0]->CentroidP() >= PC_LOW) && (Tree->theBlobs[i+1][0]->CentroidP() <= PC_HIGH) )
276  {
277  if ( Tree->theBlobs[i].size()==1 )
278  {
279  BH_Efficiency_[i]->Fill(1.5);
280 
281  // Calculation to handle polar coordinates correctly
282  double Xpre = Tree->theBlobs[i-1][0]->CentroidX(), Ypre = Tree->theBlobs[i-1][0]->CentroidY();
283  double Xpost = Tree->theBlobs[i+1][0]->CentroidX(), Ypost = Tree->theBlobs[i+1][0]->CentroidY();
284 
285  double xia, yia, xib, yib; //line must intersect the circle in 2 points
286  int found = line_circle(Xpre, Ypre, Xpost, Ypost, RadiusValue[i], xia, yia, xib, yib); //takes 9 arguments
287 
288  if (found == 2) // Nonesense unless both are found
289  {
290  double Xip = xia, Yip = yia;
291  if (yib > 0) { Xip = xib; Yip = yib; }
292 
293  double Phi_ip = atan2(Yip, Xip); // predicted/expected angle, should be around 1.58 [rad]
294  double Phi_ia = Tree->theBlobs[i][0]->CentroidP(); // actual angle
295  double Residuali = Phi_ia - Phi_ip; // difference in angles
296 
297  BH_Residual_PC_[i]->Fill(RadiusValue[i]*Residuali*1000.0); //angle * radius * conversion = arc length [microns]
298  BH_Residual_vsPhi_PC_[i]->Fill(Phi_ip, RadiusValue[i]*Residuali*1000.0); //Fill (predicted/expected angle, arc length[microns]);
299  } // End of filling actual histographs
300  } // End of one blob on 5
301  else { BH_Efficiency_[i]->Fill(0.5); }
302  } // End of one blob on 4 & 6 within angle cut
303  } // End of one blob per adjacent radii
304 
305  } // end of ARRAY of histograms
306 
307 //---------------------------------------------------------------------------------------------------------------------------------------
308 //------------------------------------------------|ALL RESOLUTIONS END|------------------------------------------------------------------
309 //---------------------------------------------------------------------------------------------------------------------------------------
310 
311 }
312 
313 
314 
315 
316 
317 
318 int line_circle(double m, double b, double r, double& Xi1, double& Yi1, double& Xi2, double& Yi2)
319 {
320  return line_circle( 0.0, b, 1.0, (m+b), r, Xi1, Yi1, Xi2, Yi2);
321 }
322 
323 
324 int line_circle(double X1, double Y1, double X2, double Y2, double r, double& Xi1, double& Yi1, double& Xi2, double& Yi2) //takes 9 arguments
325 {
326  double dx = X2 - X1;
327  double dy = Y2 - Y1;
328  double dr = sqrt(dx*dx + dy*dy);
329  double D = X1*Y2 - X2*Y1;
330 
331  double DELTA = r*r*dr*dr - D*D;
332  if (DELTA < 0) return 0;
333 
334  if (DELTA == 0)
335  {
336  double sign = (dy > 0) - (dy < 0); // Implements the sign() function
337  Xi1 = D*dy/(dr*dr);
338  Yi1 = -D*dx/(dr*dr);
339  return 1;
340  }
341 
342  double sign = (dy > 0) - (dy < 0); // Implements the sign() function
343  Xi1 = ( D*dy + sign*dx*sqrt(DELTA))/(dr*dr);
344  Yi1 = (-D*dx + abs(dy)*sqrt(DELTA))/(dr*dr);
345  Xi2 = ( D*dy - sign*dx*sqrt(DELTA))/(dr*dr);
346  Yi2 = (-D*dx - abs(dy)*sqrt(DELTA))/(dr*dr);
347  return 2;
348 }