ATLAS Offline Software
Loading...
Searching...
No Matches
ZmumuValidationExample.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
3*/
4
5#define ZMUMUVALIDATIONEXAMPLE_C
6
7#include <iostream>
8#include <list>
9#include "TMath.h"
10#include "TLorentzVector.h"
11#include "TF1.h"
12
13#include "ZmumuValidationExample.h"
14#include "TruthUtils/ParticleConstants.h"
15
16ZmumuValidationExample::ZmumuValidationExample( std::list<std::string> const & s_fileNames, string s_treeName
17 , std::string const & s_outFileName
18 , bool isMC
19 ) : m_fileNames( s_fileNames )
20 , m_outfilename( s_outFileName )
21 , m_file( s_outFileName.c_str(), "RECREATE" )
22 , m_eventChain( s_treeName.c_str() )
23 , m_truthChain( "TruthParams" )
24 , m_PrintLevel(0)
25 , m_EtaBins(1)
26 , m_PhiBins(1)
27 , m_isMC(isMC)
28 , n_iteration(0)
29
30{
31 loadChains();
32 setBranchAddresses();
33 // bookHistograms();
34
35 m_factor = 0.5; //factor to assign lambda corrections - nominal is 0.5
36 m_shift = 0; //shift in Z mass in MeV - nominal is 0
37 TempCanvasIterGaussFit = NULL;
38}
39
40
46
47
49// load TChains using the user inputed list of strings
51void ZmumuValidationExample::loadChains()
52{
53 std::cout << "Input Files:" << std::endl;
54
55 std::list<std::string>::const_iterator theEnd = m_fileNames.end();
56 for( std::list<std::string>::const_iterator itr = m_fileNames.begin()
57 ; itr != theEnd
58 ; ++itr
59 )
60 {
61 std::cout <<" "<< *itr << std::endl;
62
63 m_eventChain.Add( itr->c_str() );
64 if (m_isMC) m_truthChain.Add( itr->c_str() );
65 }
66
67 return;
68}
69
70
72// sets the branch address for each leaf to be used
74void ZmumuValidationExample::setBranchAddresses()
75{
77 // set branch addresses for event tree
79 m_eventChain.SetBranchAddress( "Positive_Px", &m_px_pos, &b_px_pos);
80 m_eventChain.SetBranchAddress( "Positive_Py", &m_py_pos, &b_py_pos);
81 m_eventChain.SetBranchAddress( "Positive_Pz", &m_pz_pos, &b_pz_pos);
82
83 m_eventChain.SetBranchAddress( "Negative_Px", &m_px_neg, &b_px_neg);
84 m_eventChain.SetBranchAddress( "Negative_Py", &m_py_neg, &b_py_neg);
85 m_eventChain.SetBranchAddress( "Negative_Pz", &m_pz_neg, &b_pz_neg);
86
87 m_eventChain.SetBranchAddress( "Positive_z0", &m_z0_pos, &b_z0_pos);
88 m_eventChain.SetBranchAddress( "Positive_d0", &m_d0_pos, &b_d0_pos);
89
90 m_eventChain.SetBranchAddress( "Negative_z0", &m_z0_neg, &b_z0_neg);
91 m_eventChain.SetBranchAddress( "Negative_d0", &m_d0_neg, &b_d0_neg);
92
93 if (m_isMC) {
94
95 m_truthChain.SetBranchAddress("Positive_Px", &m_truth_px_pos, &b_truth_px_pos);
96 m_truthChain.SetBranchAddress("Positive_Py", &m_truth_py_pos, &b_truth_py_pos);
97 m_truthChain.SetBranchAddress("Positive_Pz", &m_truth_pz_pos, &b_truth_pz_pos);
98
99 m_truthChain.SetBranchAddress("Negative_Px", &m_truth_px_neg, &b_truth_px_neg);
100 m_truthChain.SetBranchAddress("Negative_Py", &m_truth_py_neg, &b_truth_py_neg);
101 m_truthChain.SetBranchAddress("Negative_Pz", &m_truth_pz_neg, &b_truth_pz_neg);
102 }
103
104 return;
105}
106
107
109// initializes histograms
111void ZmumuValidationExample::bookHistograms()
112{
114 // initialize histograms here
116
117 const int fnEtaBins = m_EtaBins;
118 const int fnPhiBins = m_PhiBins;
119 const int fnQoverPtBins = 100;
120
121 const double eta_bound = 2.5;
122 const double phi_bound = M_PI;
123 const double p_bound = 0.5;
124 const double pt_bound = 15;
125
126 const double z_bound = 0.4; //mm
127 const double d_bound = 0.15; //mm
128
129 h_DELTA = new TH1F("h_DELTA",";#delta [GeV]",200,-50,50);
130
131 h_pt = new TH1F("Pt","p_{T} of both #mu; p_{T} [GeV]", 100, 0., 100.);
132 h_pt_pos = new TH1F("Pt_Pos","p_{T} of #mu^{+}; p_{T} [GeV]", 100, 0., 100.);
133 h_pt_neg = new TH1F("Pt_Neg","p_{T} of #mu^{-}; p_{T} [GeV]", 100, 0., 100.);
134
135 h_z0 = new TH1F("h_z0","z_{0}: longitudinal impact param.; z_{0} [mm]", 100, -150., 150.);
136 h_z0_pos = new TH1F("h_z0_Pos","z_{0} of #mu^{+};z_{0} [mm]",100, -150., 150.);
137 h_z0_neg = new TH1F("h_z0_Neg","z_{0} of #mu^{-};z_{0} [mm]",100, -150., 150.);
138 h_d0 = new TH1F("h_d0","d_{0}: transvers. impact param.; d_{0} [mm]", 100, -0.08, 0.08);
139 h_d0_pos = new TH1F("h_d0_Pos","d_{0} of #mu^{+};d_{0} [mm]",100, -0.08, 0.08);
140 h_d0_neg = new TH1F("h_d0_Neg","d_{0} of #mu^{-};d_{0} [mm]",100, -0.08, 0.08);
141
142 h_mass = new TH1F("ZMass",";Mass [GeV]",100,70,110);
143
144// pcorrected_mass = new TH1F("PCorrectedZMass",";Mass [GeV]",100,70,110);
145// lambdacorrected_mass = new TH1F("LambdaCorrectedZMass",";Mass [GeV]",100,70,110);
146
147// pcorrected_mass_vs_etaphi = new TH3F("ZMassVsEtaPhi_P",";#eta;#phi;Mass [GeV]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound, 100,70,110);
148// lambdacorrected_mass_vs_etaphi = new TH3F("ZMassVsEtaPhi_Lambda",";#eta;#phi;Mass [GeV]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound, 100,70,110);
149
150// etaphi_pos = new TH2F("MuonEtaPhi_Positive",";#eta;#phi", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
151// etaphi_neg = new TH2F("MuonEtaPhi_Negative",";#eta;#phi", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
152
153// prof_pt_vs_etaphi = new TProfile2D("AvgPtVsEtaPhi",";#eta;#phi;<p_{T}>", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
154
155 //z0 histograms ----------------------------
156
157 entries = new TH2F("entries","Entries per #eta-#phi bin;#eta;#phi;entries", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
158
159 z0delta_vs_etaphi = new TH3F("Delta_z0_VsEtaPhi",";#eta;#phi;#delta_{z_{0}} [mm]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound, 50,-z_bound,z_bound);
160
161 z0deltacorrections_vs_etaphi = new TH2D("z0CorrectionVsEtaPhi",";#eta;#phi;#delta_{z_{0}} [mm]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
162
163 z0deltacorrections_vs_etaphi_err = new TH2D("z0CorrectionVsEtaPhi_Err",";#eta;#phi;#delta_{z_{0}}", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
164
165
166 z0delta = new TH1D("Delta_z0", ";#delta_{z_{0}} [mm]", 100,-z_bound,z_bound);
167
168 z0delta_pos = new TH1D("Delta_z0_Pos",";#delta_{z_{0}} [mm]", 100,-z_bound,z_bound);
169
170 z0delta_neg = new TH1D("Delta_z0_Neg",";#delta_{z_{0}} [mm]", 100,-z_bound,z_bound);
171
172 z0delta_etaphi = new TH1D("Delta_z0_etaphi",";#delta_{z_{0}} [mm]", 100,-z_bound,z_bound);
173
174 z0delta_etaphi_pos = new TH1D("Delta_z0_etaphi_Pos",";#delta_{z_{0}} [mm]", 100,-z_bound,z_bound);
175
176 z0delta_etaphi_neg = new TH1D("Delta_z0_etaphi_Neg",";#delta_{z_{0}} [mm]", 100,-z_bound,z_bound);
177
178 //d0 Histograms------------------------------------------------------------------
179
180 d0delta_vs_etaphi = new TH3F("Delta_d0_VsEtaPhi",";#eta;#phi;#delta_{d_{0}} [mm]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound, 50,-d_bound, d_bound);
181
182 d0deltacorrections_vs_etaphi = new TH2D("d0CorrectionVsEtaPhi",";#eta;#phi;#delta_{d_{0}} [mm]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
183
184 d0deltacorrections_vs_etaphi_err = new TH2D("d0CorrectionVsEtaPhi_Err",";#eta;#phi;#delta_{d_{0}}", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
185
186 d0delta = new TH1D("Delta_d0", ";#delta_{d_{0}} [mm]", 100,-d_bound,d_bound);
187
188 d0delta_pos = new TH1D("Delta_d0_Pos",";#delta_{d_{0}} [mm]", 100,-d_bound,d_bound);
189
190 d0delta_neg = new TH1D("Delta_d0_Neg",";#delta_{d_{0}} [mm]", 100,-d_bound,d_bound);
191
192 d0delta_etaphi = new TH1D("Delta_d0_etaphi",";#delta_{d_{0}} [mm]", 100,-d_bound,d_bound);
193
194 d0delta_etaphi_pos = new TH1D("Delta_d0_etaphi_Pos",";#delta_{d_{0}} [mm]", 100,-d_bound,d_bound);
195
196 d0delta_etaphi_neg = new TH1D("Delta_d0_etaphi_Neg",";#delta_{d_{0}} [mm]", 100,-d_bound,d_bound);
197
198 //delta - charge symmetric -----------------
199 delta_vs_etaphi = new TH3F("DeltaPVsEtaPhi",";#eta;#phi;#delta_{r}", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound, 100,-p_bound,p_bound);
200
201 deltacorrections_vs_etaphi = new TH2D("PCorrectionVsEtaPhi",";#eta;#phi;#delta_{r}", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
202
203 deltacorrections_vs_etaphi_err = new TH2D("PCorrectionVsEtaPhi_Err",";#eta;#phi;#delta_{r}", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
204
205 //lambda * pT - charge anti-symmetric -----------------
206 lambda_vs_etaphi = new TH3F("LambdaVsEtaPhi",";#eta;#phi;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound, 100,-pt_bound,pt_bound);
207
208 lambdacorrections_vs_etaphi = new TH2D("LambdaCorrectionVsEtaPhi",";#eta;#phi;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
209 lambdacorrections_vs_etaphi_err = new TH2D("LambdaCorrectionVsEtaPhi_Err",";#eta;#phi;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
210 lambdacorrections_vs_etaphi_RMS = new TH2D("LambdaCorrectionVsEtaPhi_RMS",";#eta;#phi;RMS #delta_{sagitta} [GeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
211
212 //lambda * pT - charge anti-symmetric -----------------
213 lambda_vs_eta = new TH2F("LambdaVsEta",";#eta;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, 100,-pt_bound,pt_bound);
214 lambda_vs_eta_pos = new TH2F("LambdaVsEta_Pos", ";#eta;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, 100,-pt_bound,pt_bound);
215 lambda_vs_eta_neg = new TH2F("LambdaVsEta_Neg", ";#eta;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, 100,-pt_bound,pt_bound);
216
217 lambdacorrections_vs_eta = new TH1D("LambdaCorrectionVsEta",";#eta;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound);
218
219 lambda = new TH1D("Lambda", ";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
220 lambda_pos = new TH1D("Lambda_Pos",";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
221 lambda_neg = new TH1D("Lambda_Neg",";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
222
223 lambda_eta = new TH1D("Lambda_eta", ";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
224 lambda_eta_pos = new TH1D("Lambda_eta_Pos",";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
225 lambda_eta_neg = new TH1D("Lambda_eta_Neg",";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
226
227 lambda_etaphi = new TH1D("Lambda_etaphi", ";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
228 lambda_etaphi_pos = new TH1D("Lambda_etaphi_Pos",";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
229 lambda_etaphi_neg = new TH1D("Lambda_etaphi_Neg",";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
230
231 if (m_isMC) {
232
233 h_pt_truth = new TH1D("Pt_Truth", ";p_{T} [GeV]",500,0,500);
234 h_pt_pos_truth = new TH1D("Pt_Pos_Truth",";p_{T} [GeV]",500,0,500);
235 h_pt_neg_truth = new TH1D("Pt_Neg_Truth",";p_{T} [GeV]",500,0,500);
236
237 delta_vs_etaphi_truth = new TH3F("DeltaPVsEtaPhi_Truth",";#eta;#phi;#delta_{r}", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound, 100,-p_bound,p_bound);
238 deltacorrections_vs_etaphi_truth = new TH2D("PCorrectionVsEtaPhi_Truth",";#eta;#phi;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
239 deltacorrections_vs_etaphi_truth_err = new TH2D("PCorrectionVsEtaPhi_Truth_Err",";#eta;#phi;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
240
241 lambda_vs_etaphi_truth = new TH3F("LambdaVsEtaPhi_Truth",";#eta;#phi;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound, 100,-pt_bound,pt_bound);
242 lambdacorrections_vs_etaphi_truth = new TH2D("LambdaCorrectionVsEtaPhi_Truth",";#eta;#phi;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
243 lambdacorrections_vs_etaphi_truth_err = new TH2D("LambdaCorrectionVsEtaPhi_Truth_Err",";#eta;#phi;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, fnPhiBins,-phi_bound,phi_bound);
244
245 lambda_vs_eta_truth = new TH2F("LambdaVsEta_Truth",";#eta;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, 100,-pt_bound,pt_bound);
246 lambdacorrections_vs_eta_truth = new TH1D("LambdaCorrectionVsEta_Truth",";#eta;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound);
247
248 truth_mom_bias_vs_eta = new TH2F("TruthMomentumBiasVsEta",";#eta;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, 100,-pt_bound,pt_bound);
249 truth_mom_biascorrections_vs_eta = new TH1D("TruthMomentumBiasCorrectionVsEta",";#eta;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound);
250
251 truth_mass_bias_vs_eta = new TH2F("TruthMassBiasVsEta",";#eta;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound, 100,-pt_bound,pt_bound);
252 truth_mass_biascorrections_vs_eta = new TH1D("TruthMassBiasCorrectionVsEta",";#eta;#delta_{sagitta} [TeV^{-1}]", fnEtaBins,-eta_bound,eta_bound);
253
254 lambda_truth = new TH1D("Lambda_Truth", ";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
255 lambda_truth_pos = new TH1D("Lambda_Truth_Pos",";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
256 lambda_truth_neg = new TH1D("Lambda_Truth_Neg",";#delta_{sagitta} [TeV^{-1}]", 100,-pt_bound,pt_bound);
257
258 delta_phi_truth = new TH1D("DeltaPhi_Truth",";#phi_{Truth} - #phi_{Rec}", 100,-0.01,0.01);
259
260 delta_eta_truth = new TH1D("DeltaEta_Truth",";#eta_{Truth} - #eta_{Rec}", 100,-0.01,0.01);
261
262 delta_M2_vs_zpt_truth = new TH2D("delta_M2_vs_zpt_truth",";Z p_{T} [GeV]; #Delta (M^2)", 100,0,100, 100, -1.0, 1.0 );
263 delta_M2_vs_zpt = new TH2D("delta_M2_vs_zpt",";Z p_{T} [GeV]; #Delta (M^2)", 100,0,100, 100, -1.0, 1.0 );
264
265 delta_M2_vs_etaphi_pos = new TProfile2D("delta_M2_vs_etaphi_pos",";#eta;#phi; #Delta (M^2)", 20,-2.5,2.5, 20, -3.14, 3.14 );
266 delta_M2_vs_etaphi_neg = new TProfile2D("delta_M2_vs_etaphi_neg",";#eta;#phi; #Delta (M^2)", 20,-2.5,2.5, 20, -3.14, 3.14 );
267
268 }
269
270 // Histograms of q/pt
271 h_QoverPt = new TH2F("h_QoverPt","q/p_{T} whole ID", 2*fnPhiBins, -phi_bound, phi_bound, fnQoverPtBins, -1/pt_bound, 1/pt_bound);
272 h_QoverPt->GetXaxis()->SetTitle("#phi [rad]");
273 h_QoverPt->GetYaxis()->SetTitle("q/p_{T} [GeV^{-1}]");
274
275 h_QoverPt3D = new TH3F("h_QoverPt3D", "q/p_{T} whole ID",
276 fnEtaBins, -eta_bound, eta_bound,
277 2*fnPhiBins, -phi_bound, phi_bound,
278 fnQoverPtBins, -1/pt_bound, 1/pt_bound);
279 h_QoverPt3D->GetXaxis()->SetTitle("#eta");
280 h_QoverPt3D->GetYaxis()->SetTitle("#phi [rad]");
281 h_QoverPt3D->GetZaxis()->SetTitle("q/p_{T} [GeV^{-1}]");
282}
283
284
286// called by user to loop over TChains
288void ZmumuValidationExample::loop( unsigned int maxEvents )
289{
290 std::cout << " TrkValidation::loop()" << std::endl;
291 n_iteration++;
292
293 //clear histograms to be fit
294 this->ResetHistograms();
295
296 unsigned int const nEvents = m_eventChain.GetEntries();
297 unsigned int const maxItr = ( maxEvents > 0 && maxEvents < nEvents ) ? maxEvents : nEvents;
298
299 //loop through events and fill histograms
300 if (m_PrintLevel >= 1) cout << " ** ZmumuValidationExample::loop ** start loop on events ... " << endl;
301 loopThroughEvents( maxItr );
302
303 if (m_PrintLevel >= 1) cout<< " ** ZmumuValidationExample::loop ** Fitting histograms" << endl;
304
305 //Fit corrections and add to previous corrections
306 profileZwithIterativeGaussFit(lambda_vs_etaphi,lambdacorrections_vs_etaphi, lambdacorrections_vs_etaphi_RMS, 1, lambdacorrections_vs_etaphi_err);
307 profileZwithIterativeGaussFit(delta_vs_etaphi,deltacorrections_vs_etaphi, 0, 1,deltacorrections_vs_etaphi_err);
308 profileYwithIterativeGaussFit(lambda_vs_eta,lambdacorrections_vs_eta);
309
310 profileZwithIterativeGaussFit(z0delta_vs_etaphi,z0deltacorrections_vs_etaphi, 0, 1, z0deltacorrections_vs_etaphi_err);
311 profileZwithIterativeGaussFit(d0delta_vs_etaphi,d0deltacorrections_vs_etaphi, 0, 1, d0deltacorrections_vs_etaphi_err);
312
313 if (m_isMC && n_iteration == 1) {
314 profileZwithIterativeGaussFit(lambda_vs_etaphi_truth, lambdacorrections_vs_etaphi_truth, 0, 1, lambdacorrections_vs_etaphi_truth_err);
315 profileZwithIterativeGaussFit(delta_vs_etaphi_truth, deltacorrections_vs_etaphi_truth, 0, 1, deltacorrections_vs_etaphi_truth_err);
316
317 profileYwithIterativeGaussFit(lambda_vs_eta_truth, lambdacorrections_vs_eta_truth);
318 profileYwithIterativeGaussFit(truth_mom_bias_vs_eta, truth_mom_biascorrections_vs_eta);
319 profileYwithIterativeGaussFit(truth_mass_bias_vs_eta, truth_mass_biascorrections_vs_eta);
320 }
321
322 writeToFile(n_iteration);
323
324 return;
325}
326
327
329// loops over each event and calls fillHistograms() if there are tracks
331void ZmumuValidationExample::loopThroughEvents( unsigned int maxItr )
332{
333 std::cout << " Looping over " << maxItr << " events"<< std::endl;
334
335 for( unsigned int eventItr = 0
336 ; eventItr != maxItr
337 ; ++eventItr
338 )
339 {
340 if( eventItr % 10000 == 0 )
341 {
342 std::cout << " Processing event " << eventItr << " of " << maxItr << std::endl;
343 }
344
345 if( !m_eventChain.LoadTree(eventItr) || (m_isMC && !m_truthChain.LoadTree(eventItr)) )
346 continue;
347
348 if( !m_eventChain.GetEntry(eventItr) || (m_isMC && !m_truthChain.GetEntry(eventItr)) )
349 continue;
350
351 fillHistograms();
352 }
353
354 return;
355}
356
358void ZmumuValidationExample::fillHistograms()
359{
360
361 const double muon_mass = ParticleConstants::muonMassInMeV; //MeV
362
363 //create lorentz vectors for both muons
364 TLorentzVector* vec_pos = new TLorentzVector();
365 TLorentzVector* vec_neg = new TLorentzVector();
366
367 //corrected z0 for both muons
368 double corrected_z0_pos;
369 double corrected_z0_neg;
370 double corrected_d0_pos;
371 double corrected_d0_neg;
372
373 vec_pos->SetXYZM(m_px_pos, m_py_pos, m_pz_pos, muon_mass);
374 vec_neg->SetXYZM(m_px_neg, m_py_neg, m_pz_neg, muon_mass);
375
376 entries->Fill(vec_pos->Eta(),vec_pos->Phi(),1);
377 entries->Fill(vec_neg->Eta(),vec_neg->Phi(),1);
378
379 fillHistogram( lambda, vec_pos, vec_neg, 1);
380 fillHistogram( lambda_pos, vec_pos, vec_neg, 1, +1);
381 fillHistogram( lambda_neg, vec_pos, vec_neg, 1, -1);
382
383 //z0 etaphi
384 if ( m_d0_neg != 0 && m_d0_pos != 0){
385 // setting original z0values
386 corrected_z0_pos = m_z0_pos;
387 corrected_z0_neg = m_z0_neg;
388 // filling test histograms before correction
389 fillZd0Histogram( z0delta, corrected_z0_pos, corrected_z0_neg, 0);
390 fillZd0Histogram( z0delta_pos, corrected_z0_pos, corrected_z0_neg, 1);
391 fillZd0Histogram( z0delta_neg, corrected_z0_pos, corrected_z0_neg, -1);
392 // adding correction to z0_pos/neg or d0_pos/neg
393 correctZd0(z0deltacorrections_vs_etaphi, vec_pos, vec_neg, corrected_z0_pos, corrected_z0_neg);
394 // filling the main histogram to fit in the next step
395 fillZd0EtaPhiHistogram(z0delta_vs_etaphi, vec_pos, vec_neg, corrected_z0_pos, corrected_z0_neg);
396 // filling test histograms after correction
397 fillZd0Histogram( z0delta_etaphi, corrected_z0_pos, corrected_z0_neg, 0);
398 fillZd0Histogram( z0delta_etaphi_pos, corrected_z0_pos, corrected_z0_neg, 1);
399 fillZd0Histogram( z0delta_etaphi_neg, corrected_z0_pos, corrected_z0_neg, -1);
400 }
401
402 //d0 etaphi
403 if ( m_d0_neg != 0 && m_d0_pos != 0){
404 // setting original d0values
405 corrected_d0_pos = m_d0_pos;
406 corrected_d0_neg = m_d0_neg;
407 // filling test histograms before correction
408 fillZd0Histogram( d0delta, corrected_d0_pos, corrected_d0_neg, 0);
409 fillZd0Histogram( d0delta_pos, corrected_d0_pos, corrected_d0_neg, 1);
410 fillZd0Histogram( d0delta_neg, corrected_d0_pos, corrected_d0_neg, -1);
411 // adding correction to z0_pos/neg or d0_pos/neg
412 correctZd0(d0deltacorrections_vs_etaphi, vec_pos, vec_neg, corrected_d0_pos, corrected_d0_neg);
413 // filling the main histogram to fit in the next step
414 fillZd0EtaPhiHistogram(d0delta_vs_etaphi, vec_pos, vec_neg, corrected_d0_pos, corrected_d0_neg);
415 // filling test histograms after correction
416 fillZd0Histogram( d0delta_etaphi, corrected_d0_pos, corrected_d0_neg, 0);
417 fillZd0Histogram( d0delta_etaphi_pos, corrected_d0_pos, corrected_d0_neg, 1);
418 fillZd0Histogram( d0delta_etaphi_neg, corrected_d0_pos, corrected_d0_neg, -1);
419 }
420
421 //fill truth histograms on first iteration without correcting momentum
422 if (m_isMC && n_iteration == 1) {
423
424 TLorentzVector* vec_truth_pos = new TLorentzVector();
425 TLorentzVector* vec_truth_neg = new TLorentzVector();
426
427 vec_truth_pos->SetXYZM(m_truth_px_pos, m_truth_py_pos, m_truth_pz_pos, muon_mass);
428 vec_truth_neg->SetXYZM(m_truth_px_neg, m_truth_py_neg, m_truth_pz_neg, muon_mass);
429
430 double z_mass = 91187.6; //MeV
431 double mass = ((*vec_pos) + (*vec_neg)).M();
432 double mass_truth = ((*vec_truth_pos) + (*vec_truth_neg)).M();
433
434 double delta_M2 = (mass*mass - z_mass*z_mass)/(z_mass*z_mass);
435 double delta_M2_truth = (mass*mass - mass_truth*mass_truth)/(mass_truth*mass_truth);
436 double delta_M2_truth_bias = (mass_truth*mass_truth - z_mass*z_mass)/(z_mass*z_mass);
437
438 h_pt_truth->Fill(vec_truth_pos->Pt()/1000.0);
439 h_pt_truth->Fill(vec_truth_neg->Pt()/1000.0);
440 h_pt_pos_truth->Fill(vec_truth_pos->Pt()/1000.0);
441 h_pt_neg_truth->Fill(vec_truth_neg->Pt()/1000.0);
442
443 //truth momentum bias
444 double d_bias_truth_pos = +1/vec_truth_pos->Pt()*(1 - vec_truth_pos->P()/vec_pos->P())*1000000.0;
445 double d_bias_truth_neg = -1/vec_truth_neg->Pt()*(1 - vec_truth_neg->P()/vec_neg->P())*1000000.0;
446
447 truth_mom_bias_vs_eta->Fill(vec_truth_pos->Eta(), d_bias_truth_pos);
448 truth_mom_bias_vs_eta->Fill(vec_truth_neg->Eta(), d_bias_truth_neg);
449
450 //Z mass/method bias
451 double d_mass_bias_truth_pos = +1*delta_M2_truth_bias/vec_truth_pos->Pt()*1000000.0;
452 double d_mass_bias_truth_neg = -1*delta_M2_truth_bias/vec_truth_neg->Pt()*1000000.0;
453
454 truth_mass_bias_vs_eta->Fill(vec_truth_pos->Eta(), d_mass_bias_truth_pos);
455 truth_mass_bias_vs_eta->Fill(vec_truth_neg->Eta(), d_mass_bias_truth_neg);
456
457 //lambda corrections
458 double d_lambda_truth_pos = +1*delta_M2_truth/vec_truth_pos->Pt()*1000000.0;
459 double d_lambda_truth_neg = -1*delta_M2_truth/vec_truth_neg->Pt()*1000000.0;
460
461 lambda_vs_etaphi_truth->Fill(vec_truth_pos->Eta(), vec_truth_pos->Phi(), d_lambda_truth_pos);
462 lambda_vs_etaphi_truth->Fill(vec_truth_neg->Eta(), vec_truth_neg->Phi(), d_lambda_truth_neg);
463
464 lambda_vs_eta_truth->Fill(vec_truth_pos->Eta(), d_lambda_truth_pos);
465 lambda_vs_eta_truth->Fill(vec_truth_neg->Eta(), d_lambda_truth_neg);
466
467 lambda_truth->Fill( d_lambda_truth_pos );
468 lambda_truth->Fill( d_lambda_truth_neg );
469
470 lambda_truth_pos->Fill( d_lambda_truth_pos );
471 lambda_truth_neg->Fill( d_lambda_truth_neg );
472
473 //delta corrections
474 double delta_truth_pos = (vec_truth_pos->P() - vec_pos->P())/vec_truth_pos->P();
475 double delta_truth_neg = (vec_truth_neg->P() - vec_neg->P())/vec_truth_neg->P();
476
477 delta_vs_etaphi_truth->Fill(vec_truth_pos->Eta(), vec_truth_pos->Phi(), delta_truth_pos);
478 delta_vs_etaphi_truth->Fill(vec_truth_neg->Eta(), vec_truth_neg->Phi(), delta_truth_neg);
479
480 //check eta/phi biases
481 delta_phi_truth->Fill(vec_truth_pos->Phi() - vec_pos->Phi());
482 delta_phi_truth->Fill(vec_truth_neg->Phi() - vec_neg->Phi());
483
484 delta_eta_truth->Fill(vec_truth_pos->Eta() - vec_pos->Eta());
485 delta_eta_truth->Fill(vec_truth_neg->Eta() - vec_neg->Eta());
486
487 //other checks
488 double zpt_truth = ((*vec_truth_pos) + (*vec_truth_neg)).Pt();
489
490 h_DELTA->Fill(delta_M2/(vec_pos->Pt() - vec_neg->Pt())*1000000);
491 //h_DELTA->Fill(delta_M2/(vec_truth_pos->Pt() - vec_truth_neg->Pt())*1000000);
492
493 delta_M2_vs_zpt_truth->Fill(zpt_truth/1000.0,delta_M2_truth);
494 delta_M2_vs_zpt->Fill(zpt_truth/1000.0,delta_M2);
495
496 delta_M2_vs_etaphi_pos->Fill(vec_pos->Eta(), vec_pos->Phi(), delta_M2);
497 delta_M2_vs_etaphi_neg->Fill(vec_neg->Eta(), vec_neg->Phi(), delta_M2);
498
499 delete vec_truth_pos;
500 delete vec_truth_neg;
501
502 } // end if of MC and truth info
503
504 //-----------------------------------------------------------------
505 //apply lambda momentum corrections in eta/phi
506 correctMomentum( lambdacorrections_vs_etaphi, vec_pos, vec_neg, 1);
507
508 fillHistogram( lambda_etaphi, vec_pos, vec_neg, 1);
509 fillHistogram( lambda_etaphi_pos, vec_pos, vec_neg, 1, +1);
510 fillHistogram( lambda_etaphi_neg, vec_pos, vec_neg, 1, -1);
511
512 //fill histograms
513 fillEtaPhiHistogram( lambda_vs_etaphi, vec_pos, vec_neg, 1);
514
515 // fill q/pt histograms
516 fillQoverPtHistograms (vec_pos, vec_neg);
517
518 //-----------------------------------------------------------------
519
520 //reset muon vectors and use original values
521 vec_pos->SetXYZM(m_px_pos, m_py_pos, m_pz_pos, muon_mass);
522 vec_neg->SetXYZM(m_px_neg, m_py_neg, m_pz_neg, muon_mass);
523
524 //apply delta momentum corrections in eta/phi
525 correctMomentum( deltacorrections_vs_etaphi, vec_pos, vec_neg, 0);
526 //fill histograms
527 fillEtaPhiHistogram( delta_vs_etaphi, vec_pos, vec_neg, 0);
528
529 //-----------------------------------------------------------------
530
531 //reset muon vectors and use original values
532 vec_pos->SetXYZM(m_px_pos, m_py_pos, m_pz_pos, muon_mass);
533 vec_neg->SetXYZM(m_px_neg, m_py_neg, m_pz_neg, muon_mass);
534
535 //apply lambda momentum corrections in eta
536 correctMomentum( lambdacorrections_vs_eta, vec_pos, vec_neg, 1);
537
538 fillHistogram( lambda_eta, vec_pos, vec_neg, 1);
539 fillHistogram( lambda_eta_pos, vec_pos, vec_neg, 1, +1);
540 fillHistogram( lambda_eta_neg, vec_pos, vec_neg, 1, -1);
541
542 //fill histograms
543 fillEtaHistogram( lambda_vs_eta, vec_pos, vec_neg, 1);
544 fillEtaHistogram( lambda_vs_eta_pos, vec_pos, vec_neg, 1, +1);
545 fillEtaHistogram( lambda_vs_eta_neg, vec_pos, vec_neg, 1, -1);
546
547 //check corrected pt/mass distributions at each iteration
548 h_pt->Fill(vec_pos->Pt()/1000.0);
549 h_pt->Fill(vec_neg->Pt()/1000.0);
550 h_pt_pos->Fill(vec_pos->Pt()/1000.0);
551 h_pt_neg->Fill(vec_neg->Pt()/1000.0);
552
553 h_mass->Fill(((*vec_pos)+(*vec_neg)).M()/1000.0);
554
555 //check z0 distributions
556 h_z0->Fill(m_z0_pos);
557 h_z0->Fill(m_z0_neg);
558 h_z0_pos->Fill(m_z0_pos);
559 h_z0_neg->Fill(m_z0_neg);
560 //check d0 distributions
561 h_d0->Fill(m_d0_pos);
562 h_d0->Fill(m_d0_neg);
563 h_d0_pos->Fill(m_d0_pos);
564 h_d0_neg->Fill(m_d0_neg);
565 //-----------------------------------------------------------------
566
567 delete vec_pos;
568 delete vec_neg;
569
570 return;
571}
572
573//--------------------------------------------------------------------------------------------------
574void ZmumuValidationExample::correctMomentum( TH1* h_corrections, TLorentzVector* muon_pos, TLorentzVector* muon_neg, int use_lambda)
575{
576 double temp_pt_pos = muon_pos->Pt();
577 double temp_pt_neg = muon_neg->Pt();
578
579 // positive muon
580 double correction = h_corrections->GetBinContent(h_corrections->FindBin(muon_pos->Eta(), muon_pos->Phi()));
581 if (use_lambda) {
582
583 double pt_true = muon_pos->Pt();
584 //double pt_true = muon_pos->Pt()/(1+(+1)*correction*muon_pos->Pt()/1000000.0);
585
586 muon_pos->SetXYZM(muon_pos->Px()/(1+correction*pt_true/1000000.0),
587 muon_pos->Py()/(1+correction*pt_true/1000000.0),
588 muon_pos->Pz()/(1+correction*pt_true/1000000.0),
589 muon_pos->M());
590 if (m_PrintLevel >= 3) {
591 cout << " ** correctMomentum ** mu+ (eta,phi) (" << muon_pos->Eta() << ", " << muon_pos->Phi() << ") "
592 << " bin = " << h_corrections->FindBin(muon_pos->Eta(), muon_pos->Phi())
593 << " correction = " << correction
594 << endl;
595 }
596 }
597 else {
598 // mainly for delta corrections
599 muon_pos->SetXYZM(muon_pos->Px()*(1-correction),
600 muon_pos->Py()*(1-correction),
601 muon_pos->Pz()*(1-correction),
602 muon_pos->M());
603 }
604
605 // negative muon
606 correction = h_corrections->GetBinContent(h_corrections->FindBin(muon_neg->Eta(), muon_neg->Phi()));
607 if (use_lambda) {
608
609 double pt_true = muon_neg->Pt();
610 //double pt_true = muon_neg->Pt()/(1+(+1)*correction*muon_neg->Pt()/1000000.0);
611
612 muon_neg->SetXYZM(muon_neg->Px()/(1-correction*pt_true/1000000.0),
613 muon_neg->Py()/(1-correction*pt_true/1000000.0),
614 muon_neg->Pz()/(1-correction*pt_true/1000000.0),
615 muon_neg->M());
616 if (m_PrintLevel >= 3) {
617 cout << " ** correctMomentum ** mu- (eta,phi) (" << muon_neg->Eta() << ", " << muon_neg->Phi() << ") "
618 << " bin = " << h_corrections->FindBin(muon_neg->Eta(), muon_neg->Phi())
619 << " correction = " << correction
620 << endl;
621 }
622 }
623 else {
624 // mainly for delta corrections
625 muon_neg->SetXYZM(muon_neg->Px()*(1-correction),
626 muon_neg->Py()*(1-correction),
627 muon_neg->Pz()*(1-correction),
628 muon_neg->M());
629 }
630
631 if (m_PrintLevel >= 3) {
632 cout << " ** correctMomentum ** mu+ Pt: " << temp_pt_pos << " --> " << muon_pos->Pt() << " Delta = " << 100*(muon_pos->Pt()-temp_pt_pos)/temp_pt_pos << " %"<< endl
633 << " mu- Pt: " << temp_pt_neg << " --> " << muon_neg->Pt() << " Delta = " << 100*(muon_neg->Pt()-temp_pt_neg)/temp_pt_neg << " %"<< endl;
634 }
635
636
637 return;
638}
639
640
641//-------------------------------------------------------------------------------
642void ZmumuValidationExample::fillEtaPhiHistogram(TH3* hist, TLorentzVector* v_pos, TLorentzVector* v_neg, int use_lambda)
643{
644 double z_mass = 91187.6 + m_shift; //MeV
645 double mass = ((*v_pos) + (*v_neg)).M();
646 double delta_M2 = (mass*mass - z_mass*z_mass)/(z_mass*z_mass);
647
648 if (use_lambda) {
649 // to estimate the correction factor from: pt -> pt /(1+q pt delta_sagitta): q DeltaM2/pt
650 hist->Fill(v_pos->Eta(), v_pos->Phi(), +1*m_factor*delta_M2/v_pos->Pt()*1000000.0);
651 hist->Fill(v_neg->Eta(), v_neg->Phi(), -1*m_factor*delta_M2/v_neg->Pt()*1000000.0);
652
653 } else {
654 hist->Fill(v_pos->Eta(), v_pos->Phi(), +1*delta_M2/2);
655 hist->Fill(v_neg->Eta(), v_neg->Phi(), +1*delta_M2/2);
656 }
657}
658
659
660//-------------------------------------------------------------------------------
661void ZmumuValidationExample::fillEtaHistogram(TH2* hist, TLorentzVector* v_pos, TLorentzVector* v_neg, int use_lambda, int charge)
662{
663 double z_mass = 91187.6 + m_shift; //MeV
664 double mass = ((*v_pos) + (*v_neg)).M();
665 double delta_M2 = (mass*mass - z_mass*z_mass)/(z_mass*z_mass);
666
667 if (use_lambda) {
668 if (charge == 0 || charge == +1) hist->Fill(v_pos->Eta(), +1*m_factor*delta_M2/v_pos->Pt()*1000000.0 );
669 if (charge == 0 || charge == -1) hist->Fill(v_neg->Eta(), -1*m_factor*delta_M2/v_neg->Pt()*1000000.0 );
670
671 } else {
672 hist->Fill(v_pos->Eta(), +1*delta_M2/2);
673 hist->Fill(v_neg->Eta(), +1*delta_M2/2);
674 }
675}
676
677//-----------------------------------------------------------------------------------------------
678void ZmumuValidationExample::fillHistogram(TH1* hist, TLorentzVector* v_pos, TLorentzVector* v_neg, int fill_lambda, int charge)
679{
680 double z_mass = 91187.6 + m_shift; //MeV
681 double mass = ((*v_pos) + (*v_neg)).M();
682 double delta_M2 = (mass*mass - z_mass*z_mass)/(z_mass*z_mass);
683
684 if (fill_lambda) {
685 if (charge == 0 || charge == +1) hist->Fill( +1*m_factor*delta_M2/v_pos->Pt()*1000000.0 );
686 if (charge == 0 || charge == -1) hist->Fill( -1*m_factor*delta_M2/v_neg->Pt()*1000000.0 );
687
688 } else {
689 hist->Fill(+1*delta_M2/2);
690 hist->Fill(+1*delta_M2/2);
691 }
692}
693
694//-----------------------------------------------------------------------------------------------
695void ZmumuValidationExample::fillQoverPtHistograms(TLorentzVector* v_pos, TLorentzVector* v_neg)
696{
697 h_QoverPt->Fill((*v_pos).Phi(), 1000/(*v_pos).Pt()); // 1000 is for MeV to GeV
698 h_QoverPt->Fill((*v_neg).Phi(), -1000/(*v_neg).Pt());
699
700 h_QoverPt3D->Fill((*v_pos).Eta(), (*v_pos).Phi(), 1000/(*v_pos).Pt()); // 1000 is for MeV to GeV
701 h_QoverPt3D->Fill((*v_neg).Eta(), (*v_neg).Phi(), -1000/(*v_neg).Pt());
702 return;
703}
704// impact parameter biases methods------------------------------------------------------------------
705void ZmumuValidationExample::correctZd0( TH1* h_corrections, TLorentzVector* muon_pos, TLorentzVector* muon_neg, double& zd0_muon_p, double& zd0_muon_n)
706{
707 double correction = h_corrections->GetBinContent(h_corrections->FindBin(muon_pos->Eta(), muon_pos->Phi()));
708// std::cout << "Pos Cor" << correction << std::endl;
709
710 zd0_muon_p = zd0_muon_p + correction ;
711
712 correction = h_corrections->GetBinContent(h_corrections->FindBin(muon_neg->Eta(), muon_neg->Phi()));
713 zd0_muon_n = zd0_muon_n + correction ;
714// std::cout << "Neg Cor" << correction << std::endl;
715
716}
717
718void ZmumuValidationExample::fillZd0EtaPhiHistogram(TH3* hist, TLorentzVector* v_pos, TLorentzVector* v_neg, double zd0_muon_p, double zd0_muon_n)
719{
720 double deltazd0 = -(zd0_muon_p - zd0_muon_n);
721 hist->Fill(v_pos->Eta(),v_pos->Phi(), deltazd0);
722 hist->Fill(v_neg->Eta(),v_neg->Phi(), -1.*deltazd0);
723}
724
725void ZmumuValidationExample::fillZd0Histogram(TH1* hist, double zd0_muon_p, double zd0_muon_n, int pn)
726{
727 double deltazd0 = (zd0_muon_p - zd0_muon_n);
728 if (pn == +1) hist->Fill(deltazd0);
729 if (pn == -1) hist->Fill(-1.*deltazd0);
730 if (pn == 0) {
731 hist->Fill(-1.*deltazd0);
732 hist->Fill(deltazd0);
733 }
734}
735// impact parameter biases methods------------------------------------------------------------------
736
738// saves root objects to a TFile
740void ZmumuValidationExample::writeToFile(int iteration)
741{
742
743 m_file.cd();
744
745 if (m_isMC && iteration == 1) {
746
747 h_pt_truth->Write();
748 h_pt_pos_truth->Write();
749 h_pt_neg_truth->Write();
750
751 deltacorrections_vs_etaphi_truth->Write();
752 deltacorrections_vs_etaphi_truth_err->Write();
753
754 lambdacorrections_vs_etaphi_truth->Write();
755 lambdacorrections_vs_etaphi_truth_err->Write();
756
757 lambdacorrections_vs_eta_truth->Write();
758 lambda_vs_eta_truth->Write();
759
760 truth_mom_biascorrections_vs_eta->Write();
761 truth_mom_bias_vs_eta->Write();
762
763 truth_mass_biascorrections_vs_eta->Write();
764 truth_mass_bias_vs_eta->Write();
765
766 lambda_truth->Write();
767 lambda_truth_pos->Write();
768 lambda_truth_neg->Write();
769
770 delta_phi_truth->Write();
771 delta_eta_truth->Write();
772
773 delta_M2_vs_zpt_truth->Write();
774 delta_M2_vs_zpt->Write();
775
776 delta_M2_vs_etaphi_pos->Write();
777 delta_M2_vs_etaphi_neg->Write();
778
779 }
780
781 if (iteration > 0) {
782 m_file.mkdir(Form("Iteration%i",iteration));
783 m_file.cd(Form("Iteration%i",iteration));
784 }
785
786 h_pt->Write();
787 h_pt_pos->Write();
788 h_pt_neg->Write();
789 h_mass->Write();
790
791 h_z0->Write();
792 h_z0_pos->Write();
793 h_z0_neg->Write();
794 h_d0->Write();
795 h_d0_pos->Write();
796 h_d0_neg->Write();
797
798 deltacorrections_vs_etaphi->Write();
799 deltacorrections_vs_etaphi_err->Write();
800
801 entries->Write();
802
803 lambdacorrections_vs_etaphi->Write();
804 lambdacorrections_vs_etaphi_err->Write();
805 lambdacorrections_vs_etaphi_RMS->Write();
806 lambdacorrections_vs_eta->Write();
807
808 lambda_vs_eta->Write();
809 lambda_vs_etaphi->Write();
810
811 lambda_vs_eta_pos->Write();
812 lambda_vs_eta_neg->Write();
813
814 lambda->Write();
815 lambda_pos->Write();
816 lambda_neg->Write();
817
818 lambda_eta->Write();
819 lambda_eta_pos->Write();
820 lambda_eta_neg->Write();
821
822 lambda_etaphi->Write();
823 lambda_etaphi_pos->Write();
824 lambda_etaphi_neg->Write();
825
826 z0delta_vs_etaphi->Write();
827 z0deltacorrections_vs_etaphi->Write();
828 z0deltacorrections_vs_etaphi_err->Write();
829
830 z0delta->Write();
831 z0delta_pos->Write();
832 z0delta_neg->Write();
833
834
835 z0delta_etaphi->Write();
836 z0delta_etaphi_pos->Write();
837 z0delta_etaphi_neg->Write();
838
839 d0delta_vs_etaphi->Write();
840 d0deltacorrections_vs_etaphi->Write();
841 d0deltacorrections_vs_etaphi_err->Write();
842
843 d0delta->Write();
844 d0delta_pos->Write();
845 d0delta_neg->Write();
846
847 d0delta_etaphi->Write();
848 d0delta_etaphi_pos->Write();
849 d0delta_etaphi_neg->Write();
850
851 //cout<< "MEAN: " <<h_DELTA->GetMean()<<endl;
852 h_DELTA->Write();
853
854 //
855 h_QoverPt->Write();
856 h_QoverPt3D->Write();
857 // delta_vs_etaphi->Write();
858}
859
860
861//-------------------------------------------------------------
862void ZmumuValidationExample::profileZwithIterativeGaussFit(TH3* hist, TH2* mu_graph, TH2* sigma_graph, int num_bins, TH2* mu_err_graph, TH2* sigma_err_graph)
863{
864 if (!hist) {
865 cout<< "ProfileZwithIterativeGaussFit(): No histogram supplied!"<<endl;
866 return;
867 }
868
869 int minEntries = 50;
870 int fDebug = 0;
871
872 int num_bins_x = hist->GetXaxis()->GetNbins();
873 int num_bins_y = hist->GetYaxis()->GetNbins();
874
875 double num_not_converged = 0;
876 int num_skipped = 0;
877
878 double max_sigma = 0;
879 double min_sigma = 0;
880
881 double max_mu = 0;
882 double min_mu = 0;
883
884 TH1D* current_proj;
885
886 for (int i = 1; i < num_bins_x+(num_bins==1); i+=num_bins) {
887
888 for (int j = 1; j < num_bins_y+(num_bins==1); j+=num_bins) {
889
890 int index = i/num_bins;
891 int index_y = j/num_bins;
892
893 current_proj = hist->ProjectionZ(Form("%s_GaussProjection_%i_%i",hist->GetName(),index, index_y),i,i+num_bins-1,j,j+num_bins-1);
894 current_proj->SetTitle(Form("%s - Bin %i x %i",hist->GetName(), index,index_y));
895
896 double current_mu,current_err_mu, current_sigma, current_err_sigma;
897
898 if(current_proj->GetEntries() < minEntries) {
899 if (m_PrintLevel >= 1) cout << " ** profileZwithIterativeGaussFit ** fitting " << hist->GetName() << " bin (" << index << ", " << index_y << ") "
900 << " Not enough entries " << current_proj->GetEntries() << " < " << minEntries << endl;
901 //current_mu = -999;
902 current_mu = 0;
903 current_sigma = 0;
904 current_err_mu = 1;
905 current_err_sigma = 1;
906
907 if (fDebug) std::cout<<"WARNING: Only "<<current_proj->GetEntries()<<" entries in bin "<<index<<","<<index_y<< " in histogram " <<hist->GetName()<< std::endl;
908 num_skipped++;
909
910 } else {
911 if (m_PrintLevel >= 1) cout << " ** profileZwithIterativeGaussFit ** fitting " << hist->GetName() << " bin (" << index << ", " << index_y << ") " << endl;
912
913 IterativeGaussFit(current_proj, current_mu, current_err_mu, current_sigma, current_err_sigma);
914
915 if (current_sigma > max_sigma || max_sigma == 0) max_sigma = current_sigma;
916 if (current_sigma < min_sigma || min_sigma == 0) min_sigma = current_sigma;
917 if (current_mu > max_mu || max_mu == 0) max_mu = current_mu;
918 if (current_mu < min_mu || min_mu == 0) min_mu = current_mu;
919
920 }//end if entries < minEntries
921
922 float x_coord = (hist->GetXaxis()->GetBinLowEdge(i) + hist->GetXaxis()->GetBinUpEdge(i+num_bins-1))/2;
923 float y_coord = (hist->GetYaxis()->GetBinLowEdge(j) + hist->GetYaxis()->GetBinUpEdge(j+num_bins-1))/2;
924
925 if (sigma_graph) sigma_graph->Fill(x_coord,y_coord,current_sigma);
926 if (mu_graph) mu_graph->Fill(x_coord,y_coord,current_mu);
927
928 //should probably be replace bin content, not fill?
929 if (sigma_err_graph) sigma_err_graph->Fill(x_coord,y_coord,current_err_sigma);
930 if (mu_err_graph) mu_err_graph->Fill(x_coord,y_coord,current_err_mu);
931
932 delete current_proj;
933
934 } //end loop on j (y)
935 } //end loop on i (x)
936
937 if (mu_graph) {
938 mu_graph->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
939 mu_graph->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
940 mu_graph->GetYaxis()->SetTitleOffset(1);
941 mu_graph->GetZaxis()->SetTitle(hist->GetZaxis()->GetTitle());
942 mu_graph->GetZaxis()->SetTitleOffset(1.2);
943 mu_graph->SetTitle( "" );
944 //mu_graph->SetMaximum(max_mu + 0.1* (max_mu-min_mu));
945 //mu_graph->SetMinimum(min_mu - 0.1* (max_mu-min_mu));
946 }
947
948 if (sigma_graph) {
949 sigma_graph->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
950 sigma_graph->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
951 sigma_graph->GetYaxis()->SetTitleOffset(1);
952 sigma_graph->GetZaxis()->SetTitle(hist->GetZaxis()->GetTitle());
953 sigma_graph->GetZaxis()->SetTitleOffset(1.2);
954 sigma_graph->SetTitle( "" );
955 //sigma_graph->SetMaximum(max_sigma + 0.1* (max_sigma-min_sigma));
956 //sigma_graph->SetMinimum(min_sigma - 0.1* (max_sigma-min_sigma));
957 }
958
959 if (mu_err_graph) {
960 mu_err_graph->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
961 mu_err_graph->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
962 mu_err_graph->GetYaxis()->SetTitleOffset(1);
963 mu_err_graph->GetZaxis()->SetTitle(Form("Error of fit #mu: %s",hist->GetZaxis()->GetTitle()));
964 mu_err_graph->GetZaxis()->SetTitleOffset(1.2);
965 mu_err_graph->SetTitle(hist->GetTitle());
966 //mu_err_graph->SetMaximum(max_mu + 0.1* (max_mu-min_mu));
967 //mu_err_graph->SetMinimum(min_mu - 0.1* (max_mu-min_mu));
968 }
969
970 if (sigma_err_graph) {
971 sigma_err_graph->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
972 sigma_err_graph->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
973 sigma_err_graph->GetYaxis()->SetTitleOffset(1);
974 sigma_err_graph->GetZaxis()->SetTitle(Form("Error of fit #sigma: %s",hist->GetZaxis()->GetTitle()));
975 sigma_err_graph->GetZaxis()->SetTitleOffset(1.2);
976 sigma_err_graph->SetTitle(hist->GetTitle());
977 //sigma_err_graph->SetMaximum(max_mu + 0.1* (max_mu-min_mu));
978 //sigma_err_graph->SetMinimum(min_mu - 0.1* (max_mu-min_mu));
979 }
980
981
982 if (num_not_converged || num_skipped) std::cout<<"Fit Results for histogram: "<< hist->GetName()<<std::endl;
983 if (num_not_converged) std::cout<<"Non Convergent Bin Fraction: "<<num_not_converged<< " / " <<num_bins_x*num_bins_y - num_skipped<<std::endl;
984 if (num_skipped) std::cout<<"Number skipped bins: "<<num_skipped<< " / " <<num_bins_x*num_bins_y<<std::endl;
985
986 return;
987}
988
989//-----------------------------------------------------------------------------
990void ZmumuValidationExample::profileYwithIterativeGaussFit(TH2* hist, TH1* mu_graph, TH1* sigma_graph, int num_bins)
991{
992
993 if (!hist) {
994 std::cout << "Error in ProfileYwithIterativeGaussFit(): Histogram not found" <<std::endl;
995 return;
996 }
997
998 if (num_bins < 1 ) {
999 std::cout << "Error in ProfileYwithIterativeGaussFit(): Invalid number of bins to integrate over." <<std::endl;
1000 return;
1001 }
1002
1003 const int minEntries = 50;
1004 const int fDebug = 0;
1005
1006 int num_bins_x = hist->GetXaxis()->GetNbins();
1007
1008 if (mu_graph) mu_graph->Rebin(num_bins);
1009 if (sigma_graph) sigma_graph->Rebin(num_bins);
1010
1011 double* errs_mu = new double[num_bins_x/num_bins + 2]; // +2 for overflow!!
1012 double* errs_sigma = new double[num_bins_x/num_bins + 2];
1013
1014 errs_mu[0] = 0;
1015 errs_mu[num_bins_x/num_bins + 1] = 0;
1016
1017 errs_sigma[0] = 0;
1018 errs_sigma[num_bins_x/num_bins + 1] = 0;
1019
1020 double min_sigma = 0;
1021 double max_sigma = 0;
1022 double min_mu = 0;
1023 double max_mu = 0;
1024
1025 int num_skipped = 0;
1026
1027 TH1D* current_proj;
1028
1029 for (int i = 1; i < (num_bins_x + (num_bins == 1)); i+=num_bins) {
1030
1031 int index = i/num_bins;
1032 if (num_bins == 1) index--;
1033
1034 current_proj = hist->ProjectionY(Form("%s_projection_%i",hist->GetName(),index),i,i+num_bins-1);
1035
1036 double mu, mu_err, sigma, sigma_err;
1037
1038 if(current_proj->GetEntries() < minEntries) {
1039 mu = 0;
1040 mu_err = 0;
1041 sigma = 0;
1042 sigma_err = 0;
1043 num_skipped++;
1044 if ( fDebug ) std::cout<<"WARNING: Not enough entries in bin "<<index<<std::endl;
1045 } else {
1046
1047 IterativeGaussFit(current_proj, mu, mu_err, sigma, sigma_err);
1048
1049 if (sigma > max_sigma || max_sigma == 0) max_sigma = sigma;
1050 if (sigma < min_sigma || min_sigma == 0) min_sigma = sigma;
1051 if (mu > max_mu || max_mu == 0) max_mu = mu;
1052 if (mu < min_mu || min_mu == 0) min_mu = mu;
1053
1054 }
1055
1056 double value_x = (hist->GetXaxis()->GetBinLowEdge(i) + hist->GetXaxis()->GetBinUpEdge(i+num_bins-1))/2;
1057
1058 //Important!! Use Fill to increment the graph with each iteration, or SetBinContent to replace contents...
1059
1060 //if (sigma_graph) sigma_graph->SetBinContent(i, sigma);
1061 //if (mu_graph) mu_graph->SetBinContent(i, mu);
1062
1063 if (sigma_graph) sigma_graph->Fill(value_x, sigma);
1064 if (mu_graph) mu_graph->Fill(value_x, mu);
1065
1066 errs_mu[index + 1] = mu_err;
1067 errs_sigma[index + 1] = sigma_err;
1068
1069 delete current_proj;
1070 }
1071
1072 if (sigma_graph) {
1073 sigma_graph->SetError(errs_sigma);
1074 //sigma_graph->SetMaximum(max_sigma+0.15*(max_sigma - min_sigma));
1075 //sigma_graph->SetMinimum(min_sigma-0.15*(max_sigma - min_sigma));
1076 sigma_graph->GetYaxis()->SetTitleOffset(1.5);
1077 sigma_graph->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
1078 sigma_graph->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
1079 sigma_graph->SetTitle("");
1080 }
1081
1082 if (mu_graph) {
1083 mu_graph->SetError(errs_mu);
1084 //mu_graph->SetMaximum(max_mu+0.15*(max_mu - min_mu));
1085 //mu_graph->SetMinimum(min_mu-0.15*(max_mu - min_mu));
1086 mu_graph->GetYaxis()->SetTitleOffset(1.5);
1087 mu_graph->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
1088 mu_graph->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
1089 mu_graph->SetTitle("");
1090 }
1091
1092 if (fDebug && num_skipped) std::cout<<" Number of skipped bins: "<<num_skipped<<std::endl;
1093
1094 return;
1095
1096}
1097
1098//-----------------------------------------------------------------------------
1099int ZmumuValidationExample::IterativeGaussFit(TH1* hist, double &mu, double &mu_err, double &sigma, double &sigma_err)
1100{
1101
1102 //constants for fitting algorithm
1103 const int iteration_limit = 20;
1104 const float percent_limit = 0.01;
1105 const float fit_range_multiplier = 1.5;
1106 const int fDebug = 0;
1107
1108 double last_mu;
1109 double last_sigma;
1110 double current_mu;
1111 double current_sigma;
1112 double mu_percent_diff;
1113 double sigma_percent_diff;
1114
1115 if (!hist) {
1116 if (fDebug) std::cout<< "Error in Anp::IterativeGaussFit(): Histogram to be fit is missing" <<std::endl;
1117 return 1;
1118 }
1119
1120 this->HistogramConditioning(hist);
1121
1122 TF1* fit_func = new TF1("fit_func","gaus");
1123
1124 int bad_fit = hist->Fit(fit_func,"QN");
1125
1126 if (fDebug && bad_fit) std::cout <<"BAD INITIAL FIT: "<< hist->GetTitle() << std::endl;
1127
1128 last_mu = fit_func->GetParameter(1);
1129 last_sigma = fit_func->GetParameter(2);
1130
1131 if (bad_fit) last_mu = hist->GetMean();
1132
1133 // check as well that the last_mu is reasonable
1134 if (fabs(last_mu - hist->GetMean()) > 5*hist->GetBinWidth(1)) last_mu = hist->GetMean();
1135
1136 fit_func->SetRange(last_mu-fit_range_multiplier*last_sigma,last_mu+fit_range_multiplier*last_sigma);
1137
1138 int iteration = 0;
1139
1140 while ( iteration < iteration_limit ) {
1141
1142 iteration++;
1143
1144 double FitRangeLower = last_mu-fit_range_multiplier*last_sigma;
1145 double FitRangeUpper = last_mu+fit_range_multiplier*last_sigma;
1146
1147 // if range is to narrow --> broaden it
1148 if ((FitRangeUpper-FitRangeLower)/hist->GetBinWidth(1) < 4) {
1149 FitRangeLower -= hist->GetBinWidth(1);
1150 FitRangeUpper += hist->GetBinWidth(1);
1151 }
1152
1153 fit_func->SetRange(FitRangeLower, FitRangeUpper);
1154 if (m_PrintLevel >= 3) cout << " ** IterativeGaussFit ** fit iter # " << iteration
1155 << " new fit range: " << FitRangeLower << " --> " << FitRangeUpper << endl;
1156
1157
1158
1159 bad_fit = hist->Fit(fit_func, "RQN");
1160
1161 if (fDebug && bad_fit) std::cout<<" ** BAD FIT ** : bin "<< hist->GetTitle() <<" iteration "<<iteration<<std::endl;
1162
1163 current_mu = fit_func->GetParameter(1);
1164 current_sigma = fit_func->GetParameter(2);
1165
1166 //std::cout<<"Iteration: "<<iteration<<" Current: "<<current_mu<<" "<<current_sigma<<" Last: "<<last_mu<<" "<<last_sigma<<std::endl;
1167
1168 float average_mu = (last_mu+current_mu)/2;
1169 float average_sigma = (last_sigma+current_sigma)/2;
1170
1171 if (average_mu == 0) {
1172 if ( fDebug ) std::cout<<" Average mu = 0 in bin "<< hist->GetTitle() <<std::endl;
1173 average_mu = current_mu;
1174 }
1175
1176 if (average_sigma == 0) {
1177 if ( fDebug ) std::cout<<"Average sigma = 0; Fit Problem in "<< hist->GetTitle() <<". "<<last_sigma<<" "<<current_sigma<<std::endl;
1178 average_sigma = current_sigma;
1179 }
1180
1181 mu_percent_diff = fabs((last_mu-current_mu)/average_mu);
1182 sigma_percent_diff = fabs((last_sigma-current_sigma)/average_sigma);
1183
1184 if ( mu_percent_diff < percent_limit && sigma_percent_diff < percent_limit ) break;
1185
1186 if (iteration != iteration_limit) { //necessary?
1187 last_mu = current_mu;
1188 last_sigma = current_sigma;
1189 }
1190 // check as well that the last_mu is reasonable
1191 if (fabs(last_mu - hist->GetMean()) > 5*hist->GetBinWidth(1)) {
1192 if (m_PrintLevel >= 3) cout << " ** IterativeGaussFit ** fit iter # " << iteration
1193 << " ** WARNING ** last_mu looks bad: " << last_mu
1194 << " this iter mu: " << fit_func->GetParameter(1)
1195 << " proposed mu: " << hist->GetMean()
1196 << endl;
1197 last_mu = hist->GetMean();
1198 }
1199 }
1200
1201 if (iteration == iteration_limit) {
1202 if (fDebug ) std::cout<<"WARNING: Fit did not converge to < "<<percent_limit*100<<"% in "<< hist->GetTitle() <<" in "<<iteration_limit<<" iterations. Percent Diffs: "<<mu_percent_diff*100<<"% (Mu),"<<" "<<sigma_percent_diff*100<<"% (Sigma)"<<std::endl;
1203 //possibly return a value other than 0 to indicate not converged?
1204 }
1205
1206 mu = current_mu;
1207 mu_err = fit_func->GetParError(1);
1208 sigma = current_sigma;
1209 sigma_err = fit_func->GetParError(2);
1210
1211 hist->GetListOfFunctions()->Add(fit_func);
1212
1213 if (m_PrintLevel >= 1 ) {
1214 cout << " ** IterativeGaussFit ** fit result: histo name " << hist->GetName() << " title: " << hist->GetTitle() << endl
1215 << " mu = " << mu << " +- " << mu_err << endl
1216 << " sigma = " << sigma << " +- " << sigma_err
1217 << endl;
1218 if (TempCanvasIterGaussFit == NULL) {
1219 TempCanvasIterGaussFit = new TCanvas ("TempCanvasIterGaussFit","Iter Gauss fit", 400, 400);
1220 }
1221 hist->DrawCopy();
1222 TempCanvasIterGaussFit->Update();
1223 hist->Print();
1224 string input = "";
1225 cout << " ** IterativeGaussFit ** Please type RETURN to continue:\n>";
1226 getline(cin, input);
1227 }
1228
1229 return 0;
1230}
1231
1232//-----------------------------------------------------------------------------
1233void ZmumuValidationExample::SetPrintLevel (int newprintlevel)
1234{
1235 m_PrintLevel = newprintlevel;
1236 if (m_PrintLevel < 0) m_PrintLevel = 0;
1237 return;
1238}
1239
1240//-----------------------------------------------------------------------------
1241void ZmumuValidationExample::SetEtaBins (int newbins)
1242{
1243 m_EtaBins = newbins;
1244 if (m_EtaBins < 1) m_EtaBins = 1;
1245 if (m_EtaBins > 100) m_EtaBins = 100;
1246 return;
1247}
1248
1249//-----------------------------------------------------------------------------
1250void ZmumuValidationExample::SetPhiBins (int newbins)
1251{
1252 m_PhiBins = newbins;
1253 if (m_PhiBins < 1) m_PhiBins = 1;
1254 if (m_PhiBins > 100) m_PhiBins = 100;
1255 return;
1256}
1257
1258//-----------------------------------------------------------------------------
1259void ZmumuValidationExample::DrawMap ()
1260{
1261 this->SymmetrizeHisto(lambdacorrections_vs_etaphi);
1262
1263 lambdacorrections_vs_etaphi->DrawCopy("colz");
1264
1265 return;
1266}
1267
1268//-----------------------------------------------------------------------------
1269void ZmumuValidationExample::ResetHistograms ()
1270{
1271 delta_vs_etaphi->Reset();
1272
1273 //
1274 entries->Reset();
1275 z0delta_vs_etaphi->Reset();
1276
1277 z0delta->Reset();
1278 z0delta_pos->Reset();
1279 z0delta_neg->Reset();
1280 z0delta_etaphi->Reset();
1281 z0delta_etaphi_pos->Reset();
1282 z0delta_etaphi_neg->Reset();
1283
1284 d0delta_vs_etaphi->Reset();
1285
1286 d0delta->Reset();
1287 d0delta_pos->Reset();
1288 d0delta_neg->Reset();
1289 d0delta_etaphi->Reset();
1290 d0delta_etaphi_pos->Reset();
1291 d0delta_etaphi_neg->Reset();
1292 //
1293
1294 h_DELTA->Reset();
1295 h_pt->Reset();
1296 h_pt_pos->Reset();
1297 h_pt_neg->Reset();
1298 h_mass->Reset();
1299
1300 delta_vs_etaphi->Reset();
1301
1302 deltacorrections_vs_etaphi->Reset();
1303 deltacorrections_vs_etaphi_err->Reset();
1304
1305 lambda_vs_etaphi->Reset();
1306 lambdacorrections_vs_etaphi->Reset();
1307 lambdacorrections_vs_etaphi_err->Reset();
1308
1309 lambda_vs_eta->Reset();
1310 lambda_vs_eta_pos->Reset();
1311 lambda_vs_eta_neg->Reset();
1312
1313 lambdacorrections_vs_eta->Reset();
1314
1315 lambda->Reset();
1316 lambda_pos->Reset();
1317 lambda_neg->Reset();
1318
1319 lambda_eta->Reset();
1320 lambda_eta_pos->Reset();
1321 lambda_eta_neg->Reset();
1322
1323 lambda_etaphi->Reset();
1324 lambda_etaphi_pos->Reset();
1325 lambda_etaphi_neg->Reset();
1326
1327 // impact param histos
1328 h_z0->Reset();
1329 h_z0_pos->Reset();
1330 h_z0_neg->Reset();
1331 h_d0->Reset();
1332 h_d0_pos->Reset();
1333 h_d0_neg->Reset();
1334
1335 return;
1336}
1337
1338//-----------------------------------------------------------------------------
1339void ZmumuValidationExample::SymmetrizeHisto (TH2* hist)
1340{
1341 double RangeUpper = hist->GetBinContent(hist->GetMaximumBin());
1342 double RangeLower = hist->GetBinContent(hist->GetMinimumBin());
1343
1344 double NewRangeUpper = RangeUpper;
1345 double NewRangeLower = -RangeUpper;
1346
1347 if (RangeUpper < - RangeLower) {
1348 NewRangeUpper = -RangeLower;
1349 NewRangeLower = RangeLower;
1350 }
1351
1352 NewRangeUpper *= 1.01; // increase a bit the scale just to make sure everything fits in
1353 NewRangeLower *= 1.01;
1354
1355 if (m_PrintLevel >= 3) {
1356 cout << " ** SymmetrizeHisto ** old range: " << RangeLower << " --> " << RangeUpper << endl;
1357 cout << " new range: " << NewRangeLower << " --> " << NewRangeUpper << endl;
1358 }
1359
1360 hist->SetMaximum(NewRangeUpper);
1361 hist->SetMinimum(NewRangeLower);
1362
1363 return;
1364}
1365
1366//-----------------------------------------------------------------------------
1367void ZmumuValidationExample::HistogramConditioning (TH1* hist)
1368{
1369 if (m_PrintLevel>=3) cout << " ** HistogramConditioning ** START ** hist = " << hist->GetName() << endl;
1370
1371 double MinEntriesMPB = 15;
1372 Int_t NGroupBins = 2;
1373
1374 // goal:
1375 // make sure that the most populated bin has a minimum number of entries
1376 Int_t MostPopulatedBin = (hist->GetMaximumBin());
1377 double EntriesMPB = hist->GetBinContent(MostPopulatedBin);
1378 if (EntriesMPB < MinEntriesMPB) {
1379 // check the entries of the neighbour channels
1380 if ((EntriesMPB + hist->GetBinContent(MostPopulatedBin+1) + hist->GetBinContent(MostPopulatedBin-1)) > MinEntriesMPB) {
1381 NGroupBins = 2;
1382 }
1383 else {
1384 NGroupBins = 3;
1385 }
1386
1387 // now find the first divisor (factor of ) the original number of bins
1388 Bool_t DivisorFound = false;
1389 while (!DivisorFound) {
1390 if ( hist->GetNbinsX() % NGroupBins == 0) {
1391 DivisorFound = true;
1392 }
1393 else {
1394 DivisorFound = false;
1395 NGroupBins++;
1396 }
1397 }
1398 Int_t NBinsWas = hist->GetNbinsX();
1399 hist = hist->Rebin(NGroupBins);
1400 if (m_PrintLevel>=1) cout << " ** HistogramConditioning ** histogram had to be rebinned by: " << NGroupBins
1401 << " NBins was: " << NBinsWas << " and now is: " << hist->GetNbinsX() << endl;
1402
1403 }
1404
1405
1406 return;
1407}
M_PI
#define M_PI
Definition ActiveFraction.h:11
charge
double charge(const T &p)
Definition AtlasPID.h:997
ParticleConstants.h
A number of constexpr particle constants to avoid hardcoding them directly in various places.
ZmumuValidationExample::loop
void loop(unsigned maxEvents=0)
Definition ZmumuValidationExample.cxx:288
ZmumuValidationExample::correctZd0
void correctZd0(TH1 *h_corrections, TLorentzVector *muon_pos, TLorentzVector *muon_neg, double &zd0_muon_p, double &zd0_muon_n)
Definition ZmumuValidationExample.cxx:705
ZmumuValidationExample::profileZwithIterativeGaussFit
void profileZwithIterativeGaussFit(TH3 *hist, TH2 *mu_graph, TH2 *sigma_graph, int num_bins, TH2 *mu_err_graph=0, TH2 *sigma_err_graph=0)
Definition ZmumuValidationExample.cxx:862
ZmumuValidationExample::fillQoverPtHistograms
void fillQoverPtHistograms(TLorentzVector *v_pos, TLorentzVector *v_neg)
Definition ZmumuValidationExample.cxx:695
ZmumuValidationExample::correctMomentum
void correctMomentum(TH1 *h_corrections, TLorentzVector *muon_pos, TLorentzVector *muon_neg, int use_lambda)
Definition ZmumuValidationExample.cxx:574
ZmumuValidationExample::IterativeGaussFit
int IterativeGaussFit(TH1 *hist, double &mu, double &mu_err, double &sigma, double &sigma_err)
Definition ZmumuValidationExample.cxx:1099
ZmumuValidationExample::fillHistogram
void fillHistogram(TH1 *hist, TLorentzVector *v_pos, TLorentzVector *v_neg, int fill_lambda, int charge=0)
Definition ZmumuValidationExample.cxx:678
ZmumuValidationExample::fillEtaHistogram
void fillEtaHistogram(TH2 *hist, TLorentzVector *v_pos, TLorentzVector *v_neg, int use_lambda, int charge=0)
Definition ZmumuValidationExample.cxx:661
ZmumuValidationExample::m_fileNames
const std::list< std::string > m_fileNames
Definition ZmumuValidationExample.h:26
ZmumuValidationExample::ZmumuValidationExample
ZmumuValidationExample(std::list< std::string > const &s_fileNames, string s_treeName="DefaultParams", std::string const &s_outFileName="ZmumuValidationExampleOutput.root", bool isMC=false)
Definition ZmumuValidationExample.cxx:16
ZmumuValidationExample::loopThroughEvents
void loopThroughEvents(unsigned int maxItr)
Definition ZmumuValidationExample.cxx:331
ZmumuValidationExample::profileYwithIterativeGaussFit
void profileYwithIterativeGaussFit(TH2 *hist, TH1 *mu_graph=0, TH1 *sigma_graph=0, int num_bins=1)
Definition ZmumuValidationExample.cxx:990
ZmumuValidationExample::SetPrintLevel
void SetPrintLevel(int newprintlevel=0)
Definition ZmumuValidationExample.cxx:1233
ZmumuValidationExample::fillZd0EtaPhiHistogram
void fillZd0EtaPhiHistogram(TH3 *hist, TLorentzVector *v_pos, TLorentzVector *v_neg, double z0_muon_p, double z0_muon_n)
Definition ZmumuValidationExample.cxx:718
ZmumuValidationExample::fillEtaPhiHistogram
void fillEtaPhiHistogram(TH3 *hist, TLorentzVector *v_pos, TLorentzVector *v_neg, int use_lambda)
Definition ZmumuValidationExample.cxx:642
ZmumuValidationExample::fillZd0Histogram
void fillZd0Histogram(TH1 *hist, double z0_muon_p, double z0_muon_n, int pn)
Definition ZmumuValidationExample.cxx:725
nEvents
const int nEvents
Definition fbtTestBasics.cxx:78
ParticleConstants::PDG2011::muonMassInMeV
constexpr double muonMassInMeV
the mass of the muon (in MeV)
Definition ParticleConstants.h:29
index
Definition index.py:1
Generated on for ATLAS Offline Software by doxygen 1.14.0