cellShifts.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TChain.h"
3 #include "TCut.h"
4 #include "TLine.h"
5 #include "TFile.h"
6 #include "TF2.h"
7 #include "TF1.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 #include "TStyle.h"
11 #include <cmath>
12 
13 const int pmax = 200; // max number of planes
14 const int cmax = 96; // max number of cells
15 const int nbins = 200;
16 
17 TFile* file = 0;
18 TChain* al = 0;
19 TH1D *fResidV[pmax][cmax] = {};
20 TH1D *fResidZ[pmax][cmax] = {};
21 TH2D *hCellPos[pmax][cmax] = {};
22 TH2D *p1X[pmax] = {};
23 
24 void FitForMinimum(TH2 *h2, double &x0, double &y0);
25 
26 void cellPos()
27 {
28  initialise(0); // 0 = MC, 1 = Data
29  plot_residPerViewV(0); // 0 = MC, 1 = Data
31 }
32 
33 void initialise(bool IsData)
34 {
35  file = TFile::Open("/nova/app/users/gsdavies/art/align/align_1X_cosmicsMC.root");
36  file->cd("alignment");
37 
38  int tot = 0;
39 
40  // Read in histograms from file
41  for (int plane = 0; plane < pmax; plane++){
42  for (int cell =0; cell < cmax; cell++){
43  TString hname;
44  TString hname2;
45  TString hname3;
46  if(IsData) {
47  hname = Form("residualmcV_%03d_%03d",plane,cell);
48  hname2 = Form("residualmcZ_%03d_%03d",plane,cell);
49  hname3 = Form("hCellPos_%03d_%03d",plane,cell);
50  }
51  else {
52  hname = Form("residualmcV_%03d_%03d",plane,cell);
53  hname2 = Form("residualmcZ_%03d_%03d",plane,cell);
54  hname3 = Form("hCellPos_%03d_%03d",plane,cell);
55  }
56  TKey *key = gDirectory->FindKey(hname);
57  TKey *key2 = gDirectory->FindKey(hname2);
58  TKey *key3 = gDirectory->FindKey(hname3);
59  if(key==0 || key2 == 0 || key3 == 0) {
60  //std::cerr << "!!Histogram does not exist!!" << endl;
61  tot++;
62  continue;
63  }
64  else {
65  fResidV[plane][cell] = (TH1D *)gDirectory->Get(hname);
66  fResidZ[plane][cell] = (TH1D *)gDirectory->Get(hname2);
67  hCellPos[plane][cell] = (TH2D *)gDirectory->Get(hname3);
68  }
69  }
70  std::cout << "Processing plane : " << plane << std::endl;
71  }
72  std::cout << "Total histograms that do not exist: " << tot << std::endl;
73 }
74 
75 // Plot residual V for each view
76 void plot_residPerViewV(bool IsData)
77 {
78  TCanvas* mX = new TCanvas("mX","mX");
79 
80  mX->SetBorderMode(0);
81  TCanvas* mY = new TCanvas("mY","mY");
82  mY->SetBorderMode(0);
83 
84  TH2D *meanX = new TH2D("meanX",";plane;cell",nbins,0.,199,65,0.,64);
85  TH2D *meanY = new TH2D("meanY",";plane;cell",nbins,0.,199,96,0.,95);
86  for (int plane = 0; plane < pmax; plane++){
87  cout << "Current plane in v mode: " << plane << endl;
88  for (int cell = 0; cell < cmax; cell++){
89 
90  if(hCellPos[plane][cell]){
91 
92  double v0, z0;
93 
94  FitForMinimum(hCellPos[plane][cell],v0,z0);
95 
96  // hack for boundary of two vertical planes in NDOS
97  if ((plane < 155 && plane % 2 == 0) || (plane >= 155 && plane % 2 != 0)){
98  //std::cout << "Plane/cell " << plane << ", " << cell << std::endl;
99  meanX->Fill(plane,cell,v0);
100  }
101  else if ((plane <= 155 && plane % 2 != 0) || (plane > 155 && plane % 2 == 0)){
102  meanY->Fill(plane,cell,v0);
103  }
104  }
105  else {
106  }
107  }
108  }
109 
110  // Clean up channels that don't exist in each view so they show up blank!
111  for(int i =0; i< nbins+1; i++){
112  for(int j=0; j< 66;j++){
113  if(meanX->GetBinContent(i,j) == 0){
114  meanX->SetBinContent(i,j,-999);
115  }
116  }
117  for(int k=0; k< 97;k++){
118  if(meanY->GetBinContent(i,k) == 0){
119  meanY->SetBinContent(i,k,-999);
120  }
121  }
122  }
123 
124  mX->Clear();
125  mX->cd();
126 
127  meanX->Draw("colz");
128  mY->Clear();
129  mY->cd();
130 
131  meanY->Draw("colz");
132  perView_labelV(IsData);
133 }
134 
136 {
137  TCanvas* sumXV = new TCanvas("sumXV","sumXV");
138  sumXV->SetBorderMode(0);
139  TCanvas* sumYV = new TCanvas("sumYV","sumYV");
140  sumYV->SetBorderMode(0);
141  TString nameXV = Form("rsumXV_cell_%03d",cell);
142  TString nameYV = Form("rsumYV_cell_%03d",cell);
143  TString strV = Form(" sum cell %03d",cell);
144  TH1D *SumXV = new TH1D(nameXV,";Residual (cm)"+strV,1000,-50,50);
145  TH1D *SumYV = new TH1D(nameYV,";Residual (cm)"+strV,1000,-50,50);
146 
147  for (int plane = 0; plane < pmax; plane++){
148  if(fResidV[plane][cell]){
149  // hack for boundary of two vertical planes in NDOS
150  if ((plane < 155 && plane % 2 == 0) || (plane >= 155 && plane % 2 != 0)){
151  SumXV->Add(fResidV[plane][cell]);
152  }
153  else if ((plane <= 155 && plane % 2 != 0) || (plane > 155 && plane % 2 == 0)){
154  SumYV->Add(fResidV[plane][cell]);
155  }
156  }
157  else {
158  }
159  }
160  sumXV->cd();
161  sumXV->Clear();
162  SumXV->Draw();
163 
164  sumYV->cd();
165  sumYV->Clear();
166  SumYV->Draw();
167 }
168 
170 {
171  TCanvas* sumXZ = new TCanvas("sumXZ","sumXZ");
172  sumXZ->SetBorderMode(0);
173  TCanvas* sumYZ = new TCanvas("sumYZ","sumYZ");
174  sumYZ->SetBorderMode(0);
175  TString nameXZ = Form("rsumXZ_cell_%03d",cell);
176  TString nameYZ = Form("rsumYZ_cell_%03d",cell);
177  TString strZ = Form(" sum cell %03d",cell);
178  TH1D *SumXZ = new TH1D(nameXZ,";Residual (cm)"+strZ,1000,-50,50);
179  TH1D *SumYZ = new TH1D(nameYZ,";Residual (cm)"+strZ,1000,-50,50);
180 
181  for (int plane = 0; plane < pmax; plane++){
182  if(fResidZ[plane][cell]){
183  // hack for boundary of two vertical planes in NDOS
184  if ((plane < 155 && plane % 2 == 0) || (plane >= 155 && plane % 2 != 0)){
185  SumXZ->Fill(fResidZ[plane][cell]);
186  }
187  else if ((plane <= 155 && plane % 2 != 0) || (plane > 155 && plane % 2 == 0)){
188  SumYZ->Fill(fResidZ[plane][cell]);
189  }
190  }
191  else {
192  }
193  }
194  SumXZ->Draw();
195  SumYZ->Draw();
196 }
197 
198 void plot_residPerViewZ(bool IsData)
199 {
200  TCanvas* mXz = new TCanvas("mXz","mXz");
201 
202  mXz->SetBorderMode(0);
203  TCanvas* mYz = new TCanvas("mYz","mYz");
204  mYz->SetBorderMode(0);
205 
206  TH2D *meanXz = new TH2D("meanXz",";plane;cell",nbins,0.,199,65,0.,64);
207  TH2D *meanYz = new TH2D("meanYz",";plane;cell",nbins,0.,199,96,0.,95);
208 
209  for (int plane = 0; plane < pmax; plane++){
210  cout << "Current plane in z mode: " << plane << endl;
211  for (int cell = 0; cell < cmax; cell++){
212 
213  if(hCellPos[plane][cell]){
214 
215  double v0, z0;
216  FitForMinimum(hCellPos[plane][cell],v0,z0);
217 
218  // hack for boundary of two vertical planes in NDOS
219  if ((plane < 155 && plane % 2 == 0) || (plane >= 155 && plane % 2 != 0)){
220  //std::cout << "Plane/cell " << plane << ", " << cell << std::endl;
221  meanXz->Fill(plane,cell,z0);
222  }
223  else if ((plane <= 155 && plane % 2 != 0) || (plane > 155 && plane % 2 == 0)){
224  meanYz->Fill(plane,cell,z0);
225  }
226  }
227  else {
228  }
229  }
230  }
231 
232  // Clean up channels that don't exist in each view so they show up blank!
233  for(int i =0; i< nbins+1; i++){
234  for(int j=0; j< 66;j++){
235  if(meanXz->GetBinContent(i,j) == 0){
236  meanXz->SetBinContent(i,j,-999);
237  }
238  }
239  for(int k=0; k< 97;k++){
240  if(meanYz->GetBinContent(i,k) == 0){
241  meanYz->SetBinContent(i,k,-999);
242  }
243  }
244  }
245  mXz->Clear();
246  mXz->cd();
247 
248  meanXz->Draw("colz");
249  mYz->Clear();
250  mYz->cd();
251 
252  meanYz->Draw("colz");
253  perView_labelZ(IsData);
254 
255 }
256 // Can call this once histograms are loaded to plot specific plane/cell
257 void plot_resid(bool IsData, int p, int c)
258 {
259  TCanvas* h = new TCanvas("h","h");
260  h->SetBorderMode(0);
261  h->Clear();
262  h->cd();
263  fResid[p][c]->Draw();// bin 111,9 -> plane 110, cell 8
264 
265  plane_label(IsData,p,c);
266 
267 }
268 
269 
270 void plot_residV(bool IsData, int p, int c)
271 {
272  TCanvas* h1 = new TCanvas("h1","h1");
273  h1->SetBorderMode(0);
274  h1->Clear();
275  h1->cd();
276  fResidV[p][c]->Draw();// bin 111,9 -> plane 110, cell 8
277 
278  //plane_label(IsData,p,c);
279 
280 }
281 
282 void plot_residZ(bool IsData, int p, int c)
283 {
284  TCanvas* h2 = new TCanvas("h2","h2");
285  h2->SetBorderMode(0);
286  h2->Clear();
287  h2->cd();
288  fResidZ[p][c]->Draw();// bin 111,9 -> plane 110, cell 8
289 
290  //plane_label(IsData,p,c);
291 
292 }
293 void plot_CellPos(bool IsData, int p, int c)
294 {
295  TCanvas* h3 = new TCanvas("h3","h3");
296  h3->SetBorderMode(0);
297  h3->Clear();
298  h3->cd();
299  hCellPos[p][c]->Draw("colz");// bin 111,9 -> plane 110, cell 8
300 
301  //plane_label(IsData,p,c);
302 
303 }
304 void plot_residPerPlane(int Plane){
305  TCanvas* p1 = new TCanvas("p1","p1");
306  p1->SetBorderMode(0);
307  p1->Divide(10,20);
308 
309  TCanvas* p2 = new TCanvas("p1","p1");
310  p2->SetBorderMode(0);
311  p2->Divide(10,20);
312  for (int plane = 0; plane < pmax; plane++){
313  TString name = Form("plane_%03d",plane);
314  p1X[plane] = new TH2D(name,";cell;residual v (cm)",96,0.,95,1000,-50,50);
315  p2X[plane] = new TH2D(name,";cell;residual z (cm)",96,0.,95,1000,-50,50);
316  double resid;
317  for (int cell = 0; cell < cmax; cell++){
318  if(fResidV[plane][cell]){
319  if ((plane < 155 && plane % 2 == 0) || (plane >= 155 && plane % 2 != 0)){
320 
321  for (int j = 0; j < 1001; j++){
322  resid = (j/5)-100;
323  p1X[plane]->Fill(cell,resid,fResidV[plane][cell]->GetBinContent(j));
324  }
325  }
326  else if ((plane <= 155 && plane % 2 != 0) || (plane > 155 && plane % 2 == 0)){
327  for (int j = 0; j < 1001; j++){
328  resid = (j/5)-100;
329  p1X[plane]->Fill(cell,resid,fResidV[plane][cell]->GetBinContent(j));
330  }
331  }
332  }
333  if(fResidZ[plane][cell]){
334  if ((plane < 155 && plane % 2 == 0) || (plane >= 155 && plane % 2 != 0)){
335 
336  for (int j = 0; j < 1001; j++){
337  resid = (j/5)-100;
338  p2X[plane]->Fill(cell,resid,fResidZ[plane][cell]->GetBinContent(j));
339  }
340  }
341  else if ((plane <= 155 && plane % 2 != 0) || (plane > 155 && plane % 2 == 0)){
342  for (int j = 0; j < 1001; j++){
343  resid = (j/5)-100;
344  p2X[plane]->Fill(cell,resid,fResidZ[plane][cell]->GetBinContent(j));
345  }
346  }
347  }
348  else {
349  //std::cerr << "!!Histogram not filled or does not exist!!" << endl;
350  }
351  }
352  }
353  p1->cd();
354  p1X[Plane]->Draw("colz");
355  p1->cd();
356  p2X[Plane]->Draw("colz");
357 
358 }
359 // Put a "XZ- YZ-view" tag in the corner
360 void perView_labelV(bool IsData)
361 {
362 
363  mX->cd();
364  TLatex* xz;
365  if(IsData) { xz = new TLatex(.9, .95, "XZ-View [V](Data)"); }
366  else { xz = new TLatex(.9, .95, "XZ-View [V](MC)"); }
367  xz->SetTextColor(kBlue);
368  xz->SetNDC();
369  xz->SetTextSize(2/30.);
370  xz->SetTextAlign(32);
371  xz->Draw();
372 
373  mY->cd();
374  TLatex* yz;
375  if(IsData) { yz = new TLatex(.9, .95, "YZ-View [V](Data)"); }
376  else { yz = new TLatex(.9, .95, "YZ-View [V](MC)"); }
377  yz->SetTextColor(kBlue);
378  yz->SetNDC();
379  yz->SetTextSize(2/30.);
380  yz->SetTextAlign(32);
381  yz->Draw();
382 }
383 void perView_labelZ(bool IsData)
384 {
385 
386  mXz->cd();
387  TLatex* xz;
388  if(IsData) { xz = new TLatex(.9, .95, "XZ-View [Z](Data)"); }
389  else { xz = new TLatex(.9, .95, "XZ-View [Z](MC)"); }
390  xz->SetTextColor(kBlue);
391  xz->SetNDC();
392  xz->SetTextSize(2/30.);
393  xz->SetTextAlign(32);
394  xz->Draw();
395 
396  mYz->cd();
397  TLatex* yz;
398  if(IsData) { yz = new TLatex(.9, .95, "YZ-View [Z](Data)"); }
399  else { yz = new TLatex(.9, .95, "YZ-View [Z](MC)"); }
400  yz->SetTextColor(kBlue);
401  yz->SetNDC();
402  yz->SetTextSize(2/30.);
403  yz->SetTextAlign(32);
404  yz->Draw();
405 }
406 
407 void plane_label(bool IsData, int p, int c )
408 {
409  h->cd();
410  if(IsData) {
411  TString title = Form("Plane %03d, Cell %03d (Data)",p,c);
412  TLatex* hh = new TLatex(.9, .95, title);
413  }
414  else {
415  TString title = Form("Plane %03d, Cell %03d (MC)",p,c);
416  TLatex* hh = new TLatex(.9, .95, title);
417  }
418  hh->SetTextColor(kBlue);
419  hh->SetNDC();
420  hh->SetTextSize(2/30.);
421  hh->SetTextAlign(32);
422  hh->Draw();
423 
424 }
425 
426 // This method is to redraw 2D plots to white out bins with values less than or greater
427 // than the minmax ( specifically addresses the -999 values that take over the
428 // colz range
429 void SetLimits(double minmax)
430 {
431  mX->cd();
432  meanX->SetMaximum(minmax);
433  meanX->SetMinimum(-minmax);
434  meanX->Draw("colz");
435 
436  mY->cd();
437  meanY->SetMaximum(minmax);
438  meanY->SetMinimum(-minmax);
439  meanY->Draw("colz");
440 
441  mXz->cd();
442  meanXz->SetMaximum(minmax);
443  meanXz->SetMinimum(-minmax);
444  meanXz->Draw("colz");
445 
446  mYz->cd();
447  meanYz->SetMaximum(minmax);
448  meanYz->SetMinimum(-minmax);
449  meanYz->Draw("colz");
450 }
451 
452 // Save the histograms!
453 void save(bool IsData, bool Bad)
454 {
455  if(IsData) {
456  mX->SaveAs("xzView-resid_data.png");
457  mY->SaveAs("yzView-residual_data.png");
458  }
459  else {
460  if(Bad) {
461  mX->SaveAs("xzView-resid_mcbad.png");
462  mY->SaveAs("yzView-residual_mcbad.png");
463  }
464  else {
465  mX->SaveAs("xzView-resid_mc.png");
466  mY->SaveAs("yzView-residual_mc.png");
467  }
468  }
469 }
470 
471 // Bulk of the D2M method is done here! Fits for the minimum value
472 void FitForMinimum(TH2 *h2, double &x0, double &y0) {
473 
474  // for speed, only fit right around minimum
475 
476  int binX0;
477  int binY0;
478  int binZ0;
479  h2->GetMinimumBin(binX0,binY0,binZ0);
480 
481  int fitbinX1 = binX0 - 2;
482  int fitbinX2 = binX0 + 2;
483  int fitbinY1 = binY0 - 2;
484  int fitbinY2 = binY0 + 2;
485 
486  if (fitbinX1<1) {
487  fitbinX2 += (1-fitbinX1);
488  fitbinX1 += (1-fitbinX1);
489  }
490  if (fitbinX2>h2->GetNbinsX()) {
491  fitbinX1 -= (fitbinX2-h2->GetNbinsX());
492  fitbinX2 -= (fitbinX2-h2->GetNbinsX());
493  }
494  if (fitbinY1<1) {
495  fitbinY2 += (1-fitbinY1);
496  fitbinY1 += (1-fitbinY1);
497  }
498  if (fitbinY2>h2->GetNbinsY()) {
499  fitbinX1 += (fitbinX2-h2->GetNbinsY());
500  fitbinX2 += (fitbinX2-h2->GetNbinsY());
501  }
502 
503  TF2 func("func","[2]*(x-[0])*(x-[0])+[3]*(x-[0])*(y-[1])+[4]*(y-[1])*(y-[1])+[5]",
504  h2->GetXaxis()->GetBinLowEdge(fitbinX1),
505  h2->GetXaxis()->GetBinUpEdge (fitbinX2),
506  h2->GetYaxis()->GetBinLowEdge(fitbinY1),
507  h2->GetYaxis()->GetBinUpEdge (fitbinY2));
508 
509  func.SetParameters(
510  h2->GetXaxis()->GetBinCenter(binX0),
511  h2->GetYaxis()->GetBinCenter(binY0),
512  1000,-10,1000,h2->GetBinContent(binX0,binY0));
513 
514  h2->Fit(&func,"q0RN","goff");
515  x0 = func.GetParameter(0);
516  y0 = func.GetParameter(1);
517 }
const int pmax
Definition: cellShifts.C:13
const XML_Char * name
Definition: expat.h:151
void plot_CellPos(bool IsData, int p, int c)
Definition: cellShifts.C:293
TH1D * fResidZ[pmax][cmax]
Definition: cellShifts.C:20
void plot_residV(bool IsData, int p, int c)
Definition: cellShifts.C:270
TH2D * hCellPos[pmax][cmax]
Definition: cellShifts.C:21
void plot_residPerViewZ(bool IsData)
Definition: cellShifts.C:198
const int cmax
Definition: cellShifts.C:14
TH1F * h3
Definition: berger.C:36
const char * p
Definition: xmltok.h:285
void perView_labelV(bool IsData)
Definition: cellShifts.C:360
void FitForMinimum(TH2 *h2, double &x0, double &y0)
Definition: cellShifts.C:472
void initialise(bool IsData)
Definition: cellShifts.C:33
void plot_residPerPlane(int Plane)
Definition: cellShifts.C:304
const int nbins
Definition: cellShifts.C:15
void plot_residZ(bool IsData, int p, int c)
Definition: cellShifts.C:282
TH2D * p1X[pmax]
Definition: cellShifts.C:22
void perView_labelZ(bool IsData)
Definition: cellShifts.C:383
void cellPos()
Definition: cellShifts.C:26
double func(double x, double y)
const double j
Definition: BetheBloch.cxx:29
void SetLimits(double minmax)
Definition: cellShifts.C:429
void save(bool IsData, bool Bad)
Definition: cellShifts.C:453
void plot_sumcellsV(int cell)
Definition: cellShifts.C:135
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
TH1F * h1
void plot_resid(bool IsData, int p, int c)
Definition: cellShifts.C:257
TChain * al
Definition: cellShifts.C:18
void plane_label(bool IsData, int p, int c)
Definition: cellShifts.C:407
TFile * file
Definition: cellShifts.C:17
TH1D * fResidV[pmax][cmax]
Definition: cellShifts.C:19
void plot_residPerViewV(bool IsData)
Definition: cellShifts.C:76
void plot_sumcellsZ(int cell)
Definition: cellShifts.C:169
hh[ixs]
Definition: PlotSingle.C:6