Functions
DedxGamma.C File Reference
#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <string>
#include <cstring>
#include "TFile.h"
#include "TChain.h"
#include "TH2.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TProfile.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TGaxis.h"

Go to the source code of this file.

Functions

TH1D * MakeRatio (TH1D *num, TH1D *denom, const char *name)
 
TCanvas * MakeSplitRatio (TCanvas *top, TCanvas *bot, const char *name)
 
void FindTruncate (TH1D *hist, double &low, double &high)
 
double Mode (TH1D *hist)
 
double ModeRestrict (TH1D *hist, double minx, double maxx)
 
void GausFit (TH1D *hist, double &mean, double &err, double minrange, double maxrange, double lowsig, double highsig, bool useTruncate)
 
TH1D * PeakProjectionY (TH2F *hist, const char *name, double minrange=-1.0, double maxrange=-1.0, double lowsig=0.001, double highsig=0.001, bool useTruncate=true)
 
void DedxGamma ()
 

Function Documentation

void DedxGamma ( )

Definition at line 308 of file DedxGamma.C.

References demo5::c2, chisquared::c3, chisquared::c4, chi2(), om::cout, release_diff::diff, dist, allTimeWatchdog::endl, demo5::expt, make_syst_table_plots::ibin, kGreen, kRed, kViolet, MakeRatio(), MakeSplitRatio(), nsteps, and PeakProjectionY().

309 {
310 
311  //const char* filenamedata = "/local/nova/users/raddatz/S15-05-04_reco_pid_tuning/ReMId/DataMCDedx/datadedxVxFineS15-05-22.root";
312  //const char* filenamemc = "/local/nova/users/raddatz/S15-05-04_reco_pid_tuning/ReMId/DataMCDedx/mcdedxVxFineS15-05-22.root";
313 
314  const char* filenamedata = "/local/nova/users/lein/datadedxVgammaFineS15-05-22.root";
315  const char* filenamemc = "/local/nova/users/lein/mcdedxVgammaFineS15-05-22.root";
316 
317  const char* filenamedataP = "/local/nova/users/lein/datadedxVgammaFineS15-05-22Proton.root";
318  const char* filenamemcP = "/local/nova/users/lein/mcdedxVgammaFineS15-05-22Proton.root";
319 
320  TFile* filedata = new TFile(filenamedata,"READ");
321  TFile* filemc = new TFile(filenamemc,"READ");
322 
323  TFile* filedataP = new TFile(filenamedataP,"READ");
324  TFile* filemcP = new TFile(filenamemcP,"READ");
325 
326  TH2F* dedxdata = (TH2F*)filedata->Get("DedxVx");
327  TH2F* dedxmc = (TH2F*)filemc->Get("DedxVx");
328 
329  TH2F* dedxdataP = (TH2F*)filedataP->Get("DedxVx");
330  TH2F* dedxmcP = (TH2F*)filemcP->Get("DedxVx");
331 
332  //std::cout<<"This is the pointer: "<<dedxdata<<std::endl;
333 
334  gStyle->SetOptStat(00000);
335 
336  TH1D* dedxpeakdata = PeakProjectionY(dedxdata,"dedxpeakdata");
337  TH1D* dedxpeakmc = PeakProjectionY(dedxmc,"dedxpeakmc");
338 
339  TH1D* dedxpeakdataP = PeakProjectionY(dedxdataP,"dedxpeakdataP");
340  std::cout<<"About to do the one!"<<std::endl;
341  TH1D* dedxpeakmcP = PeakProjectionY(dedxmcP,"dedxpeakmcP");
342 
343  //If you want to run WITHOUT the truncation and just limit the fit range to +- 0.001 around the mode, do these two lines
344  //instead of the two lines above
345  //TH1D* dedxpeakdata = PeakProjectionY(dedxdata,"dedxpeakdata",-1.0,-1.0,0.001,0.001,false);
346  //TH1D* dedxpeakmc = PeakProjectionY(dedxmc,"dedxpeakmc",-1.0,-1.0,0.001,0.001,false);
347 
348 
349  TCanvas* c0 = new TCanvas("c0","c0");
350  dedxpeakdata->Draw();
351  dedxpeakmc->Draw("same");
352  dedxpeakmcP->Draw("same");
353  dedxpeakdataP->Draw("same");
354  dedxpeakmc->SetLineColor(kRed);
355  dedxpeakmc->SetMarkerColor(kRed);
356  dedxpeakmcP->SetLineColor(kViolet);
357  dedxpeakmcP->SetMarkerColor(kViolet);
358  dedxpeakdataP->SetLineColor(kGreen+2);
359  dedxpeakdataP->SetMarkerColor(kGreen+2);
360  dedxpeakdata->GetYaxis()->SetRangeUser(0.0,0.006);
361  dedxpeakdata->GetXaxis()->SetNoExponent();
362  //dedxpeakdata->GetYaxis()->SetNoExponent();
363  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
364  leg0->AddEntry(dedxpeakdata,"Muon Data","LE");
365  leg0->AddEntry(dedxpeakmc,"Muon MC","LE");
366  leg0->AddEntry(dedxpeakdataP,"Proton Data","LE");
367  leg0->AddEntry(dedxpeakmcP,"Proton MC","LE");
368  leg0->Draw("same");
369  c0->Print("dedxfine.pdf");
370  c0->Print("dedxfine.root");
371 
372  TCanvas* c00 = new TCanvas("c00","c00");
373  dedxdata->Draw("colz");
374  dedxpeakdata->Draw("same");
375  dedxdata->GetYaxis()->SetRangeUser(0.0,0.006);
376  dedxdata->GetXaxis()->SetNoExponent();
377  //dedxdata->GetYaxis()->SetNoExponent();
378  c00->Print("dedxDataMother.pdf");
379 
380  TCanvas* c000 = new TCanvas("c000","c000");
381  //dedxpeakmc->SetLineColor(kBlack);
382  dedxmc->Draw("colz");
383  dedxpeakmc->Draw("same");
384  dedxmc->GetYaxis()->SetRangeUser(0.0,0.006);
385  dedxmc->GetXaxis()->SetNoExponent();
386  //dedxmc->GetYaxis()->SetNoExponent();
387  c000->Print("dedxMCMother.pdf");
388 
389  TCanvas* c00P = new TCanvas("c00P","c00P");
390  dedxdataP->Draw("colz");
391  dedxpeakdataP->Draw("same");
392  dedxdataP->GetYaxis()->SetRangeUser(0.0,0.006);
393  dedxdataP->GetXaxis()->SetNoExponent();
394  //dedxdata->GetYaxis()->SetNoExponent();
395  c00P->Print("dedxDataProtonMother.pdf");
396 
397  TCanvas* c000P = new TCanvas("c000P","c000P");
398  //dedxpeakmc->SetLineColor(kBlack);
399  dedxmcP->Draw("colz");
400  dedxpeakmcP->Draw("same");
401  dedxmcP->GetYaxis()->SetRangeUser(0.0,0.006);
402  dedxmcP->GetXaxis()->SetNoExponent();
403  //dedxmc->GetYaxis()->SetNoExponent();
404  c000P->Print("dedxMCProtonMother.pdf");
405 
406 
407  int nsteps = 1000;
408  double minscale = 0.95;
409  double maxscale = 1.05;
410 
411  TH1F* chi2match = new TH1F("chi2match",";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
412  TH1F* chi2matchrestrict = new TH1F("chi2matchrestrict",";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
413 
414 
415  for(int ibin = 1; ibin <= chi2match->GetNbinsX(); ++ibin){
416  // Get the scale factor
417  double scalefactor = chi2match->GetXaxis()->GetBinCenter(ibin);
418  //calculate a chi2 between data and scaled mc
419  double chi2 = 0;
420  double chi2rest = 0;
421  for(int jbin = 1; jbin <= dedxpeakdata->GetNbinsX(); ++jbin){
422  // get the dist from the end of the track
423  double dist = dedxpeakdata->GetXaxis()->GetBinCenter(jbin);
424  double expt = dedxpeakdata->GetBinContent(jbin);
425  double errdata = dedxpeakdata->GetBinError(jbin);
426  double errmc = dedxpeakmc->GetBinError(jbin);
427  double err2 = errdata*errdata+errmc*errmc;
428  double obsunscale = dedxpeakmc->GetBinContent(jbin);
429  double obs = scalefactor*dedxpeakmc->GetBinContent(jbin);
430  double diff = expt-obs;
431  //chi2+=(diff*diff/expt);
432  // remove error cases
433  if(expt < 0 || obsunscale < 0){ continue; }
434  chi2+=(diff*diff/err2);
435  // std::cout<<"dist: "<<dist<<std::endl;
436  if(dist>100 && dist < 200){
437  chi2rest+=(diff*diff/err2);
438  }
439 
440  }
441 // // make a histogram to store the scaled version of mc dedxpeak histogram
442 // TString name;
443 // name+=ibin;
444 // TH1D* dedxpeakmcscale = new TH1D(name.Data(),"",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
445 // dedxpeakmcscale->Add(dedxpeakmc);
446 // dedxpeakmcscale->Scale(scalefactor);
447 // dedxpeakmcscale->
448 // double chi2 = dedxpeakdata->Chi2Test(dedxpeakmcscale,"CHI2");
449 // std::cout<<"scalefactor: "<<scalefactor<<" chi2: "<<chi2<<std::endl;
450  chi2match->SetBinContent(ibin,chi2);
451  chi2matchrestrict->SetBinContent(ibin,chi2rest);
452 
453  }
454 
455  TCanvas* c00a = new TCanvas("c00a","c00a");
456  chi2match->Draw();
457  c00a->Print("chi2fine.pdf");
458 
459  double optscale = chi2match->GetXaxis()->GetBinCenter(chi2match->GetMinimumBin());
460  std::cout<<"optscale: "<<optscale<<std::endl;
461 
462  TCanvas* c000a = new TCanvas("c000a","c000a");
463  chi2matchrestrict->Draw();
464  c000a->Print("chi2restrictfine.pdf");
465 
466  double optscalerestrict = chi2matchrestrict->GetXaxis()->GetBinCenter(chi2matchrestrict->GetMinimumBin());
467  std::cout<<"Restricted optscale: "<<optscalerestrict<<std::endl;
468 
469  TH1D* dedxpeakmcscale = new TH1D("OptimalMatch","",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
470  dedxpeakmcscale->Add(dedxpeakmc);
471  dedxpeakmcscale->Scale(optscale);
472 
473  TCanvas* c2 = new TCanvas("c2","c2");
474  dedxpeakdata->Draw();
475  dedxpeakmcscale->Draw("same");
476  dedxpeakmcscale->SetLineColor(kRed);
477  dedxpeakmcscale->SetMarkerColor(kRed);
478  dedxpeakdata->GetYaxis()->SetRangeUser(0.0014,0.0025);
479  dedxpeakdata->GetXaxis()->SetNoExponent();
480  //dedxpeakdata->GetYaxis()->SetNoExponent();
481  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
482  leg0->AddEntry(dedxpeakdata,"Data","LE");
483  leg0->AddEntry(dedxpeakmcscale,"MC","LE");
484  leg0->Draw("same");
485  c2->Print("dedxoptfine.pdf");
486 
487  double nonoptscale = 0.97375; //1.483/1.55;
488 
489  TH1D* dedxpeakmcnonoptscale = new TH1D("NonOptimalMatch","",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
490  dedxpeakmcnonoptscale->Add(dedxpeakmc);
491  dedxpeakmcnonoptscale->Scale(nonoptscale);
492 
493  TCanvas* c3 = new TCanvas("c3","c3");
494  dedxpeakdata->Draw();
495  dedxpeakmcnonoptscale->Draw("same");
496  dedxpeakmcnonoptscale->SetLineColor(kRed);
497  dedxpeakmcnonoptscale->SetMarkerColor(kRed);
498  dedxpeakdata->GetYaxis()->SetRangeUser(0.0014,0.0025);
499  dedxpeakdata->GetXaxis()->SetNoExponent();
500  //dedxpeakdata->GetYaxis()->SetNoExponent();
501  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
502  leg0->AddEntry(dedxpeakdata,"Data","LE");
503  leg0->AddEntry(dedxpeakmcnonoptscale,"MC","LE");
504  leg0->Draw("same");
505  c3->Print("dedxnonoptfine.pdf");
506 
507 
508  TCanvas* c4 = new TCanvas("c4","c4");
509  TH1D* dedxratio = MakeRatio(dedxpeakdata,dedxpeakmc,"dedxratio");
510  TLine* unitline = new TLine(-1,1,1001,1);
511  unitline->SetLineWidth(2);
512  unitline->SetLineStyle(7);
513  dedxratio->Draw();
514  unitline->Draw("same");
515  dedxratio->GetYaxis()->SetRangeUser(0.9,1.1);
516 
517  TCanvas* c5 = MakeSplitRatio(c0,c4,"c5");
518  c5->Print("dedxratio.pdf");
519 
520 
521 
522 
523 }
kRed
enum BeamMode kRed
MakeSplitRatio
TCanvas * MakeSplitRatio(TCanvas *top, TCanvas *bot, const char *name)
Definition: DedxGamma.C:39
nsteps
const int nsteps
Definition: plot_lowE_sig_spread_genie.C:35
chi2
double chi2()
Definition: cana_time_peak.C:532
dist
double dist
Definition: runWimpSim.h:113
demo5.c2
c2
Definition: demo5.py:33
make_syst_table_plots.ibin
ibin
Definition: make_syst_table_plots.py:246
PeakProjectionY
TH1D * PeakProjectionY(TH2F *hist, const char *name, double minrange=-1.0, double maxrange=-1.0, double lowsig=0.001, double highsig=0.001, bool useTruncate=true)
Definition: DedxGamma.C:265
demo5.expt
expt
Definition: demo5.py:34
release_diff.diff
diff
Definition: release_diff.py:80
om::cout
OStream cout
Definition: OStream.cxx:6
kViolet
enum BeamMode kViolet
MakeRatio
TH1D * MakeRatio(TH1D *num, TH1D *denom, const char *name)
Definition: DedxGamma.C:21
kGreen
enum BeamMode kGreen
void FindTruncate ( TH1D *  hist,
double &  low,
double &  high 
)

Definition at line 120 of file DedxGamma.C.

References make_syst_table_plots::ibin, and febshutoff_auto::integral.

Referenced by PeakProjectionY().

121 {
122  //std::cout<<"BEGINING OF TRUNCATE FUNCATION!"<<std::endl;
123  double percentCut = 0.05;
124  double histMin = hist->GetXaxis()->GetXmin();
125  double histMax = hist->GetXaxis()->GetXmax();
126  int minbin = hist->GetXaxis()->FindFixBin(histMin);
127  int maxbin = hist->GetXaxis()->FindFixBin(histMax)-1;
128  double integral = hist->Integral(minbin,maxbin);
129  //std::cout<<"Total integral is: "<<integral<<std::endl;
130  //if (integral == 363){
131  //TCanvas* cs = new TCanvas("cs","cs");
132  //hist->Draw();
133  //cs->Print("test.pdf");
134  //}
135  double lowCut = integral*percentCut;
136  double highCut = integral*(1-percentCut);
137 
138  int lowbin = 0;
139  int highbin = 0;
140  bool foundLow = false;
141  double ongoingIntegral = 0.0;
142 
143  for(int ibin = minbin; ibin <= maxbin; ++ibin){
144  double bincont = hist->GetBinContent(ibin);
145  //std::cout<<"bin: "<<ibin<<" with content: "<<bincont<<std::endl;
146  if (!foundLow){
147  ongoingIntegral += bincont;
148  //std::cout<<" ongoingIntegral: "<<ongoingIntegral<<std::endl;
149  if (ongoingIntegral > lowCut){
150  lowbin = ibin;
151  foundLow = true;
152  //std::cout<<"Found Low! Set lowbin : "<<lowbin<<std::endl;
153  }
154  }
155  else{
156  if ((ongoingIntegral<highCut)&&((ongoingIntegral+bincont)>=highCut)){
157  ongoingIntegral += bincont;
158  highbin = ibin;
159  //std::cout<<"Found High! Set highbin : "<<highbin<<std::endl;
160  }
161  else if (ongoingIntegral<highCut){
162  ongoingIntegral += bincont;
163  //std::cout<<" ongoingIntegral: "<<ongoingIntegral<<std::endl;
164  }
165  }//Already found low, now looking for high bin
166  //Now need to set low, high to exclude the lowbin and highbin
167  low = hist->GetXaxis()->GetBinUpEdge(lowbin);
168  high = hist->GetXaxis()->GetBinLowEdge(highbin);
169  }
170  return;
171 }
make_syst_table_plots.ibin
ibin
Definition: make_syst_table_plots.py:246
febshutoff_auto.integral
integral
Definition: febshutoff_auto.py:107
void GausFit ( TH1D *  hist,
double &  mean,
double &  err,
double  minrange,
double  maxrange,
double  lowsig,
double  highsig,
bool  useTruncate 
)

Definition at line 212 of file DedxGamma.C.

References om::cout, lem_server::cs, allTimeWatchdog::endl, febshutoff_auto::integral, submit_nova_art::mode, and ModeRestrict().

Referenced by PeakProjectionY().

213 {
214  std::cout<<" We are inside GausFit! "<<std::endl;
215  std::cout<<" we have minrange = "<<minrange<<" and maxrange = "<<maxrange<<std::endl;
216  std::cout<<" we have lowsig = "<<lowsig<<" and highsig = "<<highsig<<std::endl;
217  double mode = ModeRestrict(hist,minrange,maxrange);
218  std::cout<<" we have mode = "<<mode<<std::endl;
219  double fitminx = -1.0;
220  double fitmaxx = -1.0;
221  if (useTruncate){
222  fitminx = lowsig;
223  fitmaxx = highsig;
224  }
225  else{
226  fitminx = mode-lowsig;
227  fitmaxx = mode+highsig;
228  }
229  //std::cout<<" we have fitminx = "<<fitminx<<" and fitmaxx = "<<fitmaxx<<std::endl;
230  if(fitminx < minrange){ fitminx = minrange; }
231  if(fitmaxx > maxrange){ fitmaxx = maxrange; }
232  //std::cout<<" we have fitminx = "<<fitminx<<" and fitmaxx = "<<fitmaxx<<std::endl;
233  TF1* fgaus = new TF1("gaus","gaus(0)",fitminx,fitmaxx);
234  std::cout<<"fitminx: "<<fitminx<<" fitmaxx: "<<fitmaxx<<" mode: "<<mode<<std::endl;
235  fgaus->SetParameters(hist->GetMaximum(),mode,(fitmaxx-fitminx)/2.0);
236  int minbin = hist->GetXaxis()->FindFixBin(fitminx);
237  int maxbin = hist->GetXaxis()->FindFixBin(fitmaxx);
238  double integral = hist->Integral(minbin,maxbin);
239  mean = -1;
240  err = 0.000001;
241  if(integral > 100){
242  std::cout<<" integral is "<<integral<<std::endl;
243  if (integral == 320){
244  TCanvas* cs = new TCanvas("cs","cs");
245  hist->Fit(fgaus,"RIQ");
246  cs->Print("test.pdf");
247  }
248  int fitstatus = hist->Fit(fgaus,"RI");
249  std::cout<<"And our fit status is: "<<fitstatus<<std::endl;
250  if(fitstatus == 0){
251  mean = fgaus->GetParameter(1);
252  err = fgaus->GetParError(1);
253  if(mean > maxrange || mean < minrange){
254  mean = -1;
255  err = 0.000001;
256  }
257  //err = fgaus->GetParameter(2);
258  }
259  }
260 
261  std::cout<<"mean: "<<mean<<" err: "<<err<<std::endl;
262 
263 }
ModeRestrict
double ModeRestrict(TH1D *hist, double minx, double maxx)
Definition: DedxGamma.C:182
submit_nova_art.mode
mode
Definition: submit_nova_art.py:715
febshutoff_auto.integral
integral
Definition: febshutoff_auto.py:107
om::cout
OStream cout
Definition: OStream.cxx:6
lem_server.cs
cs
Definition: lem_server.py:582
TH1D* MakeRatio ( TH1D *  num,
TH1D *  denom,
const char *  name 
)

Definition at line 21 of file DedxGamma.C.

References PandAna.reco_validation.prod5_pid_validation::ratio().

Referenced by DedxGamma().

22 {
23 
24  TH1D* ratio = new TH1D(name,";;Data/MC",num->GetNbinsX(),num->GetXaxis()->GetXbins()->GetArray());
25 
26  ratio->GetXaxis()->SetTitle(num->GetXaxis()->GetTitle());
27 
28  ratio->Divide(num,denom);
29 
30  ratio->SetLineColor(denom->GetLineColor());
31  ratio->SetMarkerColor(denom->GetMarkerColor());
32  ratio->SetMarkerSize(denom->GetMarkerSize());
33 
34  return ratio;
35 
36 }
name
const XML_Char * name
Definition: expat.h:151
PandAna.reco_validation.prod5_pid_validation.ratio
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
Definition: prod5_pid_validation.py:554
num
int num
Definition: f2_nu.C:119
TCanvas* MakeSplitRatio ( TCanvas *  top,
TCanvas *  bot,
const char *  name 
)

Definition at line 39 of file DedxGamma.C.

References allTimeWatchdog::can, and stan::math::fabs().

Referenced by DedxGamma().

40 {
41 
42  TCanvas* can = new TCanvas(name,name,600,550);
43 
44  double pixels = 18.0;
45  double labPix = 15.0;
46  double labOff = 5.0;
47 
48  can->SetBottomMargin(0);
49  can->Divide(1,2,0,0);
50  TVirtualPad* can_1 = can->cd(1);
51  can_1->SetPad(0,0.4,1,1);
52  can_1->SetFillStyle(0);
53  top->DrawClonePad();
54  can_1->SetRightMargin(0.05);
55  can_1->SetTopMargin(0.1);
56  can_1->SetBottomMargin(0.005);
57  can_1->SetLeftMargin(0.13);
58  can_1->Update();
59  double textSizeFac = 1.0/(fabs(can_1->YtoAbsPixel(can_1->GetY2()) - can_1->YtoAbsPixel(can_1->GetY1())));
60  double marginSize = 80;
61  double textsizeA = pixels/(can_1->GetAbsHNDC()*can_1->GetWh());
62  double textsizeAl = labPix/(can_1->GetAbsHNDC()*can_1->GetWh());
63  double textsizeAo = labOff/(can_1->GetWw() * can_1->GetAbsWNDC());
64 
65  TList* primlist = can_1->GetListOfPrimitives();
66  for(int ilist = 0; ilist < primlist->GetSize(); ++ilist){
67  TObject* prim = primlist->At(ilist);
68  if(prim->InheritsFrom("TH1")){
69  TH1* primhist = (TH1*)prim;
70  primhist->GetYaxis()->SetTitleSize(textsizeA);
71  primhist->GetYaxis()->SetLabelSize(textsizeAl);
72  primhist->GetYaxis()->SetLabelOffset(textsizeAo);
73  primhist->GetYaxis()->SetTitleOffset(0.035 / textsizeA);
74  primhist->GetXaxis()->SetNoExponent(kTRUE);
75  }
76  }
77 
78  can_1->Update();
79 
80  TVirtualPad* can_2 = can->cd(2);
81  can_2->SetPad(0,0,1,0.4);
82  bot->DrawClonePad();
83  can_2->SetTopMargin(0.025);
84  can_2->SetBottomMargin(0.22);
85  can_2->SetRightMargin(0.05);
86  can_2->SetLeftMargin(0.13);
87  can_2->Update();
88 
89  double textsize = pixels/(can_2->GetAbsHNDC()*can_2->GetWh());
90  double textsizel = labPix/(can_2->GetAbsHNDC()*can_2->GetWh());
91  double textsizeo = labOff/(can_2->GetWw() * can_2->GetAbsWNDC());
92  double textsizex = pixels/(can_2->GetAbsHNDC()*can_2->GetWh());
93  double textsizelx = labPix/(can_2->GetAbsHNDC()*can_2->GetWh());
94  textsizeAl = labPix/(can_2->GetAbsHNDC()*can_2->GetWh());
95  double textsizeox = labOff/(can_2->GetWw() * can_2->GetAbsWNDC());
96 
97  TList* primlist_2 = can_2->GetListOfPrimitives();
98  for(int ilist = 0; ilist < primlist_2->GetSize(); ++ilist){
99  TObject* prim = primlist_2->At(ilist);
100  if(prim->InheritsFrom("TH1")){
101  TH1* primhist = (TH1*)prim;
102  primhist->GetYaxis()->SetTitleSize(textsize);
103  primhist->GetYaxis()->SetLabelSize(textsizel);
104  primhist->GetYaxis()->SetLabelOffset(textsizeAo);
105  primhist->GetXaxis()->SetTitleSize(textsizex);
106  primhist->GetXaxis()->SetLabelSize(textsizelx);
107  primhist->GetXaxis()->SetLabelOffset(textsizeox);
108  primhist->GetYaxis()->SetTitleOffset(0.035 / textsize);
109  primhist->GetXaxis()->SetTitleOffset(0.9);
110  primhist->GetXaxis()->SetNoExponent(kTRUE);
111  }
112  }
113 
114  can_2->Update();
115 
116  return can;
117 
118 }
name
const XML_Char * name
Definition: expat.h:151
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
top
TGeoVolume * top
Definition: make_fe_box.C:9
double Mode ( TH1D *  hist)

Definition at line 173 of file DedxGamma.C.

Referenced by ModeRestrict().

174 {
175 
176  double modex = hist->GetXaxis()->GetBinCenter(hist->GetMaximumBin());
177  //std::cout<<"In MODE with Modex: "<<modex<<std::endl;
178  //SL -- adding this line, because I think it needs something ....
179  return modex;
180 }
double ModeRestrict ( TH1D *  hist,
double  minx,
double  maxx 
)

Definition at line 182 of file DedxGamma.C.

References make_syst_table_plots::ibin, Mode(), and submit_nova_art::mode.

Referenced by GausFit().

183 {
184 
185  double mode = Mode(hist);
186  //std::cout<<"In mode restrict with Mode: "<<mode<<" maxx: "<<maxx<<" minx: "<<minx<<std::endl;
187  if( mode > maxx || mode < minx){
188  // loop over bins between maxx and minx to find the mode in the restricted range;
189  //std::cout<<"HERE HERE HERE!!!"<<std::endl;
190  int minxbin = hist->GetXaxis()->FindFixBin(minx);
191  int maxxbin = hist->GetXaxis()->FindFixBin(maxx);
192  double modeval = hist->GetBinContent(minxbin);
193  int modebin = minxbin;
194 // std::cout<<"minxbin: "<<minxbin<<" maxxbin: "<<maxxbin<<" modeval: "<<modeval<<std::endl;
195  for(int ibin = minxbin+1; ibin <= maxxbin; ++ibin){
196  double bincont = hist->GetBinContent(ibin);
197 // std::cout<<"ibin: "<<ibin<<" bincont: "<<bincont<<std::endl;
198  if(bincont > modeval){
199  modebin = ibin;
200  modeval = bincont;
201  //std::cout<<"modeval: "<<modeval<<" modebin: "<<modebin<<std::endl;
202  }
203  }
204  mode = hist->GetXaxis()->GetBinCenter(modebin);
205  }
206 
207  return mode;
208 
209 }
submit_nova_art.mode
mode
Definition: submit_nova_art.py:715
make_syst_table_plots.ibin
ibin
Definition: make_syst_table_plots.py:246
Mode
double Mode(TH1D *hist)
Definition: DedxGamma.C:173
TH1D* PeakProjectionY ( TH2F *  hist,
const char *  name,
double  minrange = -1.0,
double  maxrange = -1.0,
double  lowsig = 0.001,
double  highsig = 0.001,
bool  useTruncate = true 
)

Definition at line 265 of file DedxGamma.C.

References om::cout, allTimeWatchdog::endl, FindTruncate(), GausFit(), and extractScale::mean.

Referenced by DedxGamma().

266 {
267  std::cout<<"In peak projection!"<<std::endl;
268  TH1D* peakhist = new TH1D(name,"",hist->GetNbinsX(),hist->GetXaxis()->GetXbins()->GetArray());
269  peakhist->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
270  peakhist->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
271  //std::cout<<"And the axis title is: "<<hist->GetXaxis()->GetTitle()<<std::endl;
272  double fitmin = minrange;
273  double fitmax = maxrange;
274  //std::cout<<"Minrange: "<<fitmin<<" and maxrange: "<<fitmax<<std::endl;
275  if(minrange == -1.0){ fitmin = hist->GetXaxis()->GetXmin(); }
276  if(maxrange == -1.0){ fitmax = hist->GetXaxis()->GetXmax(); }
277  //std::cout<<"Minrange: "<<fitmin<<" and maxrange: "<<fitmax<<std::endl;
278  //std::cout<<"Number of x bins: "<<hist->GetNbinsX()<<std::endl;
279  for(int ixbin = 1; ixbin <= hist->GetNbinsX(); ++ixbin){
280  // get the projection of the 2d hist
281  TH1D* projection = hist->ProjectionY("",ixbin,ixbin,"o");
282  //std::cout<<" X bin: "<<ixbin<<" Entries: "<<projection->GetEntries()<<std::endl;
283  if(projection->GetEntries() == 0){ continue; }
284  // fit the projection to a Gaussian around the peak
285  double mean = 0.00175 ,err = 0.001;
286  double lowTruncate = -1.0;
287  double highTruncate = -1.0;
288  if (useTruncate){
289  FindTruncate(projection, lowTruncate, highTruncate);
290  if ((lowTruncate != -1.0) && (highTruncate != -1.0)){
291  lowsig = lowTruncate;
292  highsig = highTruncate;
293  }
294  }
295  //std::cout<<"Low truncate: "<<lowTruncate<<" High truncate: "<<highTruncate<<std::endl;
296  GausFit(projection,mean,err,fitmin,fitmax,lowsig,highsig,useTruncate);
297  //std::cout<<" Mean: "<<mean<<" Error: "<<err<<std::endl;
298  peakhist->SetBinContent(ixbin,mean);
299  peakhist->SetBinError(ixbin,err);
300 
301  }
302 
303  return peakhist;
304 
305 }
name
const XML_Char * name
Definition: expat.h:151
GausFit
void GausFit(TH1D *hist, double &mean, double &err, double minrange, double maxrange, double lowsig, double highsig, bool useTruncate)
Definition: DedxGamma.C:212
FindTruncate
void FindTruncate(TH1D *hist, double &low, double &high)
Definition: DedxGamma.C:120
om::cout
OStream cout
Definition: OStream.cxx:6
Generated on Sun Jan 24 2021 14:15:37 for NOvA by   doxygen 1.8.11