ATLAS Offline Software
PixelChargeInterpolationHistograms.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #ifndef PixelChargeInterpolationHistograms_C
6 #define PixelChargeInterpolationHistograms_C
7 
8 #include <iostream>
9 #include <fstream>
10 #include <cmath>
11 #include <sstream>
12 #include <string>
13 #include <vector>
14 #include <sys/types.h>
15 #include <sys/stat.h>
16 
17 #include <TProfile.h>
18 #include <TH2.h>
19 #include <TH1.h>
20 #include <TF1.h>
21 
27 
28 namespace PixelCalib{
29 
32  m_tag(tag),
33  m_parameters(0),
34  m_etaProfile(0),
35  m_phiProfile(0),
36  m_etaH(0),
37  m_phiH(0),
38  m_OmegaPhih(0),
39  m_OmegaEtah(0){
40 
41  std::vector<float> *layers = getLayersBins();
42  std::vector<float> etabins = model.getEtaBins();
43  std::vector<float> phibins = model.getAngleBins();
44  std::vector<float> clustersizeEta = model.getClusterSizeYBins();
45  std::vector<float> clustersizePhi = model.getClusterSizeXBins();
46 
47  unsigned int Neta = etabins.size()-1;
48  unsigned int Nphi = phibins.size()-1;
49  unsigned int NCSeta = clustersizeEta.size()-1;
50  unsigned int NCSphi = clustersizePhi.size()-1;
51 
52  unsigned int ntotbins = Nphi + NCSphi + Neta + NCSeta + 4;
53 
54  std::vector<float> bins;
55  bins.reserve(ntotbins);
56  for (unsigned int i=0; i<clustersizePhi.size(); i++) {
57  bins.push_back(clustersizePhi[i]);
58  }
59  for (unsigned int i=0; i<clustersizeEta.size(); i++) {
60  bins.push_back(clustersizeEta[i]);
61  }
62  for (unsigned int i=0; i<etabins.size(); i++) {
63  bins.push_back(etabins[i]);
64  }
65  for (unsigned int i=0; i<phibins.size(); i++) {
66  bins.push_back(phibins[i]);
67  }
68 
70  m_parameters->setParameters(NCSphi, NCSeta, Neta, Nphi, 0, bins);
72 
73  // eta direction
74  TProfile *Profmodel = new TProfile(("etaResVsOmega"+ tag).c_str(),
75  "#eta residual vs charge sharing", 30, 0., 1.);
76 
77  TH2F *THmodel = new TH2F(("etaResVsOmegaH"+ tag).c_str(),
78  "#eta residual vs charge sharing",
79  50, 0., 1.,50,-1000,1000);
80 
81  std::vector<std::string> binsnames(3);
82  std::vector<std::vector <float> > binsvectors(3);
83 
84  binsnames[LayerIndex] = "Layer";
85  binsnames[AngleIndex] = "#eta";
86  binsnames[ClustersizeIndex] = "ClusterSize";
87 
88  binsvectors[LayerIndex] = *layers;
89  binsvectors[AngleIndex] = etabins;
90  binsvectors[ClustersizeIndex] = clustersizeEta;
91 
92  m_etaProfile = new MultiHisto<TProfile>(*Profmodel,binsnames,binsvectors);
93  m_etaH = new MultiHisto<TH2F>(*THmodel,binsnames,binsvectors);
94  m_OmegaEtah_model = new TH1F(("m_OmegaEtah"+ tag).c_str(), "#Omega_{y} distribution",125,0,1);
95  m_OmegaEtah = new MultiHisto<TH1F>(*m_OmegaEtah_model,binsnames,binsvectors);
96 
97  // phi direction
98  Profmodel->SetName(("phiResVsOmega"+ tag).c_str());
99  Profmodel->SetTitle("#phi residual vs charge sharing");
100 
101  TH2F *THmodel1 = new TH2F(("phiResVsOmegaH"+ tag).c_str(),
102  "#phi residual vs charge balancing",
103  50, 0., 1.,50,-200,200);
104 
105  binsnames[AngleIndex] = "#phi";
106 
107  binsvectors[AngleIndex] = phibins;
108  binsvectors[ClustersizeIndex] = clustersizePhi;
109 
110  m_phiProfile = new MultiHisto<TProfile>(*Profmodel,binsnames,binsvectors);
111  m_phiH = new MultiHisto<TH2F>(*THmodel1,binsnames,binsvectors);
112 
113  m_OmegaPhih_model = new TH1F(("m_OmegaPhih"+ tag).c_str(), "#Omega_{x} distribution",125,0,1);
114  m_OmegaPhih = new MultiHisto<TH1F>(*m_OmegaPhih_model,binsnames,binsvectors);
115 
116  delete Profmodel;
117  delete THmodel;
118  delete THmodel1;
119  delete layers;
120  Profmodel = 0;
121  THmodel = 0;
122  THmodel1 = 0;
123 
124 }
125 
127 
129 
130  delete m_etaProfile;
131  delete m_phiProfile;
132  delete m_etaH;
133  delete m_phiH;
134  delete m_parameters;
135  m_etaProfile = 0;
136  m_phiProfile = 0;
137  m_etaH = 0;
138  m_phiH = 0;
139  m_parameters = 0;
140 
141  delete m_OmegaPhih; m_OmegaPhih = 0;
142  delete m_OmegaEtah; m_OmegaEtah = 0;
145 }
146 
147 
149 
151  double TrkEta, double DeltaCol, double reseta, double OmegaEta,
152  double alpha, double DeltaRow, double resphi, double OmegaPhi){
153 
154 
155  std::vector<double> Pars(3);
156  if(GeVTrkPt == 0) return -1;
157 
158  if( DeltaCol > 1){ // otherwise none to share with!
159 
160  Pars[AngleIndex] = TrkEta;
161  Pars[LayerIndex] = DetType;
162  //std::cout << DetType << std::endl;
163  Pars[ClustersizeIndex] = DeltaCol;
164  m_OmegaEtah->Fill(OmegaEta,1,Pars);
165  m_OmegaEtah_model->Fill(OmegaEta);
166 
167  if(OmegaEta > 0.1 && OmegaEta < 0.9){
168  m_etaProfile->Fill(OmegaEta,reseta,Pars);
169  m_etaH->Fill(OmegaEta,reseta,Pars);
170  }
171  }
172 
173  if( DeltaRow > 1){// otherwise none to share with!
174 
175  Pars[AngleIndex] = alpha;
176  Pars[LayerIndex] = DetType;
177  Pars[ClustersizeIndex] = DeltaRow;
178 
179  m_OmegaPhih->Fill(OmegaPhi,1,Pars);
180  m_OmegaPhih_model->Fill(OmegaPhi);
181 
182  if(OmegaPhi > 0.1 && OmegaPhi < 0.9){
183  m_phiProfile->Fill(OmegaPhi,resphi,Pars);
184  m_phiH->Fill(OmegaPhi,resphi,Pars);
185  }
186 
187  }
188 
189 
190 
191  return 0;
192 
193 }
194 
196 
198 
199  int readhistos = 0;
200  TDirectory *histodir = (TDirectory *)gDirectory->Get(m_etaProfile->GetName());
201  readhistos += m_etaProfile->FillFromFile(histodir);
202  histodir = (TDirectory *)gDirectory->Get(m_etaH->GetName());
203  readhistos += m_etaH->FillFromFile(histodir);
204  histodir = (TDirectory *)gDirectory->Get(m_phiProfile->GetName());
205  readhistos += m_phiProfile->FillFromFile(histodir);
206  histodir = (TDirectory *)gDirectory->Get(m_phiH->GetName());
207  readhistos += m_phiH->FillFromFile(histodir);
208 
209  return readhistos;
210 }
211 
213 
215 
216  m_etaProfile->Write();
217  m_phiProfile->Write();
218  m_etaH->Write();
219  m_phiH->Write();
220 
221  m_OmegaPhih->Write();
222  m_OmegaEtah->Write();
223  m_OmegaPhih_model->Write();
224  m_OmegaEtah_model->Write();
225 
228 }
229 
231 
233  std::ofstream &logfile){
234 
235  logfile << "Fitting!" << std::endl;
236 
237  SetAtlasStyle();
238  TCanvas *c1 = new TCanvas();
239  c1->UseCurrentStyle();
240 
241  //char *currpath = getcwd(nullptr,0);
242  //mkdir(m_etaProfile->GetName(),S_IRWXU | S_IRWXG | S_IRWXO);
243  //chdir(m_etaProfile->GetName());
244 
245  for(unsigned int i = 0; i < m_etaProfile->GetNhistos(); i++){
247  swap->UseCurrentStyle();
248  double value = 0, error = 0;
249  if(Fit(swap, &value, &error) ){
250  logfile << swap->GetTitle() << " --> "
251  << value << " +/- " << error << std::endl;
252  }else logfile << swap->GetTitle()
253  << " --> Failing fit!" << std::endl;
254  std::vector<int> indexes = m_etaProfile->GetDivisionsIndexes(i);
256  indexes[ClustersizeIndex],
257  indexes[LayerIndex],
258  value/1000);
260  indexes[ClustersizeIndex],
261  indexes[LayerIndex],
262  error/1000);
263  if(swap->GetEntries() < 100) continue;
264  std::string name = std::string(swap->GetName()) + std::string(".pdf");
265  //c1->Print(name.c_str());
266  TH2F *swap1 = m_etaH->GetHisto(i);
267  swap1->UseCurrentStyle();
268  swap1->SetMarkerSize(0.2);
269  swap1->GetXaxis()->SetTitle("Charge sharing");
270  swap1->GetYaxis()->SetTitle("Cluster center residuals (#mum)");
271  swap1->GetYaxis()->SetTitleOffset(1.2);
272  swap1->GetXaxis()->SetTitleOffset(1.25);
273  swap1->DrawCopy();
274  DrawTitleLatex(swap1->GetTitle(), 0.6,0.85);
275  std::string name1 = std::string(swap1->GetName()) + std::string(".pdf");
276  //c1->Print(name1.c_str());
277  }
278 
279  //chdir(currpath);
280  //mkdir(m_phiProfile->GetName(),S_IRWXU | S_IRWXG | S_IRWXO);
281  //chdir(m_phiProfile->GetName());
282  for(unsigned int i = 0; i < m_phiProfile->GetNhistos(); i++){
284  swap->UseCurrentStyle();
285  double value = 0, error = 0;
286  if(Fit(swap, &value, &error) ){
287  logfile << swap->GetTitle() << " --> "
288  << value << " +/- " << error << std::endl;
289  }else logfile << swap->GetTitle()
290  << " --> Failing fit!" << std::endl;
291  std::vector<int> indexes = m_phiProfile->GetDivisionsIndexes(i);
293  indexes[ClustersizeIndex],
294  indexes[LayerIndex],
295  value/1000);
297  indexes[ClustersizeIndex],
298  indexes[LayerIndex],
299  error/1000);
300  if(swap->GetEntries() < 100) continue;
301  std::string name = std::string(swap->GetName()) + std::string(".pdf");
302  //c1->Print(name.c_str());
303  TH2 *swap1 = m_phiH->GetHisto(i);
304  swap1->UseCurrentStyle();
305  swap1->SetMarkerSize(0.2);
306  swap1->GetXaxis()->SetTitle("Charge sharing");
307  swap1->GetYaxis()->SetTitle("Residuals from center of the cluster (#mum)");
308  swap1->GetYaxis()->SetTitleOffset(1.2);
309  swap1->GetXaxis()->SetTitleOffset(1.25);
310  swap1->DrawCopy();
311  DrawTitleLatex(swap1->GetTitle(), 0.6,0.85);
312  std::string name1 = std::string(swap1->GetName()) + std::string(".pdf");
313  //c1->Print(name1.c_str());
314  }
315 
316  //chdir(currpath);
317  //delete currpath;
318  delete c1;
319  c1 = 0;
320  //currpath = 0;
321 
322  return m_parameters;
323 }
324 
326 
328 
329  // perform fits
330  if(swap->GetEntries() > 100){
331  TF1 *fitfunc = new TF1("fitfunc","[1] * x + [0]",0.15,0.85);
332  fitfunc->SetParameter(0,0);
333  fitfunc->SetParameter(1,0);
334  //fitfunc->SetParLimits(1,lowlim,hilim);
335  if(swap->Fit("fitfunc","QR") == 0 && fitfunc->GetProb() > 0.005){
336  *value = - fitfunc->GetParameter(1);
337  *error = fitfunc->GetParError(1);
338  if(fabs(*value) < 2*(*error) || *value < 0){
339  *value = 0;
340  *error = 0;
341  }
342  swap->GetXaxis()->SetTitle("Charge sharing");
343  swap->GetYaxis()->SetTitle("Residuals from the center of the cluster (#mum)");
344  swap->GetYaxis()->SetTitleOffset(1.2);
345  swap->GetXaxis()->SetTitleOffset(1.25);
346  swap->DrawCopy();
347 
348  std::ostringstream FitString;
349  FitString.flags(std::ios::fixed);
350  FitString.precision(2);
351  FitString << "slope: " << *value << " #pm " << *error << " #mum";
352 
353  DrawATLASLabel(0.2,0.4);
354  //DrawTitleLatex(FitString.str().c_str(), 0.2,0.3);
355  DrawTitleLatex(swap->GetTitle(), 0.6,0.85);
356 
357  delete fitfunc;
358  fitfunc = 0;
359  return true;
360  }else{
361  std::cout << swap->GetTitle() << " " << swap->Fit("fitfunc","QR") << " " << fitfunc->GetProb() << std::endl;
362  *value = 0;
363  *error = 0;
364  return false;
365  }
366 
367  }else return false;
368 }
369 
370 }// end of namespace PixelCalib
371 
372 #endif // #ifdef PixelChargeInterpolationHistograms_C
PixelCalib::PixelChargeInterpolationHistograms::m_etaH
MultiHisto< TH2F > * m_etaH
Definition: PixelChargeInterpolationHistograms.h:47
MultiHisto::Write
int Write(const char *name=0, Int_t option=0, Int_t bufsize=0) const
SetAtlasStyle
void SetAtlasStyle()
Definition: InnerDetector/InDetCalibAlgs/PixelCalibAlgs/Macro/AtlasStyle.h:17
PixelCalib::PixelChargeInterpolationParameters::setDeltaY
int setDeltaY(int ieta, int iclustersize, int ilayer, float value)
Definition: PixelChargeInterpolationParameters.cxx:277
ReadCellNoiseFromCool.name1
name1
Definition: ReadCellNoiseFromCool.py:233
MultiHisto::Fill
int Fill(double xvar, double yvar, const std::vector< double > &pars)
python.App.bins
bins
Definition: App.py:410
PixelCalib::PixelChargeInterpolationHistograms::m_phiH
MultiHisto< TH2F > * m_phiH
Definition: PixelChargeInterpolationHistograms.h:48
TH2F
Definition: rootspy.cxx:420
PixelCalib::PixelChargeInterpolationHistograms::m_OmegaPhih_model
TH1F * m_OmegaPhih_model
Definition: PixelChargeInterpolationHistograms.h:52
MultiHisto::GetHisto
ht * GetHisto(int globalindex)
extractSporadic.c1
c1
Definition: extractSporadic.py:134
module_driven_slicing.layers
layers
Definition: module_driven_slicing.py:114
MultiHisto::GetNhistos
unsigned int GetNhistos() const
MultiHisto::GetDivisionsIndexes
std::vector< int > GetDivisionsIndexes(int globalindex) const
PixelCalib::PixelChargeInterpolationHistograms::Fill
int Fill(int DetType, double GeVTrkPt, double TrkEta, double DeltaCol, double reseta, double OmegaEta, double alpha, double DeltaRow, double resphi, double OmegaPhi)
Definition: PixelChargeInterpolationHistograms.cxx:150
athena.value
value
Definition: athena.py:122
MultiHisto< TProfile >
PixelCalib::PixelChargeInterpolationHistograms::m_OmegaPhih
MultiHisto< TH1F > * m_OmegaPhih
Definition: PixelChargeInterpolationHistograms.h:49
AtlasStyle.h
PixelCalib::PixelChargeInterpolationHistograms::AngleIndex
@ AngleIndex
Definition: PixelChargeInterpolationHistograms.h:56
python.TrigEgammaMonitorHelper.TH2F
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:45
PixelCalib::PixelChargeInterpolationHistograms::Write
int Write()
Definition: PixelChargeInterpolationHistograms.cxx:214
PixelCalib::PixelChargeInterpolationHistograms::Analyze
PixelChargeInterpolationParameters * Analyze(std::ofstream &logfile)
Definition: PixelChargeInterpolationHistograms.cxx:232
PixelCalib::PixelChargeInterpolationHistograms::m_OmegaEtah
MultiHisto< TH1F > * m_OmegaEtah
Definition: PixelChargeInterpolationHistograms.h:50
PixelChargeInterpolationHistograms.h
OfflineCalibUtils.icc
ParseInputs.gDirectory
gDirectory
Definition: Final2012/ParseInputs.py:133
lumiFormat.i
int i
Definition: lumiFormat.py:92
python.TrigEgammaMonitorHelper.TProfile
def TProfile(*args, **kwargs)
Definition: TrigEgammaMonitorHelper.py:81
DrawTitleLatex
void DrawTitleLatex(const char *chartitle, float x, float y, int color=1, float textsize=0.04)
Definition: InnerDetector/InDetCalibAlgs/PixelCalibAlgs/Macro/AtlasStyle.h:97
PixelCalib
Definition: PixelChargeInterpolationCalibration.h:14
PixelCalib::PixelChargeInterpolationHistograms::m_etaProfile
MultiHisto< TProfile > * m_etaProfile
Definition: PixelChargeInterpolationHistograms.h:45
WriteCalibToCool.swap
swap
Definition: WriteCalibToCool.py:94
TH2
Definition: rootspy.cxx:373
PixelChargeInterpolationParameters.h
PixelCalib::PixelChargeInterpolationParameters::setVersion
void setVersion(int version)
Definition: PixelChargeInterpolationParameters.cxx:86
PixelCalib::PixelChargeInterpolationHistograms::ClustersizeIndex
@ ClustersizeIndex
Definition: PixelChargeInterpolationHistograms.h:57
PixelCalib::PixelChargeInterpolationHistograms::m_OmegaEtah_model
TH1F * m_OmegaEtah_model
Definition: PixelChargeInterpolationHistograms.h:51
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
PixelCalib::PixelChargeInterpolationHistograms::PixelChargeInterpolationHistograms
PixelChargeInterpolationHistograms(const std::string &tag, const PixelChargeInterpolationParameters &model)
Definition: PixelChargeInterpolationHistograms.cxx:30
python.runDiffRootOnChanged.logfile
logfile
Definition: runDiffRootOnChanged.py:143
TProfile
Definition: rootspy.cxx:515
DrawATLASLabel
void DrawATLASLabel(float x, float y, bool pre=false, float textsize=0.05)
Definition: InnerDetector/InDetCalibAlgs/PixelCalibAlgs/Macro/AtlasStyle.h:85
GlobalVariables.phibins
int phibins
Definition: GlobalVariables.py:59
PixelCalib::PixelChargeInterpolationHistograms::m_parameters
PixelChargeInterpolationParameters * m_parameters
Definition: PixelChargeInterpolationHistograms.h:44
PixelCalib::PixelChargeInterpolationParameters::setErrDeltaY
int setErrDeltaY(int ieta, int iclustersize, int ilayer, float value)
Definition: PixelChargeInterpolationParameters.cxx:310
correlationModel::model
model
Definition: AsgElectronEfficiencyCorrectionTool.cxx:46
PixelCalib::PixelChargeInterpolationHistograms::m_phiProfile
MultiHisto< TProfile > * m_phiProfile
Definition: PixelChargeInterpolationHistograms.h:46
PixelCalib::PixelChargeInterpolationParameters::setDeltaX
int setDeltaX(int iangle, int iclustersize, int ilayer, float value)
Definition: PixelChargeInterpolationParameters.cxx:269
PixelCalib::PixelChargeInterpolationParameters::setParameters
void setParameters(const int ncsx, const int ncsy, const int neta, const int nalpha, int offset, std::vector< float > constants)
Definition: PixelChargeInterpolationParameters.cxx:90
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
PixelCalib::PixelChargeInterpolationHistograms::Fit
bool Fit(TProfile *swap, double *value, double *error)
Definition: PixelChargeInterpolationHistograms.cxx:327
PixelCalib::PixelChargeInterpolationParameters
Definition: PixelChargeInterpolationParameters.h:26
PixelCalib::PixelChargeInterpolationHistograms::LayerIndex
@ LayerIndex
Definition: PixelChargeInterpolationHistograms.h:55
MultiHisto::FillFromFile
int FillFromFile(TDirectory *histodir=0)
CaloCondBlobAlgs_fillNoiseFromASCII.tag
string tag
Definition: CaloCondBlobAlgs_fillNoiseFromASCII.py:24
MultiHisto.icc
DetType
Definition: DetType.h:10
PixelCalib::PixelChargeInterpolationHistograms::Read
int Read()
Definition: PixelChargeInterpolationHistograms.cxx:197
MultiHisto::GetName
virtual const char * GetName() const
error
Definition: IImpactPoint3dEstimator.h:70
PixelCalib::PixelChargeInterpolationHistograms::~PixelChargeInterpolationHistograms
virtual ~PixelChargeInterpolationHistograms()
Definition: PixelChargeInterpolationHistograms.cxx:128
PixelCalib::PixelChargeInterpolationParameters::setErrDeltaX
int setErrDeltaX(int iangle, int iclustersize, int ilayer, float value)
Definition: PixelChargeInterpolationParameters.cxx:302