]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG/FLOW/Tasks/AliJetFlowTools.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
bugfix in function which selects unfolding routine
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliJetFlowTools.cxx
... / ...
CommitLineData
1/*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16// Author: Redmer Alexander Bertens, Utrecht University, Utrecht, Netherlands
17// (rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl)
18//
19// Tools class for Jet Flow Analysis, replaces 'extractJetFlow.C' macro
20//
21// The task uses input from two analysis tasks:
22// - $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetV2.cxx
23// used to retrieve jet spectra and delta pt distributions
24// - $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx
25// used to construct the detector response function
26// and unfolds jet spectra with respect to the event plane. The user can choose
27// different alrogithms for unfolding which are available in (ali)root. RooUnfold
28// libraries must be present on the system
29// ( see http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html ).
30//
31// The weak spot of this class is the function PrepareForUnfolding, which will read
32// output from two output files and expects histograms with certain names and binning.
33// Unfolding methods itself are general and should be able to handle any input, therefore one
34// can forgo the PrepareForUnfolding method, and supply necessary input information via the
35// SetRawInput() method
36//
37// to see an example of how to use this class, see $ALICE_ROOT/PWGCF/FLOW/macros/jetFlowTools.C
38
39// root includes
40#include "TH1.h"
41#include "TF2.h"
42#include "TH2D.h"
43#include "TGraph.h"
44#include "TGraphErrors.h"
45#include "TGraphAsymmErrors.h"
46#include "TLine.h"
47#include "TCanvas.h"
48#include "TLegend.h"
49#include "TArrayD.h"
50#include "TList.h"
51#include "TMinuit.h"
52#include "TVirtualFitter.h"
53#include "TLegend.h"
54#include "TCanvas.h"
55#include "TStyle.h"
56#include "TLine.h"
57#include "TMath.h"
58#include "TVirtualFitter.h"
59#include "TFitResultPtr.h"
60#include "Minuit2/Minuit2Minimizer.h"
61#include "Math/Functor.h"
62// aliroot includes
63#include "AliUnfolding.h"
64#include "AliAnaChargedJetResponseMaker.h"
65// class includes
66#include "AliJetFlowTools.h"
67// roo unfold includes (make sure you have these available on your system)
68#include "RooUnfold.h"
69#include "RooUnfoldResponse.h"
70#include "RooUnfoldSvd.h"
71#include "RooUnfoldBayes.h"
72#include "TSVDUnfold.h"
73
74using namespace std;
75//_____________________________________________________________________________
76AliJetFlowTools::AliJetFlowTools() :
77 fResponseMaker (new AliAnaChargedJetResponseMaker()),
78 fRMS (kTRUE),
79 fSymmRMS (kTRUE),
80 fRho0 (kFALSE),
81 fBootstrap (kFALSE),
82 fPower (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
83 fSaveFull (kTRUE),
84 fActiveString (""),
85 fActiveDir (0x0),
86 fInputList (0x0),
87 fRefreshInput (kTRUE),
88 fOutputFileName ("UnfoldedSpectra.root"),
89 fOutputFile (0x0),
90 fCentralityArray (0x0),
91 fMergeBinsArray (0x0),
92 fCentralityWeights (0x0),
93 fDetectorResponse (0x0),
94 fJetFindingEff (0x0),
95 fBetaIn (.1),
96 fBetaOut (.1),
97 fBayesianIterIn (4),
98 fBayesianIterOut (4),
99 fBayesianSmoothIn (0.),
100 fBayesianSmoothOut (0.),
101 fAvoidRoundingError (kFALSE),
102 fUnfoldingAlgorithm (kChi2),
103 fPrior (kPriorMeasured),
104 fPriorUser (0x0),
105 fBinsTrue (0x0),
106 fBinsRec (0x0),
107 fBinsTruePrior (0x0),
108 fBinsRecPrior (0x0),
109 fSVDRegIn (5),
110 fSVDRegOut (5),
111 fSVDToy (kTRUE),
112 fJetRadius (0.3),
113 fEventCount (-1),
114 fNormalizeSpectra (kFALSE),
115 fSmoothenPrior (kFALSE),
116 fFitMin (60.),
117 fFitMax (300.),
118 fFitStart (75.),
119 fSmoothenCounts (kTRUE),
120 fTestMode (kFALSE),
121 fRawInputProvided (kFALSE),
122 fEventPlaneRes (.63),
123 fUseDetectorResponse(kTRUE),
124 fUseDptResponse (kTRUE),
125 fTrainPower (kTRUE),
126 fDphiUnfolding (kTRUE),
127 fDphiDptUnfolding (kFALSE),
128 fExLJDpt (kTRUE),
129 fTitleFontSize (-999.),
130 fRMSSpectrumIn (0x0),
131 fRMSSpectrumOut (0x0),
132 fRMSRatio (0x0),
133 fRMSV2 (0x0),
134 fDeltaPtDeltaPhi (0x0),
135 fJetPtDeltaPhi (0x0),
136 fSpectrumIn (0x0),
137 fSpectrumOut (0x0),
138 fDptInDist (0x0),
139 fDptOutDist (0x0),
140 fDptIn (0x0),
141 fDptOut (0x0),
142 fFullResponseIn (0x0),
143 fFullResponseOut (0x0) { // class constructor
144 // create response maker weight function (tuned to PYTHIA spectrum)
145 fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0, 200));
146 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
147}
148//_____________________________________________________________________________
149void AliJetFlowTools::Make(TH1* customIn, TH1* customOut) {
150 // core function of the class
151 if(fDphiDptUnfolding) {
152 // to extract the yield as function of Dphi, Dpt - experimental
153 MakeAU();
154 return;
155 }
156 // 1) rebin the raw output of the jet task to the desired binnings
157 // 2) calls the unfolding routine
158 // 3) writes output to file
159 // can be repeated multiple times with different configurations
160
161 // 1) manipulation of input histograms
162 // check if the input variables are present
163 if(fRefreshInput) {
164 if(!PrepareForUnfolding(customIn, customOut)) {
165 printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
166 return;
167 }
168 }
169 // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
170 // parts of the spectrum can end up in over or underflow bins
171
172 // if bootstrap mode is kTRUE, resample the underlying distributions
173 // FIXME think about resampling the rebinned results or raw results, could lead to difference
174 // in smoothness of tail of spectrum (which is probably not used in any case, but still ... )
175/*
176 if(fBootstrap) {
177 // resample but leave original spectra intact for the next unfolding round
178 fSpectrumIn = reinterpret_cast<TH1D*>(Bootstrap(fSpectrumIn, kFALSE));
179 fSpectrumOut = reinterpret_cast<TH1D*>(Bootstrap(fSpectrumOut, kFALSE));
180 }
181*/
182 TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
183 TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec, TString("resized_out_"), kFALSE);
184
185 if(fBootstrap) {
186 measuredJetSpectrumIn = reinterpret_cast<TH1D*>(Bootstrap(measuredJetSpectrumIn, kFALSE));
187 measuredJetSpectrumOut = reinterpret_cast<TH1D*>(Bootstrap(measuredJetSpectrumOut, kFALSE));
188 }
189 // for now do it BEFORE as after gives an issue in Rebin function (counts are wrong)
190
191
192
193 // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
194 // the template will be used as a prior for the chi2 unfolding
195 TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
196 TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
197 // get the full response matrix from the dpt and the detector response
198 fDetectorResponse = NormalizeTH2D(fDetectorResponse);
199 // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
200 // so that unfolding should return the initial spectrum
201 if(!fTestMode) {
202 if(fUseDptResponse && fUseDetectorResponse) {
203 fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
204 fFullResponseOut = MatrixMultiplication(fDptOut, fDetectorResponse);
205 } else if (fUseDptResponse && !fUseDetectorResponse) {
206 fFullResponseIn = fDptIn;
207 fFullResponseOut = fDptOut;
208 } else if (!fUseDptResponse && fUseDetectorResponse) {
209 fFullResponseIn = fDetectorResponse;
210 fFullResponseOut = fDetectorResponse;
211 } else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) {
212 printf(" > No response, exiting ! < \n" );
213 return;
214 }
215 } else {
216 fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
217 fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
218 }
219 // normalize each slide of the response to one
220 NormalizeTH2D(fFullResponseIn);
221 NormalizeTH2D(fFullResponseOut);
222 // resize to desired binning scheme
223 TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
224 TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
225 // get the kinematic efficiency
226 TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
227 kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
228 TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
229 kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
230 // suppress the errors
231 for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
232 kinematicEfficiencyIn->SetBinError(1+i, 0.);
233 kinematicEfficiencyOut->SetBinError(1+i, 0.);
234 }
235 TH1D* jetFindingEfficiency(0x0);
236 if(fJetFindingEff) {
237 jetFindingEfficiency = ProtectHeap(fJetFindingEff);
238 jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
239 jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
240 }
241 // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
242 TH1D* unfoldedJetSpectrumIn(0x0);
243 TH1D* unfoldedJetSpectrumOut(0x0);
244 fActiveDir->cd(); // select active dir
245 TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
246 dirIn->cd(); // select inplane subdir
247 // do the inplane unfolding
248 unfoldedJetSpectrumIn = UnfoldWrapper(
249 measuredJetSpectrumIn,
250 resizedResponseIn,
251 kinematicEfficiencyIn,
252 measuredJetSpectrumTrueBinsIn,
253 TString("in"),
254 jetFindingEfficiency);
255 resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
256 resizedResponseIn->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
257 resizedResponseIn->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
258 resizedResponseIn = ProtectHeap(resizedResponseIn);
259 resizedResponseIn->Write();
260 kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
261 kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
262 kinematicEfficiencyIn->Write();
263 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
264 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
265 fDetectorResponse->Write();
266 // optional histograms
267 if(fSaveFull) {
268 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
269 fSpectrumIn->Write();
270 fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
271 fDptInDist->Write();
272 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
273 fDptIn->Write();
274 fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
275 fFullResponseIn->Write();
276 }
277 fActiveDir->cd();
278 if(fDphiUnfolding) {
279 TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
280 dirOut->cd();
281 // do the out of plane unfolding
282 unfoldedJetSpectrumOut = UnfoldWrapper(
283 measuredJetSpectrumOut,
284 resizedResponseOut,
285 kinematicEfficiencyOut,
286 measuredJetSpectrumTrueBinsOut,
287 TString("out"),
288 jetFindingEfficiency);
289 resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
290 resizedResponseOut->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
291 resizedResponseOut->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
292 resizedResponseOut = ProtectHeap(resizedResponseOut);
293 resizedResponseOut->Write();
294 kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
295 kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
296 kinematicEfficiencyOut->Write();
297 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
298 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
299 fDetectorResponse->Write();
300 if(jetFindingEfficiency) jetFindingEfficiency->Write();
301 // optional histograms
302 if(fSaveFull) {
303 fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
304 fSpectrumOut->Write();
305 fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
306 fDptOutDist->Write();
307 fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
308 fDptOut->Write();
309 fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
310 fFullResponseOut->Write();
311 }
312
313 // write general output histograms to file
314 fActiveDir->cd();
315 if(unfoldedJetSpectrumIn && unfoldedJetSpectrumOut && unfoldedJetSpectrumIn && unfoldedJetSpectrumOut) {
316 TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out")));
317 if(ratio) {
318 ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
319 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
320 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
321 ratio = ProtectHeap(ratio);
322 ratio->Write();
323 // write histo values to RMS files if both routines converged
324 // input values are weighted by their uncertainty
325 for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
326 if(unfoldedJetSpectrumIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumIn->GetBinError(i+1), 2.));
327 if(unfoldedJetSpectrumOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), unfoldedJetSpectrumOut->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumOut->GetBinError(i+1), 2.));
328 if(unfoldedJetSpectrumOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1) / unfoldedJetSpectrumOut->GetBinContent(i+1));
329 }
330 }
331 TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
332 if(v2) {
333 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
334 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
335 v2->GetYaxis()->SetTitle("v_{2}");
336 v2 = ProtectHeap(v2);
337 v2->Write();
338 }
339 } else if (unfoldedJetSpectrumOut && unfoldedJetSpectrumIn) {
340 TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
341 if(ratio) {
342 ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
343 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
344 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
345 ratio = ProtectHeap(ratio);
346 ratio->Write();
347 }
348 TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
349 if(v2) {
350 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
351 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
352 v2->GetYaxis()->SetTitle("v_{2}");
353 v2 = ProtectHeap(v2);
354 v2->Write();
355 }
356 }
357 } // end of if(fDphiUnfolding)
358 fDeltaPtDeltaPhi->Write();
359 unfoldedJetSpectrumIn->Sumw2();
360 ProtectHeap(unfoldedJetSpectrumIn, kFALSE);
361 unfoldedJetSpectrumIn->Write();
362 unfoldedJetSpectrumOut->Sumw2();
363 ProtectHeap(unfoldedJetSpectrumOut, kFALSE);
364 unfoldedJetSpectrumOut->Write();
365 fJetPtDeltaPhi->Write();
366 // save the current state of the unfolding object
367 SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
368 TH1D* unfoldedJetSpectrumInForSub((TH1D*)unfoldedJetSpectrumIn->Clone("forSubIn"));
369 TH1D* unfoldedJetSpectrumOutForSub((TH1D*)unfoldedJetSpectrumOut->Clone("forSubOut"));
370 unfoldedJetSpectrumInForSub->Add(unfoldedJetSpectrumOutForSub, -1.);
371 unfoldedJetSpectrumInForSub= ProtectHeap(unfoldedJetSpectrumInForSub);
372 unfoldedJetSpectrumInForSub->Write();
373
374}
375//_____________________________________________________________________________
376TH1D* AliJetFlowTools::UnfoldWrapper(
377 const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
378 const TH2D* resizedResponse, // response matrix
379 const TH1D* kinematicEfficiency, // kinematic efficiency
380 const TH1D* measuredJetSpectrumTrueBins, // unfolding template: same binning is pt gen of response
381 const TString suffix, // suffix (in or out of plane)
382 const TH1D* jetFindingEfficiency) // jet finding efficiency
383{
384 // wrapper function to call specific unfolding routine
385 TH1D* (AliJetFlowTools::*myFunction)(const TH1D*, const TH2D*, const TH1D*, const TH1D*, const TString, const TH1D*);
386 // initialize functon pointer
387 if(fUnfoldingAlgorithm == kChi2) myFunction = &AliJetFlowTools::UnfoldSpectrumChi2;
388 else if(fUnfoldingAlgorithm == kBayesian) myFunction = &AliJetFlowTools::UnfoldSpectrumBayesian;
389 else if(fUnfoldingAlgorithm == kBayesianAli) myFunction = &AliJetFlowTools::UnfoldSpectrumBayesianAli;
390 else if(fUnfoldingAlgorithm == kSVD) myFunction = &AliJetFlowTools::UnfoldSpectrumSVD;
391 else if(fUnfoldingAlgorithm == kFold) myFunction = &AliJetFlowTools::FoldSpectrum;
392 else if(fUnfoldingAlgorithm == kNone) {
393 TH1D* clone((TH1D*)measuredJetSpectrum->Clone("clone"));
394 clone->SetNameTitle(Form("MeasuredJetSpectrum%s", suffix.Data()), Form("measuredJetSpectrum %s", suffix.Data()));
395 return clone;//RebinTH1D(clone, fBinsTrue, clone->GetName(), kFALSE);
396 }
397 else return 0x0;
398 // do the actual unfolding with the selected function
399 return (this->*myFunction)( measuredJetSpectrum, resizedResponse, kinematicEfficiency, measuredJetSpectrumTrueBins, suffix, jetFindingEfficiency);
400}
401//_____________________________________________________________________________
402TH1D* AliJetFlowTools::UnfoldSpectrumChi2(
403 const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
404 const TH2D* resizedResponse, // response matrix
405 const TH1D* kinematicEfficiency, // kinematic efficiency
406 const TH1D* measuredJetSpectrumTrueBins, // unfolding template: same binning is pt gen of response
407 const TString suffix, // suffix (in or out of plane)
408 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
409{
410 // unfold the spectrum using chi2 minimization
411
412 // step 0) setup the static members of AliUnfolding
413 ResetAliUnfolding(); // reset from previous iteration
414 // also deletes and re-creates the global TVirtualFitter
415 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
416 if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
417 else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
418 if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
419 else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
420 AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
421
422 // step 1) clone all input histograms. the histograms are cloned to make sure that the original histograms
423 // stay intact. a local copy of a histogram (which only exists in the scope of this function) is
424 // denoted by the suffix 'Local'
425
426 // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
427 TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
428 // unfolded local will be filled with the result of the unfolding
429 TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
430
431 // full response matrix and kinematic efficiency
432 TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
433 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
434
435 // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
436 TH1D *priorLocal = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("priorLocal_%s", suffix.Data()));
437 // optionally, the prior can be smoothened by extrapolating the spectrum using a power law fit
438 if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
439
440 // step 2) start the unfolding
441 Int_t status(-1), i(0);
442 while(status < 0 && i < 100) {
443 // i > 0 means that the first iteration didn't converge. in that case, the result of the first
444 // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
445 if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
446 status = AliUnfolding::Unfold(
447 resizedResponseLocal, // response matrix
448 kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
449 measuredJetSpectrumLocal, // measured spectrum
450 priorLocal, // initial conditions (set NULL to use measured spectrum)
451 unfoldedLocal); // results
452 // status holds the minuit fit status (where 0 means convergence)
453 i++;
454 }
455 // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
456 TH2D* hPearson(0x0);
457 if(status == 0 && gMinuit->fISW[1] == 3) {
458 // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
459 TVirtualFitter *fitter(TVirtualFitter::GetFitter());
460 if(gMinuit) gMinuit->Command("SET COV");
461 TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
462 TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
463 pearson->Print();
464 hPearson = new TH2D(*pearson);
465 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
466 hPearson = ProtectHeap(hPearson);
467 hPearson->Write();
468 if(fMergeBinsArray) unfoldedLocal = MergeSpectrumBins(fMergeBinsArray, unfoldedLocal, hPearson);
469 } else status = -1;
470
471 // step 3) refold the unfolded spectrum and save the ratio measured / refolded
472 TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
473 foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
474 unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
475 TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
476 if(ratio) {
477 ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
478 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
479 ratio = ProtectHeap(ratio);
480 ratio->Write();
481 }
482
483 // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
484 // objects are cloned using 'ProtectHeap()'
485 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
486 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
487 measuredJetSpectrumLocal->Write();
488
489 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
490 resizedResponseLocal->Write();
491
492 unfoldedLocal = ProtectHeap(unfoldedLocal);
493 if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
494 unfoldedLocal->Write();
495
496 foldedLocal = ProtectHeap(foldedLocal);
497 foldedLocal->Write();
498
499 priorLocal = ProtectHeap(priorLocal);
500 priorLocal->Write();
501
502 // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
503 TH1F* fitStatus(new TH1F(Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), 4, -0.5, 3.5));
504 fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
505 fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
506 fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
507 fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
508 fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
509 fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
510 fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
511 fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
512 fitStatus->Write();
513
514 return unfoldedLocal;
515}
516//_____________________________________________________________________________
517TH1D* AliJetFlowTools::UnfoldSpectrumSVD(
518 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
519 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
520 const TH1D* kinematicEfficiency, // kinematic efficiency
521 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
522 const TString suffix, // suffix (in, out)
523 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
524{
525
526 TH1D* priorLocal( GetPrior(
527 measuredJetSpectrum, // jet pt in pt rec bins
528 resizedResponse, // full response matrix, normalized in slides of pt true
529 kinematicEfficiency, // kinematic efficiency
530 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
531 suffix, // suffix (in, out)
532 jetFindingEfficiency)); // jet finding efficiency (optional)
533 if(!priorLocal) {
534 printf(" > couldn't find prior ! < \n");
535 return 0x0;
536 } else printf(" 1) retrieved prior \n");
537
538 // go back to the 'root' directory of this instance of the SVD unfolding routine
539 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
540
541 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
542 // measured jets in pt rec binning
543 TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
544 // local copie of the response matrix
545 TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
546 // local copy of response matrix, all true slides normalized to 1
547 // this response matrix will eventually be used in the re-folding routine
548 TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
549 resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
550 // kinematic efficiency
551 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
552 // place holder histos
553 TH1D *unfoldedLocalSVD(0x0);
554 TH1D *foldedLocalSVD(0x0);
555 cout << " 2) setup necessary input " << endl;
556 // 3) configure routine
557 RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
558 cout << " step 3) configured routine " << endl;
559
560 // 4) get transpose matrices
561 // a) get the transpose of the full response matrix
562 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
563 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
564 // normalize it with the prior. this will ensure that high statistics bins will constrain the
565 // end result most strenuously than bins with limited number of counts
566 responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
567 cout << " 4) retrieved first transpose matrix " << endl;
568
569 // 5) get response for SVD unfolding
570 RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
571 cout << " 5) retrieved roo unfold response object " << endl;
572
573 // 6) actualy unfolding loop
574 RooUnfoldSvd unfoldSVD(&responseSVD, measuredJetSpectrumLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
575 unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
576 // correct the spectrum for the kinematic efficiency
577 unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
578
579 // get the pearson coefficients from the covariance matrix
580 TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
581 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
582 TH2D* hPearson(0x0);
583 if(pearson) {
584 hPearson = new TH2D(*pearson);
585 pearson->Print();
586 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
587 hPearson = ProtectHeap(hPearson);
588 hPearson->Write();
589 if(fMergeBinsArray) unfoldedLocalSVD = MergeSpectrumBins(fMergeBinsArray, unfoldedLocalSVD, hPearson);
590 } else return 0x0; // return if unfolding didn't converge
591
592 // plot singular values and d_i vector
593 TSVDUnfold* svdUnfold(unfoldSVD.Impl());
594 TH1* hSVal(svdUnfold->GetSV());
595 TH1D* hdi(svdUnfold->GetD());
596 hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
597 hSVal->SetXTitle("singular values");
598 hSVal->Write();
599 hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
600 hdi->SetXTitle("|d_{i}^{kreg}|");
601 hdi->Write();
602 cout << " plotted singular values and d_i vector " << endl;
603
604 // 7) refold the unfolded spectrum
605 foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, resizedResponseLocalNorm, kinematicEfficiencyLocal);
606 TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalSVD, "ratio measured / re-folded", kTRUE));
607 ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
608 ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
609 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
610 ratio->Write();
611 cout << " 7) refolded the unfolded spectrum " << endl;
612
613 // write the measured, unfolded and re-folded spectra to the output directory
614 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
615 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
616 measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
617 measuredJetSpectrumLocal->Write(); // input spectrum
618 unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
619 unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
620 if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
621 unfoldedLocalSVD->Write(); // unfolded spectrum
622 foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
623 foldedLocalSVD = ProtectHeap(foldedLocalSVD);
624 foldedLocalSVD->Write(); // re-folded spectrum
625
626 // save more general bookkeeeping histograms to the output directory
627 responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
628 responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
629 responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
630 responseMatrixLocalTransposePrior->Write();
631 priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
632 priorLocal->SetXTitle("p_{t} [GeV/c]");
633 priorLocal = ProtectHeap(priorLocal);
634 priorLocal->Write();
635 resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
636 resizedResponseLocalNorm->Write();
637
638 // save some info
639 TH1F* fitStatus(new TH1F(Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), 1, -0.5, 0.5));
640 fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
641 fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
642 fitStatus->Write();
643
644 return unfoldedLocalSVD;
645}
646//_____________________________________________________________________________
647TH1D* AliJetFlowTools::UnfoldSpectrumBayesianAli(
648 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
649 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
650 const TH1D* kinematicEfficiency, // kinematic efficiency
651 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
652 const TString suffix, // suffix (in, out)
653 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
654{
655 // unfold the spectrum using the bayesian unfolding impelmented in AliUnfolding
656 // FIXME careful, not tested yet ! (06122013) FIXME
657
658 // step 0) setup the static members of AliUnfolding
659 ResetAliUnfolding(); // reset from previous iteration
660 // also deletes and re-creates the global TVirtualFitter
661 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
662 if(!strcmp("in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
663 else if(!strcmp("out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
664 else if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
665 else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
666 AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
667
668 // 1) get a prior for unfolding and clone all the input histograms
669 TH1D* priorLocal( GetPrior(
670 measuredJetSpectrum, // jet pt in pt rec bins
671 resizedResponse, // full response matrix, normalized in slides of pt true
672 kinematicEfficiency, // kinematic efficiency
673 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
674 suffix, // suffix (in, out)
675 jetFindingEfficiency)); // jet finding efficiency (optional)
676 if(!priorLocal) {
677 printf(" > couldn't find prior ! < \n");
678 return 0x0;
679 } else printf(" 1) retrieved prior \n");
680 // switch back to root dir of this unfolding procedure
681 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
682
683 // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
684 TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
685 // unfolded local will be filled with the result of the unfolding
686 TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
687
688 // full response matrix and kinematic efficiency
689 TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
690 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
691
692 // step 2) start the unfolding
693 Int_t status(-1), i(0);
694 while(status < 0 && i < 100) {
695 // i > 0 means that the first iteration didn't converge. in that case, the result of the first
696 // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
697 if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
698 status = AliUnfolding::Unfold(
699 resizedResponseLocal, // response matrix
700 kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
701 measuredJetSpectrumLocal, // measured spectrum
702 priorLocal, // initial conditions (set NULL to use measured spectrum)
703 unfoldedLocal); // results
704 // status holds the minuit fit status (where 0 means convergence)
705 i++;
706 }
707 // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
708 TH2D* hPearson(0x0);
709 if(status == 0 && gMinuit->fISW[1] == 3) {
710 // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
711 TVirtualFitter *fitter(TVirtualFitter::GetFitter());
712 if(gMinuit) gMinuit->Command("SET COV");
713 TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
714 TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
715 pearson->Print();
716 hPearson= new TH2D(*pearson);
717 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
718 hPearson = ProtectHeap(hPearson);
719 hPearson->Write();
720 } else status = -1;
721 if(fMergeBinsArray) unfoldedLocal = MergeSpectrumBins(fMergeBinsArray, unfoldedLocal, hPearson);
722
723 // step 3) refold the unfolded spectrum and save the ratio measured / refolded
724 TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
725 foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
726 unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
727 TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
728 if(ratio) {
729 ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
730 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
731 ratio = ProtectHeap(ratio);
732 ratio->Write();
733 }
734
735 // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
736 // objects are cloned using 'ProtectHeap()'
737 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
738 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
739 measuredJetSpectrumLocal->Write();
740
741 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
742 resizedResponseLocal->Write();
743
744 unfoldedLocal = ProtectHeap(unfoldedLocal);
745 if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
746 unfoldedLocal->Write();
747
748 foldedLocal = ProtectHeap(foldedLocal);
749 foldedLocal->Write();
750
751 priorLocal = ProtectHeap(priorLocal);
752 priorLocal->Write();
753
754 // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
755 TH1F* fitStatus(new TH1F(Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), 4, -0.5, 3.5));
756 fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
757 fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
758 fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
759 fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
760 fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
761 fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
762 fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
763 fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
764 fitStatus->Write();
765
766 return unfoldedLocal;
767}
768//_____________________________________________________________________________
769TH1D* AliJetFlowTools::UnfoldSpectrumBayesian(
770 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
771 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
772 const TH1D* kinematicEfficiency, // kinematic efficiency
773 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
774 const TString suffix, // suffix (in, out)
775 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
776{
777 // use bayesian unfolding from the RooUnfold package to unfold jet spectra
778
779 // 1) get a prior for unfolding.
780 TH1D* priorLocal( GetPrior(
781 measuredJetSpectrum, // jet pt in pt rec bins
782 resizedResponse, // full response matrix, normalized in slides of pt true
783 kinematicEfficiency, // kinematic efficiency
784 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
785 suffix, // suffix (in, out)
786 jetFindingEfficiency)); // jet finding efficiency (optional)
787 if(!priorLocal) {
788 printf(" > couldn't find prior ! < \n");
789 return 0x0;
790 } else printf(" 1) retrieved prior \n");
791 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
792
793 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
794 // measured jets in pt rec binning
795 TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
796 // local copie of the response matrix
797 TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
798 // local copy of response matrix, all true slides normalized to 1
799 // this response matrix will eventually be used in the re-folding routine
800 TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
801 resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
802 // kinematic efficiency
803 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
804 // place holder histos
805 TH1D *unfoldedLocalBayes(0x0);
806 TH1D *foldedLocalBayes(0x0);
807 cout << " 2) setup necessary input " << endl;
808 // 4) get transpose matrices
809 // a) get the transpose of the full response matrix
810 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
811 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
812 // normalize it with the prior. this will ensure that high statistics bins will constrain the
813 // end result most strenuously than bins with limited number of counts
814 responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
815 // 3) get response for Bayesian unfolding
816 RooUnfoldResponse responseBayes(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedBayes_%s", suffix.Data()), Form("respCombinedBayes_%s", suffix.Data()));
817
818 // 4) actualy unfolding loop
819 RooUnfoldBayes unfoldBayes(&responseBayes, measuredJetSpectrumLocal, (!strcmp("in", suffix.Data())) ? fBayesianIterIn : fBayesianIterOut);
820 RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
821 unfoldedLocalBayes = (TH1D*)unfoldBayes.Hreco(errorTreatment);
822 // correct the spectrum for the kinematic efficiency
823 unfoldedLocalBayes->Divide(kinematicEfficiencyLocal);
824 // get the pearson coefficients from the covariance matrix
825 TMatrixD covarianceMatrix = unfoldBayes.Ereco(errorTreatment);
826 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
827 TH2D* hPearson(0x0);
828 if(pearson) {
829 hPearson = new TH2D(*pearson);
830 pearson->Print();
831 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
832 hPearson = ProtectHeap(hPearson);
833 hPearson->Write();
834 if(fMergeBinsArray) unfoldedLocalBayes = MergeSpectrumBins(fMergeBinsArray, unfoldedLocalBayes, hPearson);
835 } else return 0x0; // return if unfolding didn't converge
836
837 // 5) refold the unfolded spectrum
838 foldedLocalBayes = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalBayes, resizedResponseLocalNorm, kinematicEfficiencyLocal);
839 TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalBayes, "ratio measured / re-folded", kTRUE));
840 ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
841 ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
842 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
843 ratio->Write();
844 cout << " 7) refolded the unfolded spectrum " << endl;
845
846 // write the measured, unfolded and re-folded spectra to the output directory
847 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
848 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
849 measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
850 measuredJetSpectrumLocal->Write(); // input spectrum
851 unfoldedLocalBayes->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
852 unfoldedLocalBayes = ProtectHeap(unfoldedLocalBayes);
853 if(jetFindingEfficiency) unfoldedLocalBayes->Divide(jetFindingEfficiency);
854 unfoldedLocalBayes->Write(); // unfolded spectrum
855 foldedLocalBayes->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
856 foldedLocalBayes = ProtectHeap(foldedLocalBayes);
857 foldedLocalBayes->Write(); // re-folded spectrum
858
859 // save more general bookkeeeping histograms to the output directory
860 responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
861 responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
862 responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
863 responseMatrixLocalTransposePrior->Write();
864 priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
865 priorLocal->SetXTitle("p_{t} [GeV/c]");
866 priorLocal = ProtectHeap(priorLocal);
867 priorLocal->Write();
868 resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
869 resizedResponseLocalNorm->Write();
870
871 // save some info
872 TH1F* fitStatus(new TH1F(Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), 1, -0.5, 0.5));
873 fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fBayesianIterIn : fBayesianIterOut);
874 fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fBayesianIterIn" : "fBayesianIterOut");
875 fitStatus->Write();
876
877 return unfoldedLocalBayes;
878}
879//_____________________________________________________________________________
880TH1D* AliJetFlowTools::FoldSpectrum(
881 const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
882 const TH2D* resizedResponse, // response matrix
883 const TH1D* kinematicEfficiency, // kinematic efficiency
884 const TH1D* measuredJetSpectrumTrueBins, // unfolding template: same binning is pt gen of response
885 const TString suffix, // suffix (in or out of plane)
886 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
887{
888 // simple function to fold the given spectrum with the in-plane and out-of-plane response
889
890 // 0) for consistency with the other methods, keep the same nomenclature which admittedly is a bit confusing
891 // what is 'unfolded' here, is just a clone of the input spectrum, binned to the 'unfolded' binning
892 TH1D* unfoldedLocal((TH1D*)measuredJetSpectrum->Clone(Form("unfoldedLocal_%s", suffix.Data())));
893
894 // 1) full response matrix and kinematic efficiency
895 TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
896 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
897
898 // step 2) fold the 'unfolded' spectrum and save the ratio measured / refolded
899 TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
900 foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
901 unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
902
903 // step 3) write histograms to file. to ensure that these have unique identifiers on the heap,
904 // objects are cloned using 'ProtectHeap()'
905 TH1D* measuredJetSpectrumLocal((TH1D*)(measuredJetSpectrum->Clone("tempObject")));
906 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
907 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
908 measuredJetSpectrumLocal->Write();
909
910 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
911 resizedResponseLocal->Write();
912
913 unfoldedLocal = ProtectHeap(unfoldedLocal);
914 if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
915 unfoldedLocal->Write();
916
917 foldedLocal = ProtectHeap(foldedLocal);
918 foldedLocal->Write();
919
920 // return the folded result
921 return foldedLocal;
922}
923//_____________________________________________________________________________
924Bool_t AliJetFlowTools::PrepareForUnfolding(TH1* customIn, TH1* customOut)
925{
926 // prepare for unfolding
927 if(fRawInputProvided) return kTRUE;
928 if(!fInputList) {
929 printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
930 return kFALSE;
931 }
932 if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
933 // check if the pt bin for true and rec have been set
934 if(!fBinsTrue || !fBinsRec) {
935 printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
936 return kFALSE;
937 }
938 if(!fRMSSpectrumIn && fDphiUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
939 // procedures, these profiles will be nonsensical, user is responsible
940 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
941 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
942 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
943 }
944 if(!fTrainPower) {
945 // clear minuit state to avoid constraining the fit with the results of the previous iteration
946 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
947 }
948 // extract the spectra
949 TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
950 if(fRho0) spectrumName = Form("fHistJetPsi2PtRho0_%i", fCentralityArray->At(0));
951 if(!fInputList->FindObject(spectrumName.Data())) {
952 printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
953 return kFALSE;
954 }
955
956 // get the first scaled spectrum
957 fJetPtDeltaPhi = (TH2D*)fInputList->FindObject(spectrumName.Data());
958 // clone the spectrum on the heap. this is necessary since scale or add change the
959 // contents of the original histogram
960 fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
961 fJetPtDeltaPhi->Scale(fCentralityWeights->At(0));
962 printf("Extracted %s wight weight %.2f \n", spectrumName.Data(), fCentralityWeights->At(0));
963 // merge subsequent bins (if any)
964 for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
965 spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i));
966 printf( " Merging with %s with weight %.4f \n", spectrumName.Data(), fCentralityWeights->At(i));
967 fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())), fCentralityWeights->At(i));
968 }
969
970 // in plane spectrum
971 if(!fDphiUnfolding) {
972 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e");
973 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40, "e");
974 } else {
975 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10, "e");
976 fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40, "e"));
977 fSpectrumIn = ProtectHeap(fSpectrumIn);
978 // out of plane spectrum
979 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30, "e");
980 fSpectrumOut = ProtectHeap(fSpectrumOut);
981 }
982 // if a custom input is passed, overwrite existing histograms
983 if(customIn) fSpectrumIn = dynamic_cast<TH1D*>(customIn);
984 if(customOut) fSpectrumOut = dynamic_cast<TH1D*>(customOut);
985
986 // normalize spectra to event count if requested
987 if(fNormalizeSpectra) {
988 TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(0))));
989 rho->Scale(fCentralityWeights->At(0));
990 for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
991 rho->Add((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(i))), fCentralityWeights->At(i));
992 }
993 if(!rho) return 0x0;
994 Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
995 if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
996 if(fEventCount > 0) {
997 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
998 fSpectrumOut->Sumw2();
999 Double_t pt(0);
1000 Double_t error(0); // lots of issues with the errors here ...
1001 for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
1002 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount; // normalized count
1003 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
1004 fSpectrumIn->SetBinContent(1+i, pt);
1005 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
1006 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
1007 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
1008 }
1009 for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
1010 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount; // normalized count
1011 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
1012 fSpectrumOut->SetBinContent(1+i, pt);
1013 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
1014 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
1015 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
1016 }
1017 }
1018 if(normalizeToFullSpectrum) fEventCount = -1;
1019 }
1020 // extract the delta pt matrices
1021 TString deltaptName("");
1022 if(!fRho0) {
1023 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
1024 } else {
1025 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJRho0_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2Rho0_%i", fCentralityArray->At(0));
1026 }
1027 fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1028 if(!fDeltaPtDeltaPhi) {
1029 printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
1030 printf(" > may be ok, depending no what you want to do < \n");
1031 fRefreshInput = kTRUE;
1032 return kTRUE;
1033 }
1034
1035 // clone the distribution on the heap and if requested merge with other centralities
1036 fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
1037 fDeltaPtDeltaPhi->Scale(fCentralityWeights->At(0));
1038 printf("Extracted %s with weight %.2f \n", deltaptName.Data(), fCentralityWeights->At(0));
1039 for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1040 deltaptName = (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i));
1041 printf(" Merging with %s with weight %.4f \n", deltaptName.Data(), fCentralityWeights->At(i));
1042 fDeltaPtDeltaPhi->Add((TH2D*)fInputList->FindObject(deltaptName.Data()), fCentralityWeights->At(i));
1043 }
1044
1045 // in plane delta pt distribution
1046 if(!fDphiUnfolding) {
1047 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e");
1048 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40, "e");
1049 } else {
1050 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10, "e");
1051 fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40, "e"));
1052 // out of plane delta pt distribution
1053 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30, "e");
1054 fDptInDist = ProtectHeap(fDptInDist);
1055 fDptOutDist = ProtectHeap(fDptOutDist);
1056 // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
1057 }
1058
1059 // create a rec - true smeared response matrix
1060 TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1061 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
1062 Bool_t skip = kFALSE;
1063 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
1064 (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1065 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1066 }
1067 }
1068 TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
1069 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
1070 Bool_t skip = kFALSE;
1071 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
1072 (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
1073 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1074 }
1075 }
1076 fDptIn = new TH2D(*rfIn);
1077 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
1078 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1079 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1080 fDptIn = ProtectHeap(fDptIn);
1081 fDptOut = new TH2D(*rfOut);
1082 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
1083 fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1084 fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1085 fDptOut = ProtectHeap(fDptOut);
1086
1087 fRefreshInput = kTRUE; // force cloning of the input
1088 return kTRUE;
1089}
1090//_____________________________________________________________________________
1091Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) {
1092 // prepare for unfoldingUA - more robust method to extract input spectra from file
1093 // will replace PrepareForUnfolding eventually (09012014)
1094 if(!fInputList) {
1095 printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
1096 return kFALSE;
1097 }
1098 if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
1099 // check if the pt bin for true and rec have been set
1100 if(!fBinsTrue || !fBinsRec) {
1101 printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
1102 return kFALSE;
1103 }
1104 if(!fTrainPower) {
1105 // clear minuit state to avoid constraining the fit with the results of the previous iteration
1106 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
1107 }
1108 // extract the spectra
1109 TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
1110 if(fRho0) spectrumName = Form("fHistJetPsi2PtRho0_%i", fCentralityArray->At(0));
1111 fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
1112 if(!fJetPtDeltaPhi) {
1113 printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
1114 return kFALSE;
1115 }
1116 if(fCentralityArray) {
1117 for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1118 spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i));
1119 fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())));
1120 }
1121 }
1122 fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
1123 // in plane spectrum
1124 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e");
1125 // extract the delta pt matrices
1126 TString deltaptName("");
1127 if(!fRho0) {
1128 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
1129 } else {
1130 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJRho0_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2Rho0_%i", fCentralityArray->At(0));
1131 }
1132 fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1133 if(!fDeltaPtDeltaPhi) {
1134 printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
1135 printf(" > may be ok, depending no what you want to do < \n");
1136 fRefreshInput = kTRUE;
1137 return kTRUE;
1138 }
1139 if(fCentralityArray) {
1140 for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1141 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i));
1142 fDeltaPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(deltaptName.Data())));
1143 }
1144 }
1145
1146 fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
1147 // in plane delta pt distribution
1148 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), low, up, "e");
1149 // create a rec - true smeared response matrix
1150 TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1151 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
1152 Bool_t skip = kFALSE;
1153 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
1154 (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1155 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1156 }
1157 }
1158 fDptIn = new TH2D(*rfIn);
1159 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
1160 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1161 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1162 fDptIn = ProtectHeap(fDptIn);
1163
1164 return kTRUE;
1165}
1166//_____________________________________________________________________________
1167TH1D* AliJetFlowTools::GetPrior(
1168 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
1169 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
1170 const TH1D* kinematicEfficiency, // kinematic efficiency
1171 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
1172 const TString suffix, // suffix (in, out)
1173 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
1174{
1175 // 1) get a prior for unfolding.
1176 // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
1177 TH1D* unfolded(0x0);
1178 TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
1179 dirOut->cd();
1180 switch (fPrior) { // select the prior for unfolding
1181 case kPriorTF1 : {
1182 if(!fPriorUser) {
1183 printf("> GetPrior:: FATAL ERROR! TF1 prior requested but prior has not been set ! < \n");
1184 return 0x0;
1185 }
1186 return fPriorUser;
1187 } break;
1188 case kPriorPythia : {
1189 if(!fPriorUser) {
1190 printf("> GetPrior:: FATAL ERROR! pythia prior requested but prior has not been set ! < \n");
1191 return 0x0;
1192 }
1193 // rebin the given prior to the true spectrum (creates a new histo)
1194 return RebinTH1D(fPriorUser, fBinsTrue, Form("kPriorPythia_%s", suffix.Data()), kFALSE);
1195 } break;
1196 case kPriorChi2 : {
1197 TArrayI* placeHolder(0x0);
1198 if(fMergeBinsArray) {
1199 placeHolder = fMergeBinsArray;
1200 fMergeBinsArray = 0x0;
1201 }
1202 if(fBinsTruePrior && fBinsRecPrior) { // if set, use different binning for the prior
1203 TArrayD* tempArrayTrue(fBinsTrue); // temporarily cache the original binning
1204 fBinsTrue = fBinsTruePrior; // switch binning schemes (will be used in UnfoldSpectrumChi2())
1205 TArrayD* tempArrayRec(fBinsRec);
1206 fBinsRec = fBinsRecPrior;
1207 // for the prior, do not re-bin the output
1208 TH1D* measuredJetSpectrumChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
1209 TH1D* measuredJetSpectrumTrueBinsChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
1210 TH2D* resizedResponseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
1211 TH1D* kinematicEfficiencyChi2(resizedResponseChi2->ProjectionX());
1212 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
1213 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
1214 unfolded= UnfoldSpectrumChi2(
1215 measuredJetSpectrumChi2,
1216 resizedResponseChi2,
1217 kinematicEfficiencyChi2,
1218 measuredJetSpectrumTrueBinsChi2, // prior for chi2 unfolding (measured)
1219 TString(Form("prior_%s", suffix.Data())));
1220 if(!unfolded) {
1221 printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1222 printf(" probably Chi2 unfolding did not converge < \n");
1223 if(placeHolder) fMergeBinsArray = placeHolder;
1224 return 0x0;
1225 }
1226 fBinsTrue = tempArrayTrue; // reset bins borders
1227 fBinsRec = tempArrayRec;
1228 // if the chi2 prior has a different binning, rebin to the true binning for the unfolding
1229 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data()))); // rebin unfolded
1230 } else {
1231 unfolded = UnfoldSpectrumChi2(
1232 measuredJetSpectrum,
1233 resizedResponse,
1234 kinematicEfficiency,
1235 measuredJetSpectrumTrueBins, // prior for chi2 unfolding (measured)
1236 TString(Form("prior_%s", suffix.Data())));
1237 if(!unfolded) {
1238 printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1239 printf(" probably Chi2 unfolding did not converge < \n");
1240 if(placeHolder) fMergeBinsArray = placeHolder;
1241 return 0x0;
1242 }
1243 if(placeHolder) fMergeBinsArray = placeHolder;
1244 }
1245 break;
1246
1247 }
1248 case kPriorMeasured : {
1249 unfolded = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("kPriorMeasured_%s", suffix.Data())); // copy template to unfolded to use as prior
1250 }
1251 default : break;
1252 }
1253 // it can be important that the prior is smooth, this can be achieved by
1254 // extrapolating the spectrum with a fitted power law when bins are sparsely filed
1255 if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
1256 TH1D *priorLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
1257 if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
1258 return priorLocal;
1259}
1260//_____________________________________________________________________________
1261TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
1262 // resize the x-axis of a th1d
1263 if(!histo) {
1264 printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1265 return NULL;
1266 }
1267 // see how many bins we need to copy
1268 TH1D* resized = new TH1D(Form("%s_resized_%s", histo->GetName(), suffix.Data()), Form("%s_resized_%s", histo->GetName(), suffix.Data()), up-low, (double)low, (double)up);
1269 // low is the bin number of the first new bin
1270 Int_t l = histo->GetXaxis()->FindBin(low);
1271 // set the values
1272 Double_t _x(0), _xx(0);
1273 for(Int_t i(0); i < up-low; i++) {
1274 _x = histo->GetBinContent(l+i);
1275 _xx=histo->GetBinError(l+i);
1276 resized->SetBinContent(i+1, _x);
1277 resized->SetBinError(i+1, _xx);
1278 }
1279 return resized;
1280}
1281//_____________________________________________________________________________
1282TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
1283 // resize the y-axis of a th2d
1284 if(!histo) {
1285 printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1286 return NULL;
1287 }
1288 // see how many bins we need to copy
1289 TH2D* resized = new TH2D(Form("%s_resized_%s", histo->GetName(), suffix.Data()), Form("%s_resized_%s", histo->GetName(), suffix.Data()), x->GetSize()-1, x->GetArray(), y->GetSize()-1, y->GetArray());
1290 // assume only the y-axis has changed
1291 // low is the bin number of the first new bin
1292 Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1293 // set the values
1294 Double_t _x(0), _xx(0);
1295 for(Int_t i(0); i < x->GetSize(); i++) {
1296 for(Int_t j(0); j < y->GetSize(); j++) {
1297 _x = histo->GetBinContent(i, low+j);
1298 _xx=histo->GetBinError(i, low+1+j);
1299 resized->SetBinContent(i, j, _x);
1300 resized->SetBinError(i, j, _xx);
1301 }
1302 }
1303 return resized;
1304}
1305//_____________________________________________________________________________
1306TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo, Bool_t noError) {
1307 // general method to normalize all vertical slices of a th2 to unity
1308 // i.e. get a probability matrix
1309 if(!histo) {
1310 printf(" > NormalizeTH2D:: NULL pointer passed, returning NULL < \n");
1311 return NULL;
1312 }
1313 Int_t binsX = histo->GetXaxis()->GetNbins();
1314 Int_t binsY = histo->GetYaxis()->GetNbins();
1315
1316 // normalize all slices in x
1317 for(Int_t i(0); i < binsX; i++) { // for each vertical slice
1318 Double_t weight = 0;
1319 for(Int_t j(0); j < binsY; j++) { // loop over all the horizontal components
1320 weight+=histo->GetBinContent(i+1, j+1);
1321 } // now we know the total weight
1322 for(Int_t j(0); j < binsY; j++) {
1323 if (weight <= 0 ) continue;
1324 histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
1325 if(noError) histo->SetBinError( 1+i, j+1, 0.);
1326 else histo->SetBinError( 1+i, j+1, histo->GetBinError( 1+i, j+1)/weight);
1327 }
1328 }
1329 return histo;
1330}
1331//_____________________________________________________________________________
1332TH1* AliJetFlowTools::Bootstrap(TH1* hist, Bool_t kill) {
1333 // resample a TH1
1334 // the returned histogram is new, the original is deleted from the heap if kill is true
1335 if(!hist) {
1336 printf(" > Bootstrap:: fatal error,NULL pointer passed! \n");
1337 return 0x0;
1338 }
1339 else printf(" > Bootstrap:: resampling, may take some time \n");
1340 // clone input histo
1341 TH1* bootstrapped((TH1*)(hist->Clone(Form("%s_bootstrapped", hist->GetName()))));
1342 bootstrapped->Reset();
1343
1344 /* OLD method - slightly underestimates fluctuations
1345 // reset the content
1346 bootstrapped->Reset();
1347 // resample the input histogram
1348 for(Int_t i(0); i < hist->GetEntries(); i++) bootstrapped->Fill(hist->GetRandom()); */
1349
1350 // new method
1351 Double_t mean(0), sigma(0);
1352 Int_t sampledMean(0), entries(0);
1353
1354 for(Int_t i(0); i < hist->GetXaxis()->GetNbins(); i++) {
1355 // for each bin, get the value
1356 mean = hist->GetBinContent(i+1);
1357 sigma = hist->GetBinError(i+1);
1358 // draw a new mean
1359 sampledMean = TMath::Nint(gRandom->Gaus(mean, sigma));
1360 printf(" sampled %i from original number %.2f \n", sampledMean, mean);
1361 // set the new bin content
1362 bootstrapped->SetBinContent(i+1, sampledMean);
1363 if(sampledMean > 0) bootstrapped->SetBinError(i+1, TMath::Sqrt(sampledMean));
1364 entries += sampledMean;
1365 }
1366 printf(" Done bootstrapping, setting number of entries to %i \n", entries);
1367 bootstrapped->SetEntries((double)entries);
1368
1369
1370 // if requested kill input histo
1371 if(kill) delete hist;
1372 // return resampled histogram
1373 return bootstrapped;
1374}
1375//_____________________________________________________________________________
1376TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
1377 // return a TH1D with the supplied histogram rebinned to the supplied bins
1378 // the returned histogram is new, the original is deleted from the heap if kill is true
1379 if(!histo || !bins) {
1380 printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
1381 return NULL;
1382 }
1383 // create the output histo
1384 TString name = histo->GetName();
1385 name+="_template";
1386 name+=suffix;
1387 TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1388 for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
1389 // loop over the bins of the old histo and fill the new one with its data
1390 rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
1391 }
1392 if(kill) delete histo;
1393 return rebinned;
1394}
1395//_____________________________________________________________________________
1396TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1397 // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
1398 // not static as it is just a wrapper for the response maker object
1399 if(!fResponseMaker || !binsTrue || !binsRec) {
1400 printf(" > RebinTH2D:: function called with NULL arguments < \n");
1401 return 0x0;
1402 }
1403 TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1404 return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
1405}
1406//_____________________________________________________________________________
1407TH2D* AliJetFlowTools::MatrixMultiplication(TH2D* a, TH2D* b, TString name)
1408{
1409 // multiply two matrices
1410 if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1411 TH2D* c = (TH2D*)a->Clone("c");
1412 for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1413 for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1414 Double_t val = 0;
1415 for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1416 Int_t y2 = x1;
1417 val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1418 }
1419 c->SetBinContent(x2, y1, val);
1420 c->SetBinError(x2, y1, 0.);
1421 }
1422 }
1423 if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1424 return c;
1425}
1426//_____________________________________________________________________________
1427TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale)
1428{
1429 // normalize a th1d to a certain scale
1430 histo->Sumw2();
1431 Double_t integral = histo->Integral()*scale;
1432 if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1433 else if (scale != 1.) histo->Scale(1./scale, "width");
1434 else printf(" > Histogram integral < 0, cannot normalize \n");
1435 return histo;
1436}
1437//_____________________________________________________________________________
1438TH1D* AliJetFlowTools::MergeSpectrumBins(TArrayI* bins, TH1D* spectrum, TH2D* corr)
1439{
1440 // merge a spectrum histogram taking into account the correlation terms
1441 if(!(bins&&spectrum)) {
1442 printf(" > NULL pointer passed as argument in MergeSpectrumBins ! < \n");
1443 return 0x0;
1444 }
1445 Double_t sum(0), error(0), pearson(0);
1446 // take the sum of the bin content
1447 sum += spectrum->GetBinContent(bins->At(0));
1448 sum += spectrum->GetBinContent(bins->At(1));
1449 // quadratically sum the errors
1450 error += TMath::Power(spectrum->GetBinError(bins->At(0)), 2);
1451 error += TMath::Power(spectrum->GetBinError(bins->At(1)), 2);
1452 // add the covariance term
1453 pearson = corr->GetBinContent(bins->At(0), bins->At(1));
1454 if(!corr) {
1455 printf(" > PANIC ! something is wrong with the covariance matrix, assuming full correlation ! < \n ");
1456 pearson = 1;
1457 }
1458 error += 2.*spectrum->GetBinError(bins->At(0))*spectrum->GetBinError(bins->At(1))*pearson;
1459 spectrum->SetBinContent(1, sum);
1460 if(error <= 0) return spectrum;
1461 spectrum->SetBinError(1, TMath::Sqrt(error));
1462 printf(" > sum is %.2f \t error is %.8f < \n", sum, error);
1463 printf(" > REPLACING BIN CONTENT OF FIRST BIN (%i) WITH MERGED BINS (%i) and (%i) < \n", 1, bins->At(0), bins->At(1));
1464 return spectrum;
1465}
1466//_____________________________________________________________________________
1467TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
1468{
1469 // Calculate pearson coefficients from covariance matrix
1470 TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1471 Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
1472 Double_t pearson(0.);
1473 if(nrows==0 && ncols==0) return 0x0;
1474 for(Int_t row = 0; row < nrows; row++) {
1475 for(Int_t col = 0; col<ncols; col++) {
1476 if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1477 (*pearsonCoefficients)(row,col) = pearson;
1478 }
1479 }
1480 return pearsonCoefficients;
1481}
1482//_____________________________________________________________________________
1483TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
1484 // smoothen the spectrum using a user defined function
1485 // returns a clone of the original spectrum if fitting failed
1486 // if kill is kTRUE the input spectrum will be deleted from the heap
1487 // if 'count' is selected, bins are filled with integers (necessary if the
1488 // histogram is interpreted in a routine which accepts only counts)
1489 if(!spectrum || !function) return 0x0;
1490 if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
1491 TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1492 temp->Sumw2(); // if already called on the original, this will give off a warning but do nothing
1493 TFitResultPtr r = temp->Fit(function, "", "", min, max);
1494 if((int)r == 0) { // MINUIT status
1495 for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1496 if(temp->GetBinCenter(i) > start) { // from this pt value use extrapolation
1497 (counts) ? temp->SetBinContent(i, (int)(function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i))) : temp->SetBinContent(i, function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
1498 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1499 }
1500 }
1501 }
1502 if(kill) delete spectrum;
1503 return temp;
1504}
1505//_____________________________________________________________________________
1506void AliJetFlowTools::Style(Bool_t legacy)
1507{
1508 // set global style for your current aliroot session
1509 if(!gStyle) return;
1510 // legacy style is pleasing to the eye, default is the formal ALICE style
1511 if(legacy) {
1512 gStyle->SetCanvasColor(-1);
1513 gStyle->SetPadColor(-1);
1514 gStyle->SetFrameFillColor(-1);
1515 gStyle->SetHistFillColor(-1);
1516 gStyle->SetTitleFillColor(-1);
1517 gStyle->SetFillColor(-1);
1518 gStyle->SetFillStyle(4000);
1519 gStyle->SetStatStyle(0);
1520 gStyle->SetTitleStyle(0);
1521 gStyle->SetCanvasBorderSize(0);
1522 gStyle->SetFrameBorderSize(0);
1523 gStyle->SetLegendBorderSize(0);
1524 gStyle->SetStatBorderSize(0);
1525 gStyle->SetTitleBorderSize(0);
1526 } else {
1527 gStyle->Reset("Plain");
1528 gStyle->SetOptTitle(0);
1529 gStyle->SetOptStat(0);
1530 gStyle->SetPalette(1);
1531 gStyle->SetCanvasColor(10);
1532 gStyle->SetCanvasBorderMode(0);
1533 gStyle->SetFrameLineWidth(1);
1534 gStyle->SetFrameFillColor(kWhite);
1535 gStyle->SetPadColor(10);
1536 gStyle->SetPadTickX(1);
1537 gStyle->SetPadTickY(1);
1538 gStyle->SetPadBottomMargin(0.15);
1539 gStyle->SetPadLeftMargin(0.15);
1540 gStyle->SetHistLineWidth(1);
1541 gStyle->SetHistLineColor(kRed);
1542 gStyle->SetFuncWidth(2);
1543 gStyle->SetFuncColor(kGreen);
1544 gStyle->SetLineWidth(2);
1545 gStyle->SetLabelSize(0.045,"xyz");
1546 gStyle->SetLabelOffset(0.01,"y");
1547 gStyle->SetLabelOffset(0.01,"x");
1548 gStyle->SetLabelColor(kBlack,"xyz");
1549 gStyle->SetTitleSize(0.05,"xyz");
1550 gStyle->SetTitleOffset(1.25,"y");
1551 gStyle->SetTitleOffset(1.2,"x");
1552 gStyle->SetTitleFillColor(kWhite);
1553 gStyle->SetTextSizePixels(26);
1554 gStyle->SetTextFont(42);
1555 gStyle->SetLegendBorderSize(0);
1556 gStyle->SetLegendFillColor(kWhite);
1557 gStyle->SetLegendFont(42);
1558 }
1559}
1560//_____________________________________________________________________________
1561void AliJetFlowTools::Style(TCanvas* c, TString style)
1562{
1563 // set a default style for a canvas
1564 if(!strcmp(style.Data(), "PEARSON")) {
1565 printf(" > style PEARSON canvas < \n");
1566 gStyle->SetOptStat(0);
1567 c->SetGridx();
1568 c->SetGridy();
1569 c->SetTicks();
1570 return;
1571 } else if(!strcmp(style.Data(), "SPECTRUM")) {
1572 printf(" > style SPECTRUM canvas < \n");
1573 gStyle->SetOptStat(0);
1574 c->SetLogy();
1575 c->SetGridx();
1576 c->SetGridy();
1577 c->SetTicks();
1578 return;
1579 } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1580}
1581//_____________________________________________________________________________
1582void AliJetFlowTools::Style(TVirtualPad* c, TString style, Bool_t legacy)
1583{
1584 // set a default style for a canva
1585
1586 if(legacy) {
1587 c->SetLeftMargin(.25);
1588 c->SetBottomMargin(.25);
1589 }
1590 else Style();
1591 if(!strcmp(style.Data(), "PEARSON")) {
1592 printf(" > style PEARSON pad < \n");
1593 gStyle->SetOptStat(0);
1594 c->SetGridx();
1595 c->SetGridy();
1596 c->SetTicks();
1597 return;
1598 } else if(!strcmp(style.Data(), "SPECTRUM")) {
1599 printf(" > style SPECTRUM pad < \n");
1600 gStyle->SetOptStat(0);
1601 c->SetLogy();
1602 c->SetGridx();
1603 c->SetGridy();
1604 c->SetTicks();
1605 return;
1606 } else if (!strcmp(style.Data(), "GRID")) {
1607 printf(" > style GRID pad < \n");
1608 gStyle->SetOptStat(0);
1609 c->SetGridx();
1610 c->SetGridy();
1611 c->SetTicks();
1612 } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1613}
1614//_____________________________________________________________________________
1615void AliJetFlowTools::Style(TLegend* l)
1616{
1617 // set a default style for a legend
1618 l->SetFillColor(0);
1619 l->SetBorderSize(0);
1620 if(gStyle) l->SetTextSize(gStyle->GetTextSize()*.08);
1621}
1622//_____________________________________________________________________________
1623void AliJetFlowTools::Style(TH1* h, EColor col, histoType type, Bool_t legacy)
1624{
1625 // style a histo
1626 h->GetYaxis()->SetNdivisions(505);
1627 h->SetLineColor(col);
1628 h->SetMarkerColor(col);
1629 h->SetLineWidth(2);
1630 h->SetMarkerSize(1);
1631 if(legacy) {
1632 h->SetTitle("");
1633 h->GetYaxis()->SetLabelSize(0.05);
1634 h->GetXaxis()->SetLabelSize(0.05);
1635 h->GetYaxis()->SetTitleOffset(1.5);
1636 h->GetXaxis()->SetTitleOffset(1.5);
1637 h->GetYaxis()->SetTitleSize(.05);
1638 h->GetXaxis()->SetTitleSize(.05);
1639 } else Style();
1640 switch (type) {
1641 case kInPlaneSpectrum : {
1642 h->SetTitle("IN PLANE");
1643 h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1644 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1645 } break;
1646 case kOutPlaneSpectrum : {
1647 h->SetTitle("OUT OF PLANE");
1648 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1649 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1650 } break;
1651 case kUnfoldedSpectrum : {
1652 h->SetTitle("UNFOLDED");
1653 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1654 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1655 } break;
1656 case kFoldedSpectrum : {
1657 h->SetTitle("FOLDED");
1658 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1659 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1660 } break;
1661 case kMeasuredSpectrum : {
1662 h->SetTitle("MEASURED");
1663 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1664 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{cht} (GeV/#it{c})");
1665 } break;
1666 case kBar : {
1667 h->SetFillColor(col);
1668 h->SetBarWidth(.6);
1669 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1670 h->SetBarOffset(0.2);
1671 }
1672 case kRatio : {
1673 h->SetMarkerStyle(8);
1674 h->SetMarkerSize(1);
1675 } break;
1676 case kDeltaPhi : {
1677 h->GetYaxis()->SetTitle("[counts]");
1678 h->GetXaxis()->SetTitle("#Delta #phi");
1679 }
1680 default : break;
1681 }
1682}
1683//_____________________________________________________________________________
1684void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type, Bool_t legacy)
1685{
1686 // style a tgraph
1687 h->GetYaxis()->SetNdivisions(505);
1688 h->SetLineColor(col);
1689 h->SetMarkerColor(col);
1690 h->SetLineWidth(2);
1691 h->SetMarkerSize(1);
1692 h->SetTitle("");
1693 h->SetFillColor(kCyan);
1694 if(legacy) {
1695 h->GetYaxis()->SetLabelSize(0.05);
1696 h->GetXaxis()->SetLabelSize(0.05);
1697 h->GetYaxis()->SetTitleOffset(1.6);
1698 h->GetXaxis()->SetTitleOffset(1.6);
1699 h->GetYaxis()->SetTitleSize(.05);
1700 h->GetXaxis()->SetTitleSize(.05);
1701 } else Style();
1702 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1703 switch (type) {
1704 case kInPlaneSpectrum : {
1705 h->SetTitle("IN PLANE");
1706 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1707 } break;
1708 case kOutPlaneSpectrum : {
1709 h->SetTitle("OUT OF PLANE");
1710 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1711 } break;
1712 case kUnfoldedSpectrum : {
1713 h->SetTitle("UNFOLDED");
1714 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1715 } break;
1716 case kFoldedSpectrum : {
1717 h->SetTitle("FOLDED");
1718 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1719 } break;
1720 case kMeasuredSpectrum : {
1721 h->SetTitle("MEASURED");
1722 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1723 } break;
1724 case kRatio : {
1725// h->GetYaxis()->SetTitle("#frac{d#it{N_{in plane}^{jet}}}{d#it{p}_{T}} / #frac{d#it{N_{out of plane}^{jet}}}{d#it{p}_{T}}");
1726 h->GetYaxis()->SetTitle("(d#it{N}^{ch, jet}_{in plane}/(d#it{p}_{T}d#eta))/(d#it{N}^{ch,jet}_{out of plane}/(d#it{p}_{T}d#eta))");
1727 } break;
1728 case kV2 : {
1729// h->GetYaxis()->SetTitle("#it{v}_{2} = #frac{1}{#it{R}} #frac{#pi}{4} #frac{#it{N_{in plane}} - #it{N_{out of plane}}}{#it{N_{in plane}} + #it{N_{out of plane}}}");
1730 h->GetYaxis()->SetTitle("#it{v}_{2}^{ch, jet} \{EP, |#Delta#eta|>0.9 \} ");
1731 h->GetYaxis()->SetRangeUser(-.5, 1.);
1732 h->SetMarkerStyle(8);
1733 h->SetMarkerSize(1);
1734 } break;
1735 default : break;
1736 }
1737}
1738//_____________________________________________________________________________
1739void AliJetFlowTools::GetNominalValues(
1740 TH1D*& ratio, // pointer reference, output of this function
1741 TGraphErrors*& v2, // pointer reference, as output of this function
1742 TArrayI* in,
1743 TArrayI* out,
1744 TString inFile,
1745 TString outFile) const
1746{
1747 // pass clones of the nominal points and nominal v2 values
1748 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close(); // if for some weird reason the unfolding output is still mutable
1749 TFile* readMe(new TFile(inFile.Data(), "READ")); // open unfolding output read-only
1750 if(readMe->IsZombie()) {
1751 printf(" > Fatal error, couldn't read %s for post processing ! < \n", inFile.Data());
1752 return;
1753 }
1754 printf("\n\n\n\t\t GetNominalValues \n > Recovered the following file structure : \n <");
1755 readMe->ls();
1756 TFile* output(new TFile(outFile.Data(), "RECREATE")); // create new output file
1757 // get some placeholders, have to be initialized but will be deleted
1758 ratio = new TH1D("nominal", "nominal", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1759 TH1D* nominalIn(new TH1D("nominal in", "nominal in", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1760 TH1D* nominalOut(new TH1D("nominal out", "nominal out", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1761 Int_t iIn[] = {in->At(0), in->At(0)}; // first value is the nominal point
1762 Int_t iOut[] = {out->At(0), out->At(0)};
1763
1764 // call the functions and set the relevant pointer references
1765 TH1D* dud(0x0);
1766 DoIntermediateSystematics(
1767 new TArrayI(2, iIn),
1768 new TArrayI(2, iOut),
1769 dud, dud, dud, dud, dud, dud,
1770 ratio, // pointer reference, output of this function
1771 nominalIn,
1772 nominalOut,
1773 1,
1774 fBinsTrue->At(0),
1775 fBinsTrue->At(fBinsTrue->GetSize()-1),
1776 readMe,
1777 "nominal_values");
1778 v2 = GetV2(nominalIn, nominalOut, fEventPlaneRes, "nominal v_{2}");
1779
1780 // close the open files, reclaim ownership of histograms which are necessary outside of the file scope
1781 ratio->SetDirectory(0); // disassociate from current gDirectory
1782 readMe->Close();
1783}
1784//_____________________________________________________________________________
1785void AliJetFlowTools::GetCorrelatedUncertainty(
1786 TGraphAsymmErrors*& corrRatio, // correlated uncertainty pointer reference
1787 TGraphAsymmErrors*& corrV2, // correlated uncertainty pointer reference
1788 TArrayI* variationsIn, // variantions in plane
1789 TArrayI* variationsOut, // variantions out of plane
1790 Bool_t sym, // treat as symmmetric
1791 TArrayI* variations2ndIn, // second source of variations
1792 TArrayI* variations2ndOut, // second source of variations
1793 Bool_t sym2nd, // treat as symmetric
1794 TString type, // type of uncertaitny
1795 TString type2,
1796 Int_t columns, // divide the output canvasses in this many columns
1797 Float_t rangeLow, // lower pt range
1798 Float_t rangeUp, // upper pt range
1799 Float_t corr, // correlation strength
1800 TString in, // input file name (created by this unfolding class)
1801 TString out // output file name (which will hold results of the systematic test)
1802 ) const
1803{
1804 // do full systematics
1805 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close(); // if for some weird reason the unfolding output is still mutable
1806 TFile* readMe(new TFile(in.Data(), "READ")); // open unfolding output read-only
1807 if(readMe->IsZombie()) {
1808 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1809 return;
1810 }
1811 printf("\n\n\n\t\t GetCorrelatedUncertainty \n > Recovered the following file structure : \n <");
1812 readMe->ls();
1813 TFile* output(new TFile(out.Data(), "RECREATE")); // create new output file
1814
1815 // create some null placeholder pointers
1816 TH1D* relativeErrorVariationInLow(0x0);
1817 TH1D* relativeErrorVariationInUp(0x0);
1818 TH1D* relativeErrorVariationOutLow(0x0);
1819 TH1D* relativeErrorVariationOutUp(0x0);
1820 TH1D* relativeError2ndVariationInLow(0x0);
1821 TH1D* relativeError2ndVariationInUp(0x0);
1822 TH1D* relativeError2ndVariationOutLow(0x0);
1823 TH1D* relativeError2ndVariationOutUp(0x0);
1824 TH1D* relativeStatisticalErrorIn(0x0);
1825 TH1D* relativeStatisticalErrorOut(0x0);
1826 // histo for the nominal ratio in / out
1827 TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1828 TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1829 TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1830
1831 // call the functions
1832 if(variationsIn && variationsOut) {
1833 DoIntermediateSystematics(
1834 variationsIn,
1835 variationsOut,
1836 relativeErrorVariationInUp, // pointer reference
1837 relativeErrorVariationInLow, // pointer reference
1838 relativeErrorVariationOutUp, // pointer reference
1839 relativeErrorVariationOutLow, // pointer reference
1840 relativeStatisticalErrorIn, // pointer reference
1841 relativeStatisticalErrorOut, // pointer reference
1842 nominal,
1843 nominalIn,
1844 nominalOut,
1845 columns,
1846 rangeLow,
1847 rangeUp,
1848 readMe,
1849 type);
1850 if(relativeErrorVariationInUp) {
1851 // canvas with the error from variation strength
1852 TCanvas* relativeErrorVariation(new TCanvas(Form("relativeError_%s", type.Data()), Form("relativeError_%s", type.Data())));
1853 relativeErrorVariation->Divide(2);
1854 relativeErrorVariation->cd(1);
1855 Style(gPad, "GRID");
1856 relativeErrorVariationInUp->DrawCopy("b");
1857 relativeErrorVariationInLow->DrawCopy("same b");
1858 Style(AddLegend(gPad));
1859 relativeErrorVariation->cd(2);
1860 Style(gPad, "GRID");
1861 relativeErrorVariationOutUp->DrawCopy("b");
1862 relativeErrorVariationOutLow->DrawCopy("same b");
1863 SavePadToPDF(relativeErrorVariation);
1864 Style(AddLegend(gPad));
1865 relativeErrorVariation->Write();
1866
1867 // now smoothen the diced response error (as it is expected to be flat)
1868 // this is done by fitting a constant to the diced resonse error histo
1869 //
1870 TF1* lin = new TF1("lin", "[0]", rangeLow, rangeUp);
1871 relativeErrorVariationInUp->Fit(lin, "L", "", rangeLow, rangeUp);
1872 if(!gMinuit->fISW[1] == 3) printf(" fit is NOT ok ! " );
1873 for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1874 relativeErrorVariationInUp->SetBinContent(i+1, lin->GetParameter(0));
1875 }
1876 relativeErrorVariationInLow->Fit(lin, "L", "", rangeLow, rangeUp);
1877 printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
1878 for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1879 relativeErrorVariationInLow->SetBinContent(i+1, 0);//lin->GetParameter(0));
1880 }
1881 relativeErrorVariationOutUp->Fit(lin, "L", "", rangeLow, rangeUp);
1882 printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
1883 for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1884 relativeErrorVariationOutUp->SetBinContent(i+1, lin->GetParameter(0));
1885 }
1886 relativeErrorVariationOutLow->Fit(lin, "L", "", rangeLow, rangeUp);
1887 printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
1888 for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1889 relativeErrorVariationOutLow->SetBinContent(i+1, 0);//lin->GetParameter(0));
1890 }
1891 }
1892 }
1893 // call the functions for a second set of systematic sources
1894 if(variations2ndIn && variations2ndOut) {
1895 DoIntermediateSystematics(
1896 variations2ndIn,
1897 variations2ndOut,
1898 relativeError2ndVariationInUp, // pointer reference
1899 relativeError2ndVariationInLow, // pointer reference
1900 relativeError2ndVariationOutUp, // pointer reference
1901 relativeError2ndVariationOutLow, // pointer reference
1902 relativeStatisticalErrorIn, // pointer reference
1903 relativeStatisticalErrorOut, // pointer reference
1904 nominal,
1905 nominalIn,
1906 nominalOut,
1907 columns,
1908 rangeLow,
1909 rangeUp,
1910 readMe,
1911 type2);
1912 if(relativeError2ndVariationInUp) {
1913 // canvas with the error from variation strength
1914 TCanvas* relativeError2ndVariation(new TCanvas(Form("relativeError2nd_%s", type2.Data()), Form("relativeError2nd_%s", type2.Data())));
1915 relativeError2ndVariation->Divide(2);
1916 relativeError2ndVariation->cd(1);
1917 Style(gPad, "GRID");
1918 relativeError2ndVariationInUp->DrawCopy("b");
1919 relativeError2ndVariationInLow->DrawCopy("same b");
1920 Style(AddLegend(gPad));
1921 relativeError2ndVariation->cd(2);
1922 Style(gPad, "GRID");
1923 relativeError2ndVariationOutUp->DrawCopy("b");
1924 relativeError2ndVariationOutLow->DrawCopy("same b");
1925 SavePadToPDF(relativeError2ndVariation);
1926 Style(AddLegend(gPad));
1927 relativeError2ndVariation->Write();
1928 }
1929
1930 }
1931
1932 // and the placeholder for the final systematic
1933 TH1D* relativeErrorInUp(new TH1D("max correlated uncertainty in plane", "max correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1934 TH1D* relativeErrorOutUp(new TH1D("max correlated uncertainty out of plane", "max correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1935 TH1D* relativeErrorInLow(new TH1D("min correlated uncertainty in plane", "min correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1936 TH1D* relativeErrorOutLow(new TH1D("min correlated uncertainty out of plane", "min correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1937 relativeErrorInUp->SetYTitle("relative uncertainty");
1938 relativeErrorOutUp->SetYTitle("relative uncertainty");
1939 relativeErrorInLow->SetYTitle("relative uncertainty");
1940 relativeErrorOutLow->SetYTitle("relative uncertainty");
1941
1942 // merge the correlated errors (FIXME) trivial for one set
1943 Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.);
1944 Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.);
1945 Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.);
1946 Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.);
1947
1948 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
1949 // for the upper bound
1950 if(relativeErrorVariationInUp) aInUp = relativeErrorVariationInUp->GetBinContent(b+1);
1951 if(relativeErrorVariationOutUp) aOutUp = relativeErrorVariationOutUp->GetBinContent(b+1);
1952 if(relativeError2ndVariationInUp) bInUp = relativeError2ndVariationInUp->GetBinContent(b+1);
1953 if(relativeError2ndVariationOutUp) bInLow = relativeError2ndVariationOutUp->GetBinContent(b+1);
1954 dInUp = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp;
1955 // for a symmetric first variation
1956 if(sym) dInUp += aInLow*aInLow;
1957 if(dInUp > 0) relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(dInUp));
1958 dOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp;
1959 if(sym) dOutUp += aOutLow*aOutLow;
1960 if(dOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(dOutUp));
1961 // for the lower bound
1962 if(relativeErrorVariationInLow) aInLow = relativeErrorVariationInLow->GetBinContent(b+1);
1963 if(relativeErrorVariationOutLow) aOutLow = relativeErrorVariationOutLow->GetBinContent(b+1);
1964 if(relativeError2ndVariationInLow) bInLow = relativeError2ndVariationInLow->GetBinContent(b+1);
1965 if(relativeError2ndVariationOutLow) bOutLow = relativeError2ndVariationOutLow->GetBinContent(b+1);
1966 dInLow = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow;
1967 if(sym) dInLow += aInUp*aInUp;
1968 if(dInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1*TMath::Sqrt(dInLow));
1969 dOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow;
1970 if(sym) dOutLow += aOutUp*aOutUp;
1971 if(dOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(dOutLow));
1972 }
1973 // project the estimated errors on the nominal ratio
1974 if(nominal) {
1975 Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
1976 Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
1977 Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
1978 Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
1979 Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
1980 Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
1981 Double_t lowErr(0.), upErr(0.);
1982 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
1983 // add the in and out of plane errors in quadrature
1984 // propagation is tricky because we should consider two cases:
1985 // [1] the error is uncorrelated, and we propagate upper 1 with lower 2 and vice versa
1986 // [2] the error is correlated - in this case upper 1 should be propagated with upper 2
1987 // as the fluctuations are bound to always to of same sign
1988 if(corr <= 0) {
1989 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(1+i)*relativeErrorOutUp->GetBinContent(i+1);
1990 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
1991 } else {
1992 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(1+i)*relativeErrorOutLow->GetBinContent(i+1) -2.*corr*relativeErrorInLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
1993 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1) - 2.*corr*relativeErrorInUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1);
1994 }
1995 ayl[i] = TMath::Sqrt(TMath::Abs(lowErr))*nominal->GetBinContent(i+1);
1996 ayh[i] = TMath::Sqrt(TMath::Abs(upErr))*nominal->GetBinContent(i+1);
1997 // get the bin width (which is the 'error' on x
1998 Double_t binWidth(nominal->GetBinWidth(i+1));
1999 axl[i] = binWidth/2.;
2000 axh[i] = binWidth/2.;
2001 // now get the coordinate for the point
2002 ax[i] = nominal->GetBinCenter(i+1);
2003 ay[i] = nominal->GetBinContent(i+1);
2004 }
2005 // save the nominal ratio
2006 TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
2007 corrRatio = (TGraphAsymmErrors*)nominalError->Clone();
2008 nominalError->SetName("correlated uncertainty");
2009 TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
2010 nominalCanvas->Divide(2);
2011 nominalCanvas->cd(1);
2012 Style(nominal, kBlack);
2013 Style(nominalError, kCyan, kRatio);
2014 nominalError->SetLineColor(kCyan);
2015 nominalError->SetMarkerColor(kCyan);
2016 nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2017 nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
2018 nominalError->DrawClone("a2");
2019 nominal->DrawCopy("same E1");
2020 Style(AddLegend(gPad));
2021 Style(gPad, "GRID");
2022 Style(nominalCanvas);
2023 // save nominal v2 and systematic errors
2024 TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
2025 nominalIn,
2026 nominalOut,
2027 fEventPlaneRes,
2028 "v_{2} with shape uncertainty",
2029 relativeErrorInUp,
2030 relativeErrorInLow,
2031 relativeErrorOutUp,
2032 relativeErrorOutLow,
2033 corr));
2034 // pass the nominal values to the pointer references
2035 corrV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
2036 TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
2037 nominalCanvas->cd(2);
2038 Style(nominalV2, kBlack);
2039 Style(nominalV2Error, kCyan, kV2);
2040 nominalV2Error->SetLineColor(kCyan);
2041 nominalV2Error->SetMarkerColor(kCyan);
2042 nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2043 nominalV2Error->DrawClone("a2");
2044 nominalV2->DrawClone("same E1");
2045 Style(AddLegend(gPad));
2046 Style(gPad, "GRID");
2047 Style(nominalCanvas);
2048 SavePadToPDF(nominalCanvas);
2049 nominalCanvas->Write();
2050 }
2051
2052 TCanvas* relativeError(new TCanvas("relativeCorrelatedError"," relativeCorrelatedError"));
2053 relativeError->Divide(2);
2054 relativeError->cd(1);
2055 Style(gPad, "GRID");
2056 relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2057 Style(relativeErrorInUp, kBlue, kBar);
2058 Style(relativeErrorInLow, kGreen, kBar);
2059 relativeErrorInUp->DrawCopy("b");
2060 relativeErrorInLow->DrawCopy("same b");
2061 Style(relativeStatisticalErrorIn, kRed);
2062 relativeStatisticalErrorIn->DrawCopy("same");
2063 Style(AddLegend(gPad));
2064 relativeError->cd(2);
2065 Style(gPad, "GRID");
2066 relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2067 Style(relativeErrorOutUp, kBlue, kBar);
2068 Style(relativeErrorOutLow, kGreen, kBar);
2069 relativeErrorOutUp->DrawCopy("b");
2070 relativeErrorOutLow->DrawCopy("same b");
2071 Style(relativeStatisticalErrorOut, kRed);
2072 relativeStatisticalErrorOut->DrawCopy("same");
2073 Style(AddLegend(gPad));
2074
2075 // write the buffered file to disk and close the file
2076 SavePadToPDF(relativeError);
2077 relativeError->Write();
2078 output->Write();
2079// output->Close();
2080}
2081//_____________________________________________________________________________
2082void AliJetFlowTools::GetShapeUncertainty(
2083 TGraphAsymmErrors*& shapeRatio, // pointer reference to final shape uncertainty of ratio
2084 TGraphAsymmErrors*& shapeV2, // pointer reference to final shape uncertainty of v2
2085 TArrayI* regularizationIn, // regularization strength systematics, in plane
2086 TArrayI* regularizationOut, // regularization strenght systematics, out of plane
2087 TArrayI* trueBinIn, // true bin ranges, in plane
2088 TArrayI* trueBinOut, // true bin ranges, out of plane
2089 TArrayI* recBinIn, // rec bin ranges, in plane
2090 TArrayI* recBinOut, // rec bin ranges, out of plane
2091 TArrayI* methodIn, // method in
2092 TArrayI* methodOut, // method out
2093 Int_t columns, // divide the output canvasses in this many columns
2094 Float_t rangeLow, // lower pt range
2095 Float_t rangeUp, // upper pt range
2096 Float_t corr, // correlation strength
2097 TString in, // input file name (created by this unfolding class)
2098 TString out // output file name (which will hold results of the systematic test)
2099 ) const
2100{
2101 // do full systematics
2102 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close(); // if for some weird reason the unfolding output is still mutable
2103 TFile* readMe(new TFile(in.Data(), "READ")); // open unfolding output read-only
2104 if(readMe->IsZombie()) {
2105 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
2106 return;
2107 }
2108 printf("\n\n\n\t\t DOSYSTEMATICS \n > Recovered the following file structure : \n <");
2109 readMe->ls();
2110 TFile* output(new TFile(out.Data(), "RECREATE")); // create new output file
2111
2112 // create some null placeholder pointers
2113 TH1D* relativeErrorRegularizationInLow(0x0);
2114 TH1D* relativeErrorRegularizationInUp(0x0);
2115 TH1D* relativeErrorTrueBinInLow(0x0);
2116 TH1D* relativeErrorTrueBinInUp(0x0);
2117 TH1D* relativeErrorRecBinInLow(0x0);
2118 TH1D* relativeErrorRecBinInUp(0x0);
2119 TH1D* relativeErrorMethodInLow(0x0);
2120 TH1D* relativeErrorMethodInUp(0x0);
2121 TH1D* relativeErrorRegularizationOutLow(0x0);
2122 TH1D* relativeErrorRegularizationOutUp(0x0);
2123 TH1D* relativeErrorTrueBinOutLow(0x0);
2124 TH1D* relativeErrorTrueBinOutUp(0x0);
2125 TH1D* relativeErrorRecBinOutLow(0x0);
2126 TH1D* relativeErrorRecBinOutUp(0x0);
2127 TH1D* relativeStatisticalErrorIn(0x0);
2128 TH1D* relativeStatisticalErrorOut(0x0);
2129 TH1D* relativeErrorMethodOutLow(0x0);
2130 TH1D* relativeErrorMethodOutUp(0x0);
2131 // histo for the nominal ratio in / out
2132 TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2133 TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2134 TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2135
2136 // call the functions
2137 if(regularizationIn && regularizationOut) {
2138 DoIntermediateSystematics(
2139 regularizationIn,
2140 regularizationOut,
2141 relativeErrorRegularizationInUp, // pointer reference
2142 relativeErrorRegularizationInLow, // pointer reference
2143 relativeErrorRegularizationOutUp, // pointer reference
2144 relativeErrorRegularizationOutLow, // pointer reference
2145 relativeStatisticalErrorIn, // pointer reference
2146 relativeStatisticalErrorOut, // pointer reference
2147 nominal,
2148 nominalIn,
2149 nominalOut,
2150 columns,
2151 rangeLow,
2152 rangeUp,
2153 readMe,
2154 "regularization",
2155 fRMS);
2156 if(relativeErrorRegularizationInUp) {
2157 // canvas with the error from regularization strength
2158 TCanvas* relativeErrorRegularization(new TCanvas("relativeErrorRegularization", "relativeErrorRegularization"));
2159 relativeErrorRegularization->Divide(2);
2160 relativeErrorRegularization->cd(1);
2161 Style(gPad, "GRID");
2162 relativeErrorRegularizationInUp->DrawCopy("b");
2163 relativeErrorRegularizationInLow->DrawCopy("same b");
2164 Style(AddLegend(gPad));
2165 relativeErrorRegularization->cd(2);
2166 Style(gPad, "GRID");
2167 relativeErrorRegularizationOutUp->DrawCopy("b");
2168 relativeErrorRegularizationOutLow->DrawCopy("same b");
2169 SavePadToPDF(relativeErrorRegularization);
2170 Style(AddLegend(gPad));
2171 relativeErrorRegularization->Write();
2172 }
2173 }
2174 if(trueBinIn && trueBinOut) {
2175 DoIntermediateSystematics(
2176 trueBinIn,
2177 trueBinOut,
2178 relativeErrorTrueBinInUp, // pointer reference
2179 relativeErrorTrueBinInLow, // pointer reference
2180 relativeErrorTrueBinOutUp, // pointer reference
2181 relativeErrorTrueBinOutLow, // pointer reference
2182 relativeStatisticalErrorIn,
2183 relativeStatisticalErrorOut,
2184 nominal,
2185 nominalIn,
2186 nominalOut,
2187 columns,
2188 rangeLow,
2189 rangeUp,
2190 readMe,
2191 "trueBin");
2192 if(relativeErrorTrueBinInUp) {
2193 TCanvas* relativeErrorTrueBin(new TCanvas("relativeErrorTrueBin", "relativeErrorTrueBin"));
2194 relativeErrorTrueBin->Divide(2);
2195 relativeErrorTrueBin->cd(1);
2196 Style(gPad, "GRID");
2197 relativeErrorTrueBinInUp->DrawCopy("b");
2198 relativeErrorTrueBinInLow->DrawCopy("same b");
2199 Style(AddLegend(gPad));
2200 relativeErrorTrueBin->cd(2);
2201 Style(gPad, "GRID");
2202 relativeErrorTrueBinOutUp->DrawCopy("b");
2203 relativeErrorTrueBinOutLow->DrawCopy("same b");
2204 SavePadToPDF(relativeErrorTrueBin);
2205 Style(AddLegend(gPad));
2206 relativeErrorTrueBin->Write();
2207 }
2208 }
2209 if(recBinIn && recBinOut) {
2210 DoIntermediateSystematics(
2211 recBinIn,
2212 recBinOut,
2213 relativeErrorRecBinInUp, // pointer reference
2214 relativeErrorRecBinInLow, // pointer reference
2215 relativeErrorRecBinOutUp, // pointer reference
2216 relativeErrorRecBinOutLow, // pointer reference
2217 relativeStatisticalErrorIn,
2218 relativeStatisticalErrorOut,
2219 nominal,
2220 nominalIn,
2221 nominalOut,
2222 columns,
2223 rangeLow,
2224 rangeUp,
2225 readMe,
2226 "recBin",
2227 fRMS);
2228 if(relativeErrorRecBinOutUp) {
2229 // canvas with the error from regularization strength
2230 TCanvas* relativeErrorRecBin(new TCanvas("relativeErrorRecBin"," relativeErrorRecBin"));
2231 relativeErrorRecBin->Divide(2);
2232 relativeErrorRecBin->cd(1);
2233 Style(gPad, "GRID");
2234 relativeErrorRecBinInUp->DrawCopy("b");
2235 relativeErrorRecBinInLow->DrawCopy("same b");
2236 Style(AddLegend(gPad));
2237 relativeErrorRecBin->cd(2);
2238 Style(gPad, "GRID");
2239 relativeErrorRecBinOutUp->DrawCopy("b");
2240 relativeErrorRecBinOutLow->DrawCopy("same b");
2241 SavePadToPDF(relativeErrorRecBin);
2242 Style(AddLegend(gPad));
2243 relativeErrorRecBin->Write();
2244 }
2245 }
2246 if(methodIn && methodOut) {
2247 DoIntermediateSystematics(
2248 methodIn,
2249 methodOut,
2250 relativeErrorMethodInUp, // pointer reference
2251 relativeErrorMethodInLow, // pointer reference
2252 relativeErrorMethodOutUp, // pointer reference
2253 relativeErrorMethodOutLow, // pointer reference
2254 relativeStatisticalErrorIn,
2255 relativeStatisticalErrorOut,
2256 nominal,
2257 nominalIn,
2258 nominalOut,
2259 columns,
2260 rangeLow,
2261 rangeUp,
2262 readMe,
2263 "method"
2264 );
2265 if(relativeErrorMethodInUp) {
2266 TCanvas* relativeErrorMethod(new TCanvas("relativeErrorMethod", "relativeErrorMethod"));
2267 relativeErrorMethod->Divide(2);
2268 relativeErrorMethod->cd(1);
2269 Style(gPad, "GRID");
2270 relativeErrorMethodInUp->DrawCopy("b");
2271 relativeErrorMethodInLow->DrawCopy("same b");
2272 Style(AddLegend(gPad));
2273 relativeErrorMethod->cd(2);
2274 Style(gPad, "GRID");
2275 relativeErrorMethodOutUp->DrawCopy("b");
2276 relativeErrorMethodOutLow->DrawCopy("same b");
2277 SavePadToPDF(relativeErrorMethod);
2278 Style(AddLegend(gPad));
2279 relativeErrorMethod->Write();
2280 }
2281 }
2282
2283 // and the placeholder for the final systematic
2284 TH1D* relativeErrorInUp(new TH1D("max shape uncertainty in plane", "max shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2285 TH1D* relativeErrorOutUp(new TH1D("max shape uncertainty out of plane", "max shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2286 TH1D* relativeErrorInLow(new TH1D("min shape uncertainty in plane", "min shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2287 TH1D* relativeErrorOutLow(new TH1D("min shape uncertainty out of plane", "min shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2288 relativeErrorInUp->SetYTitle("relative uncertainty");
2289 relativeErrorOutUp->SetYTitle("relative uncertainty");
2290 relativeErrorInLow->SetYTitle("relative uncertainty");
2291 relativeErrorOutLow->SetYTitle("relative uncertainty");
2292
2293 // sum of squares for the total systematic error
2294 Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.), eInUp(0.);
2295 Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.), eOutUp(0.);
2296 Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.), eInLow(0.);
2297 Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.), eOutLow(0.);
2298
2299 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2300 // for the upper bound
2301 if(relativeErrorRegularizationInUp) aInUp = relativeErrorRegularizationInUp->GetBinContent(b+1);
2302 if(relativeErrorRegularizationOutUp) aOutUp = relativeErrorRegularizationOutUp->GetBinContent(b+1);
2303 if(relativeErrorTrueBinInUp) bInUp = relativeErrorTrueBinInUp->GetBinContent(b+1);
2304 if(relativeErrorTrueBinOutUp) bOutUp = relativeErrorTrueBinOutUp->GetBinContent(b+1);
2305 if(relativeErrorRecBinInUp) cInUp = relativeErrorRecBinInUp->GetBinContent(b+1);
2306 if(relativeErrorRecBinOutUp) cOutUp = relativeErrorRecBinOutUp->GetBinContent(b+1);
2307 if(relativeErrorMethodInUp) dInUp = relativeErrorMethodInUp->GetBinContent(b+1);
2308 if(relativeErrorMethodOutUp) dOutUp = relativeErrorMethodOutUp->GetBinContent(b+1);
2309 eInUp = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp + dInUp*dInUp;
2310 if(eInUp > 0) relativeErrorInUp->SetBinContent(b+1, 1.*TMath::Sqrt(eInUp));
2311 eOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp + dOutUp*dOutUp;
2312 if(eOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, 1.*TMath::Sqrt(eOutUp));
2313 // for the lower bound
2314 if(relativeErrorRegularizationInLow) aInLow = relativeErrorRegularizationInLow->GetBinContent(b+1);
2315 if(relativeErrorRegularizationOutLow) aOutLow = relativeErrorRegularizationOutLow->GetBinContent(b+1);
2316 if(relativeErrorTrueBinInLow) bInLow = relativeErrorTrueBinInLow->GetBinContent(b+1);
2317 if(relativeErrorTrueBinOutLow) bOutLow = relativeErrorTrueBinOutLow->GetBinContent(b+1);
2318 if(relativeErrorRecBinInLow) cInLow = relativeErrorRecBinInLow->GetBinContent(b+1);
2319 if(relativeErrorRecBinOutLow) cOutLow = relativeErrorRecBinOutLow->GetBinContent(b+1);
2320 if(relativeErrorMethodInLow) dInLow = relativeErrorMethodInLow->GetBinContent(b+1);
2321 if(relativeErrorMethodOutLow) dOutLow = relativeErrorMethodOutLow->GetBinContent(b+1);
2322 if(fSymmRMS) { // take first category as symmetric
2323 aInLow = aInUp;
2324 aOutLow = aOutUp;
2325 cInLow = cInUp;
2326 cOutLow = cOutUp; // other sources
2327 if(dInLow < dInUp) dInLow = dInUp;
2328 if(dOutLow < dOutUp) dOutLow = dOutUp;
2329 }
2330
2331 eInLow = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow + dInLow*dInLow;
2332 if(eInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(eInLow));
2333 eOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow + dOutLow*dOutLow;
2334 if(eOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(eOutLow));
2335 }
2336 // project the estimated errors on the nominal ratio
2337 if(nominal) {
2338 Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
2339 Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
2340 Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
2341 Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
2342 Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
2343 Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
2344 Double_t lowErr(0.), upErr(0.);
2345 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
2346 // add the in and out of plane errors in quadrature
2347 // take special care here: to propagate the assymetric error, we need to correlate the
2348 // InLow with OutUp (minimum value of ratio) and InUp with OutLow (maximum value of ratio)
2349 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(1+i)*relativeErrorOutUp->GetBinContent(i+1);
2350 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
2351 // set the errors
2352 ayl[i] = TMath::Sqrt(lowErr)*nominal->GetBinContent(i+1);
2353 ayh[i] = TMath::Sqrt(upErr)*nominal->GetBinContent(i+1);
2354 // get the bin width (which is the 'error' on x
2355 Double_t binWidth(nominal->GetBinWidth(i+1));
2356 axl[i] = binWidth/2.;
2357 axh[i] = binWidth/2.;
2358 // now get the coordinate for the point
2359 ax[i] = nominal->GetBinCenter(i+1);
2360 ay[i] = nominal->GetBinContent(i+1);
2361 }
2362 // save the nominal ratio
2363 TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
2364 shapeRatio = (TGraphAsymmErrors*)nominalError->Clone();
2365 nominalError->SetName("shape uncertainty");
2366 TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
2367 nominalCanvas->Divide(2);
2368 nominalCanvas->cd(1);
2369 Style(nominal, kBlack);
2370 Style(nominalError, kCyan, kRatio);
2371 nominalError->SetLineColor(kCyan);
2372 nominalError->SetMarkerColor(kCyan);
2373 nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2374 nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
2375 nominalError->DrawClone("a2");
2376 nominal->DrawCopy("same E1");
2377 Style(AddLegend(gPad));
2378 Style(gPad, "GRID");
2379 Style(nominalCanvas);
2380 // save nominal v2 and systematic errors
2381 TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
2382 nominalIn,
2383 nominalOut,
2384 fEventPlaneRes,
2385 "v_{2} with shape uncertainty",
2386 relativeErrorInUp,
2387 relativeErrorInLow,
2388 relativeErrorOutUp,
2389 relativeErrorOutLow,
2390 corr));
2391 shapeV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
2392 TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
2393 nominalCanvas->cd(2);
2394 Style(nominalV2, kBlack);
2395 Style(nominalV2Error, kCyan, kV2);
2396 nominalV2Error->SetLineColor(kCyan);
2397 nominalV2Error->SetMarkerColor(kCyan);
2398 nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2399 nominalV2Error->DrawClone("a2");
2400 nominalV2->DrawClone("same E1");
2401 Style(AddLegend(gPad));
2402 Style(gPad, "GRID");
2403 Style(nominalCanvas);
2404 SavePadToPDF(nominalCanvas);
2405 nominalCanvas->Write();
2406 }
2407
2408 TCanvas* relativeError(new TCanvas("relativeError"," relativeError"));
2409 relativeError->Divide(2);
2410 relativeError->cd(1);
2411 Style(gPad, "GRID");
2412 relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2413 Style(relativeErrorInUp, kBlue, kBar);
2414 Style(relativeErrorInLow, kGreen, kBar);
2415 relativeErrorInUp->DrawCopy("b");
2416 relativeErrorInLow->DrawCopy("same b");
2417 Style(relativeStatisticalErrorIn, kRed);
2418 relativeStatisticalErrorIn->DrawCopy("same");
2419 Style(AddLegend(gPad));
2420 relativeError->cd(2);
2421 Style(gPad, "GRID");
2422 relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2423 Style(relativeErrorOutUp, kBlue, kBar);
2424 Style(relativeErrorOutLow, kGreen, kBar);
2425 relativeErrorOutUp->DrawCopy("b");
2426 relativeErrorOutLow->DrawCopy("same b");
2427 Style(relativeStatisticalErrorOut, kRed);
2428 relativeStatisticalErrorOut->DrawCopy("same");
2429 Style(AddLegend(gPad));
2430
2431 // write the buffered file to disk and close the file
2432 SavePadToPDF(relativeError);
2433 relativeError->Write();
2434 output->Write();
2435// output->Close();
2436}
2437//_____________________________________________________________________________
2438 void AliJetFlowTools::DoIntermediateSystematics(
2439 TArrayI* variationsIn, // variantions in plane
2440 TArrayI* variationsOut, // variantions out of plane
2441 TH1D*& relativeErrorInUp, // pointer reference to minimum relative error histogram in plane
2442 TH1D*& relativeErrorInLow, // pointer reference to maximum relative error histogram in plane
2443 TH1D*& relativeErrorOutUp, // pointer reference to minimum relative error histogram out of plane
2444 TH1D*& relativeErrorOutLow, // pointer reference to maximum relative error histogram out of plane
2445 TH1D*& relativeStatisticalErrorIn, // relative systematic error on ratio
2446 TH1D*& relativeStatisticalErrorOut, // relative systematic error on ratio
2447 TH1D*& nominal, // clone of the nominal ratio in / out of plane
2448 TH1D*& nominalIn, // clone of the nominal in plane yield
2449 TH1D*& nominalOut, // clone of the nominal out of plane yield
2450 Int_t columns, // divide the output canvasses in this many columns
2451 Float_t rangeLow, // lower pt range
2452 Float_t rangeUp, // upper pt range
2453 TFile* readMe, // input file name (created by this unfolding class)
2454 TString source, // source of the variation
2455 Bool_t RMS // return RMS of distribution of variations as error
2456 ) const
2457{
2458 // intermediate systematic check function. first index of supplied array is nominal value
2459 TList* listOfKeys((TList*)readMe->GetListOfKeys());
2460 if(!listOfKeys) {
2461 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2462 return;
2463 }
2464 // check input params
2465 if(variationsIn->GetSize() != variationsOut->GetSize()) {
2466 printf(" > DoSystematics: fatal error, input arrays have different sizes ! < \n ");
2467 return;
2468 }
2469 TDirectoryFile* defRootDirIn(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(0))->GetName())));
2470 TDirectoryFile* defRootDirOut(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(0))->GetName())));
2471 if(!(defRootDirIn && defRootDirOut)) {
2472 printf(" > DoSystematics: fatal error, couldn't retrieve nominal values ! < \n ");
2473 return;
2474 }
2475 TString defIn(defRootDirIn->GetName());
2476 TString defOut(defRootDirOut->GetName());
2477
2478 // define lines to make the output prettier
2479 TLine* lineLow(new TLine(rangeLow, 0., rangeLow, 2.));
2480 TLine* lineUp(new TLine(rangeUp, 0., rangeUp, 2.));
2481 lineLow->SetLineColor(11);
2482 lineUp->SetLineColor(11);
2483 lineLow->SetLineWidth(3);
2484 lineUp->SetLineWidth(3);
2485
2486 // define an output histogram with the maximum relative error from this systematic constribution
2487 // if the option RMS is set to false, sigma is not really a standard deviation but holds the maximum (or minimum) relative value that the data has
2488 // reached in this function call.
2489 // if the option RMS is set to true, sigma holds the RMS value (equal to sigma as the mean is zero for relative errors) of the distribution of variations
2490 // which should correspond to a 68% confidence level
2491 relativeErrorInUp = new TH1D(Form("relative error (up) from %s", source.Data()), Form("relative error (up) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2492 relativeErrorInLow = new TH1D(Form("relative error (low) from %s", source.Data()), Form("relative error (low) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2493 relativeErrorOutUp = new TH1D(Form("relative error (up) from %s", source.Data()), Form("relative error (up) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2494 relativeErrorOutLow = new TH1D(Form("relative error (low) from %s", source.Data()), Form("relative error (low) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2495 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2496 if(!RMS) {
2497 relativeErrorInUp->SetBinContent(b+1, 1.);
2498 relativeErrorInUp->SetBinError(b+1, 0.);
2499 relativeErrorOutUp->SetBinContent(b+1, 1.);
2500 relativeErrorOutUp->SetBinError(b+1, .0);
2501 relativeErrorInLow->SetBinContent(b+1, 1.);
2502 relativeErrorInLow->SetBinError(b+1, 0.);
2503 relativeErrorOutLow->SetBinContent(b+1, 1.);
2504 relativeErrorOutLow->SetBinError(b+1, .0);
2505 } else if(RMS) {
2506 relativeErrorInUp->SetBinContent(b+1, 0.);
2507 relativeErrorInUp->SetBinError(b+1, 0.);
2508 relativeErrorOutUp->SetBinContent(b+1, 0.);
2509 relativeErrorOutUp->SetBinError(b+1, 0.);
2510 relativeErrorInLow->SetBinContent(b+1, 0.);
2511 relativeErrorInLow->SetBinError(b+1, 0.);
2512 relativeErrorOutLow->SetBinContent(b+1, 0.);
2513 relativeErrorOutLow->SetBinError(b+1, 0.);
2514 }
2515 }
2516 Int_t relativeErrorInUpN[100] = {0};
2517 Int_t relativeErrorOutUpN[100] = {0};
2518 Int_t relativeErrorInLowN[100] = {0};
2519 Int_t relativeErrorOutLowN[100] = {0};
2520
2521 // define an output histogram with the systematic error from this systematic constribution
2522 if(!relativeStatisticalErrorIn && !relativeStatisticalErrorOut) {
2523 relativeStatisticalErrorIn = new TH1D("relative statistical error, in plane", "relative statistical error, in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2524 relativeStatisticalErrorOut = new TH1D("relative statistical error, out of plane", "relative statistital error, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2525 }
2526
2527 // prepare necessary canvasses
2528 TCanvas* canvasIn(new TCanvas(Form("SYST_%s_PearsonIn", source.Data()), Form("SYST_%s_PearsonIn", source.Data())));
2529 TCanvas* canvasOut(0x0);
2530 if(fDphiUnfolding) canvasOut = new TCanvas(Form("SYST_%s_PearsonOut", source.Data()), Form("SYST_%s_PearsonOut", source.Data()));
2531 TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas(Form("SYST_%s_RefoldedIn", source.Data()), Form("SYST_%s_RefoldedIn", source.Data())));
2532 TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
2533 if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas(Form("SYST_%s_RefoldedOut", source.Data()), Form("SYST_%s_RefoldedOut", source.Data()));
2534 TCanvas* canvasSpectraIn(new TCanvas(Form("SYST_%s_SpectraIn", source.Data()), Form("SYST_%s_SpectraIn", source.Data())));
2535 TCanvas* canvasSpectraOut(0x0);
2536 if(fDphiUnfolding) canvasSpectraOut = new TCanvas(Form("SYST_%s_SpectraOut", source.Data()), Form("SYST_%s_SpectraOut", source.Data()));
2537 TCanvas* canvasRatio(0x0);
2538 if(fDphiUnfolding) canvasRatio = new TCanvas(Form("SYST_%s_Ratio", source.Data()), Form("SYST_%s_Ratio", source.Data()));
2539 TCanvas* canvasV2(0x0);
2540 if(fDphiUnfolding) canvasV2 = new TCanvas(Form("SYST_%s_V2", source.Data()), Form("SYST_%s_V2", source.Data()));
2541 TCanvas* canvasMISC(new TCanvas(Form("SYST_%s_MISC", source.Data()), Form("SYST_%s_MISC", source.Data())));
2542 TCanvas* canvasMasterIn(new TCanvas(Form("SYST_%s_defaultIn", source.Data()), Form("SYST_%s_defaultIn", source.Data())));
2543 TCanvas* canvasMasterOut(0x0);
2544 if(fDphiUnfolding) canvasMasterOut = new TCanvas(Form("SYST_%s_defaultOut", source.Data()), Form("SYST_%s_defaultOut", source.Data()));
2545 (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
2546
2547 TCanvas* canvasProfiles(new TCanvas(Form("SYST_%s_canvasProfiles", source.Data()), Form("SYST_%s_canvasProfiles", source.Data())));
2548 canvasProfiles->Divide(2);
2549 TProfile* ratioProfile(new TProfile(Form("SYST_%s_ratioProfile", source.Data()), Form("SYST_%s_ratioProfile", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2550 TProfile* v2Profile(new TProfile(Form("SYST_%s_v2Profile", source.Data()), Form("SYST_%s_v2Profile", source.Data()),fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2551 // get an estimate of the number of outputs and find the default set
2552
2553 Int_t rows = 1;
2554 columns = variationsIn->GetSize()-1;
2555 (TMath::Floor(variationsIn->GetSize()/(float)columns)+((variationsIn->GetSize()%columns)>0));
2556 canvasIn->Divide(columns, rows);
2557 if(canvasOut) canvasOut->Divide(columns, rows);
2558 canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
2559 if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
2560 canvasSpectraIn->Divide(columns, rows);
2561 if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
2562 if(canvasRatio) canvasRatio->Divide(columns, rows);
2563 if(canvasV2) canvasV2->Divide(columns, rows);
2564 canvasMasterIn->Divide(columns, rows);
2565 if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
2566
2567 // prepare a separate set of canvases to hold the nominal points
2568 TCanvas* canvasNominalIn(new TCanvas(Form("Nominal_%s_PearsonIn", source.Data()), Form("Nominal_%s_PearsonIn", source.Data())));
2569 TCanvas* canvasNominalOut(0x0);
2570 if(fDphiUnfolding) canvasNominalOut = new TCanvas(Form("Nominal_%s_PearsonOut", source.Data()), Form("Nominal_%s_PearsonOut", source.Data()));
2571 TCanvas* canvasNominalRatioMeasuredRefoldedIn(new TCanvas(Form("Nominal_%s_RefoldedIn", source.Data()), Form("Nominal_%s_RefoldedIn", source.Data())));
2572 TCanvas* canvasNominalRatioMeasuredRefoldedOut(0x0);
2573 if(fDphiUnfolding) canvasNominalRatioMeasuredRefoldedOut = new TCanvas(Form("Nominal_%s_RefoldedOut", source.Data()), Form("Nominal_%s_RefoldedOut", source.Data()));
2574 TCanvas* canvasNominalSpectraIn(new TCanvas(Form("Nominal_%s_SpectraIn", source.Data()), Form("Nominal_%s_SpectraIn", source.Data())));
2575 TCanvas* canvasNominalSpectraOut(0x0);
2576 if(fDphiUnfolding) canvasNominalSpectraOut = new TCanvas(Form("Nominal_%s_SpectraOut", source.Data()), Form("Nominal_%s_SpectraOut", source.Data()));
2577 TCanvas* canvasNominalRatio(0x0);
2578 if(fDphiUnfolding) canvasNominalRatio = new TCanvas(Form("Nominal_%s_Ratio", source.Data()), Form("Nominal_%s_Ratio", source.Data()));
2579 TCanvas* canvasNominalV2(0x0);
2580 if(fDphiUnfolding) canvasNominalV2 = new TCanvas(Form("Nominal_%s_V2", source.Data()), Form("Nominal_%s_V2", source.Data()));
2581 TCanvas* canvasNominalMISC(new TCanvas(Form("Nominal_%s_MISC", source.Data()), Form("Nominal_%s_MISC", source.Data())));
2582 TCanvas* canvasNominalMasterIn(new TCanvas(Form("Nominal_%s_defaultIn", source.Data()), Form("Nominal_%s_defaultIn", source.Data())));
2583 TCanvas* canvasNominalMasterOut(0x0);
2584 if(fDphiUnfolding) canvasNominalMasterOut = new TCanvas(Form("Nominal_%s_defaultOut", source.Data()), Form("Nominal_%s_defaultOut", source.Data()));
2585 (fDphiUnfolding) ? canvasNominalMISC->Divide(4, 2) : canvasNominalMISC->Divide(4, 1);
2586
2587 canvasNominalSpectraIn->Divide(2);
2588 if(canvasNominalSpectraOut) canvasNominalSpectraOut->Divide(2);
2589
2590 canvasNominalMasterIn->Divide(2);
2591 if(canvasNominalMasterOut) canvasNominalMasterOut->Divide(2);
2592
2593 // extract the default output
2594 TH1D* defaultUnfoldedJetSpectrumIn(0x0);
2595 TH1D* defaultUnfoldedJetSpectrumOut(0x0);
2596 TDirectoryFile* defDirIn = (TDirectoryFile*)defRootDirIn->Get(Form("InPlane___%s", defIn.Data()));
2597 TDirectoryFile* defDirOut = (TDirectoryFile*)defRootDirOut->Get(Form("OutOfPlane___%s", defOut.Data()));
2598 if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", defIn.Data()));
2599 if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", defOut.Data()));
2600 printf(" > succesfully extracted default results < \n");
2601
2602 // loop through the directories, only plot the graphs if the deconvolution converged
2603 TDirectoryFile* tempDirIn(0x0);
2604 TDirectoryFile* tempDirOut(0x0);
2605 TDirectoryFile* tempIn(0x0);
2606 TDirectoryFile* tempOut(0x0);
2607 TH1D* unfoldedSpectrumInForRatio(0x0);
2608 TH1D* unfoldedSpectrumOutForRatio(0x0);
2609 for(Int_t i(0), j(-1); i < variationsIn->GetSize(); i++) {
2610 tempDirIn = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(i))->GetName())));
2611 tempDirOut = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(i))->GetName())));
2612 if(!(tempDirIn && tempDirOut)) {
2613 printf(" > DoSystematics: couldn't get a set of variations < \n");
2614 continue;
2615 }
2616 TString dirNameIn(tempDirIn->GetName());
2617 TString dirNameOut(tempDirOut->GetName());
2618 // try to read the in- and out of plane subdirs
2619 tempIn = (TDirectoryFile*)tempDirIn->Get(Form("InPlane___%s", dirNameIn.Data()));
2620 tempOut = (TDirectoryFile*)tempDirOut->Get(Form("OutOfPlane___%s", dirNameOut.Data()));
2621 j++;
2622 if(tempIn) {
2623 // to see if the unfolding converged try to extract the pearson coefficients
2624 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirNameIn.Data())));
2625 if(pIn) {
2626 printf(" - %s in plane converged \n", dirNameIn.Data());
2627 canvasIn->cd(j);
2628 if(i==0) canvasNominalIn->cd(j);
2629 Style(gPad, "PEARSON");
2630 pIn->DrawCopy("colz");
2631 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirNameIn.Data())));
2632 if(rIn) {
2633 Style(rIn);
2634 printf(" > found RatioRefoldedMeasured < \n");
2635 canvasRatioMeasuredRefoldedIn->cd(j);
2636 if(i==0) canvasNominalRatioMeasuredRefoldedIn->cd(j);
2637 Style(gPad, "GRID");
2638 rIn->SetFillColor(kRed);
2639 rIn->Draw("ap");
2640 }
2641 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2642 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2643 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirNameIn.Data())));
2644 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirNameIn.Data())));
2645 if(dvector && avalue && rm && eff) {
2646 Style(dvector);
2647 Style(avalue);
2648 Style(rm);
2649 Style(eff);
2650 canvasMISC->cd(1);
2651 if(i==0) canvasNominalMISC->cd(1);
2652 Style(gPad, "SPECTRUM");
2653 dvector->DrawCopy();
2654 canvasMISC->cd(2);
2655 if(i==0) canvasNominalMISC->cd(2);
2656 Style(gPad, "SPECTRUM");
2657 avalue->DrawCopy();
2658 canvasMISC->cd(3);
2659 if(i==0) canvasNominalMISC->cd(3);
2660 Style(gPad, "PEARSON");
2661 rm->DrawCopy("colz");
2662 canvasMISC->cd(4);
2663 if(i==0) canvasNominalMISC->cd(4);
2664 Style(gPad, "GRID");
2665 eff->DrawCopy();
2666 } else if(rm && eff) {
2667 Style(rm);
2668 Style(eff);
2669 canvasMISC->cd(3);
2670 if(i==0) canvasNominalMISC->cd(3);
2671 Style(gPad, "PEARSON");
2672 rm->DrawCopy("colz");
2673 canvasMISC->cd(4);
2674 if(i==0) canvasNominalMISC->cd(4);
2675 Style(gPad, "GRID");
2676 eff->DrawCopy();
2677 }
2678 }
2679 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirNameIn.Data())));
2680 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirNameIn.Data())));
2681 unfoldedSpectrumInForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2682 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirNameIn.Data())));
2683 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2684 if(defaultUnfoldedJetSpectrumIn) {
2685 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2686 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirNameIn.Data())));
2687 temp->Divide(unfoldedSpectrum);
2688 // get the absolute relative error
2689 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2690 if(!RMS) { // save the maximum deviation that a variation can cause
2691 // the variation is HIGHER than the nominal point, so the bar goes UP
2692 if( temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorInUp->GetBinContent(b+1)) {
2693 relativeErrorInUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2694 relativeErrorInUp->SetBinError(b+1, 0.);
2695 }
2696 // the variation is LOWER than the nominal point, so the bar goes DOWN
2697 else if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorInLow->GetBinContent(b+1)) {
2698 relativeErrorInLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2699 relativeErrorInLow->SetBinError(b+1, 0.);
2700 }
2701 } else if (RMS && !fSymmRMS) { // save info necessary for evaluating the RMS of a distribution of variations
2702 printf(" oops shouldnt be here \n " );
2703 if(temp->GetBinContent(b+1) < 1) {
2704 relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2705 relativeErrorInUpN[b]++;
2706 }
2707 // the variation is LOWER than the nominal point, so the bar goes DOWN
2708 else if(temp->GetBinContent(b+1) > 1) {
2709 relativeErrorInLow->SetBinContent(b+1, relativeErrorInLow->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2710 relativeErrorInLowN[b]++;
2711 }
2712 } else if (fSymmRMS) {
2713 // save symmetric sum of square to get a symmetric rms
2714 relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(temp->GetBinContent(b+1)-1., 2));
2715 relativeErrorInUpN[b]++;
2716 }
2717 if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorIn->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2718 }
2719 temp->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2720 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2721 temp->GetYaxis()->SetTitle("ratio");
2722 canvasMasterIn->cd(j);
2723 temp->GetYaxis()->SetRangeUser(0., 2);
2724 Style(gPad, "GRID");
2725 temp->DrawCopy();
2726 canvasNominalMasterIn->cd(1);
2727 Style(gPad, "GRID");
2728 if(i > 0 ) {
2729 TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2730 tempSyst->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2731 Style(tempSyst, (EColor)(i+2));
2732 if(i==1) tempSyst->DrawCopy();
2733 else tempSyst->DrawCopy("same");
2734 }
2735 }
2736 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirNameIn.Data())));
2737 canvasSpectraIn->cd(j);
2738 if(i==0) canvasNominalSpectraIn->cd(1);
2739 Style(gPad);
2740 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2741 unfoldedSpectrum->DrawCopy();
2742 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2743 inputSpectrum->DrawCopy("same");
2744 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2745 refoldedSpectrum->DrawCopy("same");
2746 TLegend* l(AddLegend(gPad));
2747 Style(l);
2748 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2749 Float_t chi(fitStatus->GetBinContent(1));
2750 Float_t pen(fitStatus->GetBinContent(2));
2751 Int_t dof((int)fitStatus->GetBinContent(3));
2752 Float_t beta(fitStatus->GetBinContent(4));
2753 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2754 } else if (fitStatus) { // only available in SVD fit
2755 Int_t reg((int)fitStatus->GetBinContent(1));
2756 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2757 }
2758 canvasNominalSpectraIn->cd(2);
2759 TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2760 tempSyst->SetTitle(Form("[%s]", dirNameIn.Data()));
2761 Style(tempSyst, (EColor)(i+2));
2762 Style(gPad, "SPECTRUM");
2763 if(i==0) tempSyst->DrawCopy();
2764 else tempSyst->DrawCopy("same");
2765 }
2766 }
2767 if(tempOut) {
2768 TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirNameOut.Data())));
2769 if(pOut) {
2770 printf(" - %s out of plane converged \n", dirNameOut.Data());
2771 canvasOut->cd(j);
2772 if(i==0) canvasNominalOut->cd(j);
2773 Style(gPad, "PEARSON");
2774 pOut->DrawCopy("colz");
2775 TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirNameOut.Data())));
2776 if(rOut) {
2777 Style(rOut);
2778 printf(" > found RatioRefoldedMeasured < \n");
2779 canvasRatioMeasuredRefoldedOut->cd(j);
2780 if(i==0) canvasNominalRatioMeasuredRefoldedOut->cd(j);
2781 Style(gPad, "GRID");
2782 rOut->SetFillColor(kRed);
2783 rOut->Draw("ap");
2784 }
2785 TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2786 TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2787 TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirNameOut.Data())));
2788 TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirNameOut.Data())));
2789 if(dvector && avalue && rm && eff) {
2790 Style(dvector);
2791 Style(avalue);
2792 Style(rm);
2793 Style(eff);
2794 canvasMISC->cd(5);
2795 if(i==0) canvasNominalMISC->cd(5);
2796 Style(gPad, "SPECTRUM");
2797 dvector->DrawCopy();
2798 canvasMISC->cd(6);
2799 if(i==0) canvasNominalMISC->cd(6);
2800 Style(gPad, "SPECTRUM");
2801 avalue->DrawCopy();
2802 canvasMISC->cd(7);
2803 if(i==0) canvasNominalMISC->cd(7);
2804 Style(gPad, "PEARSON");
2805 rm->DrawCopy("colz");
2806 canvasMISC->cd(8);
2807 if(i==0) canvasNominalMISC->cd(8);
2808 Style(gPad, "GRID");
2809 eff->DrawCopy();
2810 } else if(rm && eff) {
2811 Style(rm);
2812 Style(eff);
2813 canvasMISC->cd(7);
2814 if(i==0) canvasNominalMISC->cd(7);
2815 Style(gPad, "PEARSON");
2816 rm->DrawCopy("colz");
2817 canvasMISC->cd(8);
2818 if(i==0) canvasNominalMISC->cd(8);
2819 Style(gPad, "GRID");
2820 eff->DrawCopy();
2821 }
2822 }
2823 TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirNameOut.Data())));
2824 TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirNameOut.Data())));
2825 unfoldedSpectrumOutForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2826 TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirNameOut.Data())));
2827 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2828 if(defaultUnfoldedJetSpectrumOut) {
2829 Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
2830 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirNameOut.Data())));
2831 temp->Divide(unfoldedSpectrum);
2832 // get the absolute relative error
2833 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2834 if(!RMS) {
2835 // check if the error is larger than the current maximum
2836 if(temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorOutUp->GetBinContent(b+1)) {
2837 relativeErrorOutUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2838 relativeErrorOutUp->SetBinError(b+1, 0.);
2839 }
2840 // check if the error is smaller than the current minimum
2841 else if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorOutLow->GetBinContent(b+1)) {
2842 relativeErrorOutLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2843 relativeErrorOutLow->SetBinError(b+1, 0.);
2844 }
2845 } else if (RMS && !fSymmRMS) {
2846 printf(" OOps \n ");
2847 if(temp->GetBinContent(b+1) < 1) {
2848 relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2849 relativeErrorOutUpN[b]++;
2850 }
2851 else if(temp->GetBinContent(b+1) > 1) {
2852 relativeErrorOutLow->SetBinContent(b+1, relativeErrorOutLow->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2853 relativeErrorOutLowN[b]++;
2854 }
2855 } else if (fSymmRMS) {
2856 // save symmetric rms value
2857 relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(temp->GetBinContent(b+1)-1., 2));
2858 relativeErrorOutUpN[b]++;
2859 }
2860 if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorOut->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2861 }
2862 temp->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2863 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2864 temp->GetYaxis()->SetTitle("ratio");
2865 canvasMasterOut->cd(j);
2866 temp->GetYaxis()->SetRangeUser(0., 2);
2867 Style(gPad, "GRID");
2868 temp->DrawCopy();
2869 canvasNominalMasterOut->cd(1);
2870 Style(gPad, "GRID");
2871 if(i > 0 ) {
2872 TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2873 tempSyst->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2874 Style(tempSyst, (EColor)(i+2));
2875 if(i==1) tempSyst->DrawCopy();
2876 else tempSyst->DrawCopy("same");
2877 }
2878 }
2879 TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirNameOut.Data())));
2880 canvasSpectraOut->cd(j);
2881 if(i==0) canvasNominalSpectraOut->cd(1);
2882 Style(gPad);
2883 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2884 unfoldedSpectrum->DrawCopy();
2885 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2886 inputSpectrum->DrawCopy("same");
2887 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2888 refoldedSpectrum->DrawCopy("same");
2889 TLegend* l(AddLegend(gPad));
2890 Style(l);
2891 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2892 Float_t chi(fitStatus->GetBinContent(1));
2893 Float_t pen(fitStatus->GetBinContent(2));
2894 Int_t dof((int)fitStatus->GetBinContent(3));
2895 Float_t beta(fitStatus->GetBinContent(4));
2896 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2897 } else if (fitStatus) { // only available in SVD fit
2898 Int_t reg((int)fitStatus->GetBinContent(1));
2899 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2900 }
2901 canvasNominalSpectraOut->cd(2);
2902 TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2903 tempSyst->SetTitle(Form("[%s]", dirNameOut.Data()));
2904 Style(tempSyst, (EColor)(i+2));
2905 Style(gPad, "SPECTRUM");
2906 if(i==0) tempSyst->DrawCopy();
2907 else tempSyst->DrawCopy("same");
2908 }
2909 }
2910 if(canvasRatio && canvasV2) {
2911 canvasRatio->cd(j);
2912 if(i==0) {
2913 canvasNominalRatio->cd(j);
2914 if(nominal && nominalIn && nominalOut) {
2915 // if a nominal ratio is requested, delete the placeholder and update the nominal point
2916 delete nominal;
2917 delete nominalIn;
2918 delete nominalOut;
2919 nominalIn = (TH1D*)unfoldedSpectrumInForRatio->Clone("in plane jet yield");
2920 nominalOut = (TH1D*)unfoldedSpectrumOutForRatio->Clone("out of plane jet yield");
2921 nominal = (TH1D*)unfoldedSpectrumInForRatio->Clone("ratio in plane out of plane");
2922 nominal->Divide(nominal, unfoldedSpectrumOutForRatio); // standard root divide for uncorrelated histos
2923 }
2924 }
2925 TGraphErrors* ratioYield(GetRatio(unfoldedSpectrumInForRatio, unfoldedSpectrumOutForRatio, TString(Form("ratio [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2926 Double_t _tempx(0), _tempy(0);
2927 if(ratioYield) {
2928 Style(ratioYield);
2929 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2930 ratioYield->GetPoint(b, _tempx, _tempy);
2931 ratioProfile->Fill(_tempx, _tempy);
2932 }
2933 ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
2934 ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2935 ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
2936 ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2937 ratioYield->SetFillColor(kRed);
2938 ratioYield->Draw("ap");
2939 }
2940 canvasV2->cd(j);
2941 if(i==0) canvasNominalV2->cd(j);
2942 TGraphErrors* ratioV2(GetV2(unfoldedSpectrumInForRatio,unfoldedSpectrumOutForRatio, fEventPlaneRes, TString(Form("v_{2} [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2943 if(ratioV2) {
2944 Style(ratioV2);
2945 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2946 ratioV2->GetPoint(b, _tempx, _tempy);
2947 v2Profile->Fill(_tempx, _tempy);
2948
2949 }
2950 v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
2951 v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2952 ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
2953 ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2954 ratioV2->SetFillColor(kRed);
2955 ratioV2->Draw("ap");
2956 }
2957 }
2958 delete unfoldedSpectrumInForRatio;
2959 delete unfoldedSpectrumOutForRatio;
2960 }
2961 // save the canvasses
2962 canvasProfiles->cd(1);
2963 Style(ratioProfile);
2964 ratioProfile->DrawCopy();
2965 canvasProfiles->cd(2);
2966 Style(v2Profile);
2967 v2Profile->DrawCopy();
2968 SavePadToPDF(canvasProfiles);
2969 canvasProfiles->Write();
2970 SavePadToPDF(canvasIn);
2971 canvasIn->Write();
2972 if(canvasOut) {
2973 SavePadToPDF(canvasOut);
2974 canvasOut->Write();
2975 }
2976 SavePadToPDF(canvasRatioMeasuredRefoldedIn);
2977 canvasRatioMeasuredRefoldedIn->Write();
2978 if(canvasRatioMeasuredRefoldedOut) {
2979 SavePadToPDF(canvasRatioMeasuredRefoldedOut);
2980 canvasRatioMeasuredRefoldedOut->Write();
2981 }
2982 SavePadToPDF(canvasSpectraIn);
2983 canvasSpectraIn->Write();
2984 if(canvasSpectraOut) {
2985 SavePadToPDF(canvasSpectraOut);
2986 canvasSpectraOut->Write();
2987 SavePadToPDF(canvasRatio);
2988 canvasRatio->Write();
2989 SavePadToPDF(canvasV2);
2990 canvasV2->Write();
2991 }
2992 SavePadToPDF(canvasMasterIn);
2993 canvasMasterIn->Write();
2994 if(canvasMasterOut) {
2995 SavePadToPDF(canvasMasterOut);
2996 canvasMasterOut->Write();
2997 }
2998 SavePadToPDF(canvasMISC);
2999 canvasMISC->Write();
3000 // save the nomial canvasses
3001 SavePadToPDF(canvasNominalIn);
3002 canvasNominalIn->Write();
3003 if(canvasNominalOut) {
3004 SavePadToPDF(canvasNominalOut);
3005 canvasNominalOut->Write();
3006 }
3007 SavePadToPDF(canvasNominalRatioMeasuredRefoldedIn);
3008 canvasNominalRatioMeasuredRefoldedIn->Write();
3009 if(canvasNominalRatioMeasuredRefoldedOut) {
3010 SavePadToPDF(canvasNominalRatioMeasuredRefoldedOut);
3011 canvasNominalRatioMeasuredRefoldedOut->Write();
3012 }
3013 canvasNominalSpectraIn->cd(2);
3014 Style(AddLegend(gPad));
3015 SavePadToPDF(canvasNominalSpectraIn);
3016 canvasNominalSpectraIn->Write();
3017 if(canvasNominalSpectraOut) {
3018 canvasNominalSpectraOut->cd(2);
3019 Style(AddLegend(gPad));
3020 SavePadToPDF(canvasNominalSpectraOut);
3021 canvasNominalSpectraOut->Write();
3022 SavePadToPDF(canvasNominalRatio);
3023 canvasNominalRatio->Write();
3024 SavePadToPDF(canvasNominalV2);
3025 canvasNominalV2->Write();
3026 }
3027 canvasNominalMasterIn->cd(1);
3028 Style(AddLegend(gPad));
3029 lineUp->DrawClone("same");
3030 lineLow->DrawClone("same");
3031 SavePadToPDF(canvasNominalMasterIn);
3032 canvasNominalMasterIn->Write();
3033 if(canvasNominalMasterOut) {
3034 canvasNominalMasterOut->cd(1);
3035 Style(AddLegend(gPad));
3036 lineUp->DrawClone("same");
3037 lineLow->DrawClone("same");
3038 SavePadToPDF(canvasNominalMasterOut);
3039 canvasNominalMasterOut->Write();
3040 }
3041 SavePadToPDF(canvasNominalMISC);
3042 canvasNominalMISC->Write();
3043
3044 // save the relative errors
3045 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
3046 // to arrive at a min and max from here, combine in up and out low
3047 if(!RMS) {
3048 relativeErrorInUp->SetBinContent(b+1, -1.*(relativeErrorInUp->GetBinContent(b+1)-1));
3049 relativeErrorInUp->SetBinError(b+1, 0.);
3050 relativeErrorOutUp->SetBinContent(b+1, -1.*(relativeErrorOutUp->GetBinContent(b+1)-1));
3051 relativeErrorOutUp->SetBinError(b+1, .0);
3052 relativeErrorInLow->SetBinContent(b+1, -1.*(relativeErrorInLow->GetBinContent(b+1)-1));
3053 relativeErrorInLow->SetBinError(b+1, 0.);
3054 relativeErrorOutLow->SetBinContent(b+1, -1.*(relativeErrorOutLow->GetBinContent(b+1)-1));
3055 relativeErrorOutLow->SetBinError(b+1, .0);
3056 } else if (RMS) {
3057 // these guys are already stored as percentages, so no need to remove the offset of 1
3058 // RMS is defined as sqrt(sum(squared))/N
3059 // min is defined as negative, max is defined as positive
3060 if(!fSymmRMS) {
3061 if(relativeErrorInUpN[b] < 1) relativeErrorInUpN[b] = 1;
3062 if(relativeErrorInLowN[b] < 1) relativeErrorInLowN[b] = 1;
3063 if(relativeErrorOutUpN[b] < 1) relativeErrorOutUpN[b] = 1;
3064 if(relativeErrorOutLowN[b] < 1) relativeErrorOutLowN[b] = 1;
3065 relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(b+1)/relativeErrorInUpN[b]));
3066 relativeErrorInUp->SetBinError(b+1, 0.);
3067 relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorOutUp->GetBinContent(b+1)/relativeErrorOutUpN[b]));
3068 relativeErrorOutUp->SetBinError(b+1, .0);
3069 relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(relativeErrorInLow->GetBinContent(b+1)/relativeErrorInLowN[b]));
3070 relativeErrorInLow->SetBinError(b+1, 0.);
3071 relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(relativeErrorOutLow->GetBinContent(b+1)/relativeErrorOutLowN[b]));
3072 relativeErrorOutLow->SetBinError(b+1, .0);
3073 } else if (fSymmRMS) {
3074 if(relativeErrorInUpN[b] < 1) relativeErrorInUpN[b] = 1;
3075 if(relativeErrorOutUpN[b] < 1) relativeErrorOutUpN[b] = 1;
3076 relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(b+1)/relativeErrorInUpN[b]));
3077 relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorOutUp->GetBinContent(b+1)/relativeErrorOutUpN[b]));
3078 }
3079 }
3080 }
3081 relativeErrorInUp->SetYTitle("relative uncertainty");
3082 relativeErrorOutUp->SetYTitle("relative uncertainty");
3083 relativeErrorInLow->SetYTitle("relative uncertainty");
3084 relativeErrorOutLow->SetYTitle("relative uncertainty");
3085 relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
3086 relativeErrorInLow->GetYaxis()->SetRangeUser(-1.5, 3.);
3087 relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
3088 relativeErrorOutLow->GetYaxis()->SetRangeUser(-1.5, 3.);
3089
3090 canvasNominalMasterIn->cd(2);
3091 Style(gPad, "GRID");
3092 Style(relativeErrorInUp, kBlue, kBar);
3093 Style(relativeErrorInLow, kGreen, kBar);
3094 relativeErrorInUp->DrawCopy("b");
3095 relativeErrorInLow->DrawCopy("same b");
3096 Style(AddLegend(gPad));
3097 SavePadToPDF(canvasNominalMasterIn);
3098 canvasNominalMasterIn->Write();
3099 canvasNominalMasterOut->cd(2);
3100 Style(gPad, "GRID");
3101 Style(relativeErrorOutUp, kBlue, kBar);
3102 Style(relativeErrorOutLow, kGreen, kBar);
3103 relativeErrorOutUp->DrawCopy("b");
3104 relativeErrorOutLow->DrawCopy("same b");
3105 Style(AddLegend(gPad));
3106 SavePadToPDF(canvasNominalMasterOut);
3107 canvasNominalMasterOut->Write();
3108}
3109//_____________________________________________________________________________
3110void AliJetFlowTools::PostProcess(TString def, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
3111{
3112 // go through the output file and perform post processing routines
3113 // can either be performed in one go with the unfolding, or at a later stage
3114 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
3115 TFile readMe(in.Data(), "READ"); // open file read-only
3116 if(readMe.IsZombie()) {
3117 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
3118 return;
3119 }
3120 printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
3121 readMe.ls();
3122 TList* listOfKeys((TList*)readMe.GetListOfKeys());
3123 if(!listOfKeys) {
3124 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
3125 return;
3126 }
3127 // prepare necessary canvasses
3128 TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
3129 TCanvas* canvasOut(0x0);
3130 if(fDphiUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
3131 TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
3132 TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
3133 if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
3134 TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
3135 TCanvas* canvasSpectraOut(0x0);
3136 if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
3137 TCanvas* canvasRatio(0x0);
3138 if(fDphiUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
3139 TCanvas* canvasV2(0x0);
3140 if(fDphiUnfolding) canvasV2 = new TCanvas("V2", "V2");
3141 TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
3142 TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
3143 TCanvas* canvasMasterOut(0x0);
3144 if(fDphiUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
3145 (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
3146 TDirectoryFile* defDir(0x0);
3147 TCanvas* canvasProfiles(new TCanvas("canvasProfiles", "canvasProfiles"));
3148 canvasProfiles->Divide(2);
3149 TProfile* ratioProfile(new TProfile("ratioProfile", "ratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3150 TProfile* v2Profile(new TProfile("v2Profile", "v2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3151 // get an estimate of the number of outputs and find the default set
3152 Int_t cacheMe(0);
3153 for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
3154 if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
3155 if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3156 cacheMe++;
3157 }
3158 }
3159 Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%columns)>0));
3160 canvasIn->Divide(columns, rows);
3161 if(canvasOut) canvasOut->Divide(columns, rows);
3162 canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
3163 if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
3164 canvasSpectraIn->Divide(columns, rows);
3165 if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
3166 if(canvasRatio) canvasRatio->Divide(columns, rows);
3167 if(canvasV2) canvasV2->Divide(columns, rows);
3168
3169 canvasMasterIn->Divide(columns, rows);
3170 if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
3171 // extract the default output
3172 TH1D* defaultUnfoldedJetSpectrumIn(0x0);
3173 TH1D* defaultUnfoldedJetSpectrumOut(0x0);
3174 if(defDir) {
3175 TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
3176 TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
3177 if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
3178 if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
3179 printf(" > succesfully extracted default results < \n");
3180 }
3181
3182 // loop through the directories, only plot the graphs if the deconvolution converged
3183 TDirectoryFile* tempDir(0x0);
3184 TDirectoryFile* tempIn(0x0);
3185 TDirectoryFile* tempOut(0x0);
3186 for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
3187 // read the directory by its name
3188 tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3189 if(!tempDir) continue;
3190 TString dirName(tempDir->GetName());
3191 // try to read the in- and out of plane subdirs
3192 tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
3193 tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
3194 j++;
3195 if(!tempIn) { // attempt to get MakeAU output
3196 TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
3197 TCanvas* tempPearson(new TCanvas(Form("pearson_%s", dirName.Data()), Form("pearson_%s", dirName.Data())));
3198 tempPearson->Divide(4,2);
3199 TCanvas* tempRatio(new TCanvas(Form("ratio_%s", dirName.Data()), Form("ratio_%s", dirName.Data())));
3200 tempRatio->Divide(4,2);
3201 TCanvas* tempResult(new TCanvas(Form("result_%s", dirName.Data()), Form("result_%s", dirName.Data())));
3202 tempResult->Divide(4,2);
3203 for(Int_t q(0); q < 8; q++) {
3204 tempIn = (TDirectoryFile*)tempDir->Get(Form("%s___%s", stringArray[q].Data(), dirName.Data()));
3205 if(tempIn) {
3206 // to see if the unfolding converged try to extract the pearson coefficients
3207 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
3208 if(pIn) {
3209 printf(" - %s in plane converged \n", dirName.Data());
3210 tempPearson->cd(1+q);
3211 Style(gPad, "PEARSON");
3212 pIn->DrawCopy("colz");
3213 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3214 if(rIn) {
3215 Style(rIn);
3216 printf(" > found RatioRefoldedMeasured < \n");
3217 tempRatio->cd(q+1);
3218 rIn->SetFillColor(kRed);
3219 rIn->Draw("ap");
3220 }
3221 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
3222 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
3223 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
3224 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
3225 if(dvector && avalue && rm && eff) {
3226 Style(dvector);
3227 Style(avalue);
3228 Style(rm);
3229 Style(eff);
3230 canvasMISC->cd(1);
3231 Style(gPad, "SPECTRUM");
3232 dvector->DrawCopy();
3233 canvasMISC->cd(2);
3234 Style(gPad, "SPECTRUM");
3235 avalue->DrawCopy();
3236 canvasMISC->cd(3);
3237 Style(gPad, "PEARSON");
3238 rm->DrawCopy("colz");
3239 canvasMISC->cd(4);
3240 Style(gPad, "GRID");
3241 eff->DrawCopy();
3242 } else if(rm && eff) {
3243 Style(rm);
3244 Style(eff);
3245 canvasMISC->cd(3);
3246 Style(gPad, "PEARSON");
3247 rm->DrawCopy("colz");
3248 canvasMISC->cd(4);
3249 Style(gPad, "GRID");
3250 eff->DrawCopy();
3251 }
3252 }
3253 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
3254 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
3255 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
3256 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3257 if(defaultUnfoldedJetSpectrumIn) {
3258 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
3259 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
3260 temp->Divide(unfoldedSpectrum);
3261 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3262 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3263 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3264 canvasMasterIn->cd(j);
3265 temp->GetYaxis()->SetRangeUser(0., 2);
3266 temp->DrawCopy();
3267 }
3268 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
3269 tempResult->cd(q+1);
3270 Style(gPad);
3271 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3272 unfoldedSpectrum->DrawCopy();
3273 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3274 inputSpectrum->DrawCopy("same");
3275 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3276 refoldedSpectrum->DrawCopy("same");
3277 TLegend* l(AddLegend(gPad));
3278 Style(l);
3279 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3280 Float_t chi(fitStatus->GetBinContent(1));
3281 Float_t pen(fitStatus->GetBinContent(2));
3282 Int_t dof((int)fitStatus->GetBinContent(3));
3283 Float_t beta(fitStatus->GetBinContent(4));
3284 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3285 } else if (fitStatus) { // only available in SVD fit
3286 Int_t reg((int)fitStatus->GetBinContent(1));
3287 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3288 }
3289 }
3290 }
3291 }
3292 }
3293 if(tempIn) {
3294 // to see if the unfolding converged try to extract the pearson coefficients
3295 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
3296 if(pIn) {
3297 printf(" - %s in plane converged \n", dirName.Data());
3298 canvasIn->cd(j);
3299 Style(gPad, "PEARSON");
3300 pIn->DrawCopy("colz");
3301 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3302 if(rIn) {
3303 Style(rIn);
3304 printf(" > found RatioRefoldedMeasured < \n");
3305 canvasRatioMeasuredRefoldedIn->cd(j);
3306 rIn->SetFillColor(kRed);
3307 rIn->Draw("ap");
3308 }
3309 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
3310 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
3311 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
3312 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
3313 if(dvector && avalue && rm && eff) {
3314 Style(dvector);
3315 Style(avalue);
3316 Style(rm);
3317 Style(eff);
3318 canvasMISC->cd(1);
3319 Style(gPad, "SPECTRUM");
3320 dvector->DrawCopy();
3321 canvasMISC->cd(2);
3322 Style(gPad, "SPECTRUM");
3323 avalue->DrawCopy();
3324 canvasMISC->cd(3);
3325 Style(gPad, "PEARSON");
3326 rm->DrawCopy("colz");
3327 canvasMISC->cd(4);
3328 Style(gPad, "GRID");
3329 eff->DrawCopy();
3330 } else if(rm && eff) {
3331 Style(rm);
3332 Style(eff);
3333 canvasMISC->cd(3);
3334 Style(gPad, "PEARSON");
3335 rm->DrawCopy("colz");
3336 canvasMISC->cd(4);
3337 Style(gPad, "GRID");
3338 eff->DrawCopy();
3339 }
3340 }
3341 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
3342 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
3343 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
3344 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3345 if(defaultUnfoldedJetSpectrumIn) {
3346 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
3347 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
3348 temp->Divide(unfoldedSpectrum);
3349 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3350 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3351 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3352 canvasMasterIn->cd(j);
3353 temp->GetYaxis()->SetRangeUser(0., 2);
3354 temp->DrawCopy();
3355 }
3356 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
3357 canvasSpectraIn->cd(j);
3358 Style(gPad);
3359 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3360 unfoldedSpectrum->DrawCopy();
3361 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3362 inputSpectrum->DrawCopy("same");
3363 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3364 refoldedSpectrum->DrawCopy("same");
3365 TLegend* l(AddLegend(gPad));
3366 Style(l);
3367 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3368 Float_t chi(fitStatus->GetBinContent(1));
3369 Float_t pen(fitStatus->GetBinContent(2));
3370 Int_t dof((int)fitStatus->GetBinContent(3));
3371 Float_t beta(fitStatus->GetBinContent(4));
3372 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3373 } else if (fitStatus) { // only available in SVD fit
3374 Int_t reg((int)fitStatus->GetBinContent(1));
3375 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3376 }
3377 }
3378 }
3379 if(tempOut) {
3380 TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
3381 if(pOut) {
3382 printf(" - %s out of plane converged \n", dirName.Data());
3383 canvasOut->cd(j);
3384 Style(gPad, "PEARSON");
3385 pOut->DrawCopy("colz");
3386 TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3387 if(rOut) {
3388 Style(rOut);
3389 printf(" > found RatioRefoldedMeasured < \n");
3390 canvasRatioMeasuredRefoldedOut->cd(j);
3391 rOut->SetFillColor(kRed);
3392 rOut->Draw("ap");
3393 }
3394 TH1D* dvector((TH1D*)tempOut->Get("dVector"));
3395 TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
3396 TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
3397 TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
3398 if(dvector && avalue && rm && eff) {
3399 Style(dvector);
3400 Style(avalue);
3401 Style(rm);
3402 Style(eff);
3403 canvasMISC->cd(5);
3404 Style(gPad, "SPECTRUM");
3405 dvector->DrawCopy();
3406 canvasMISC->cd(6);
3407 Style(gPad, "SPECTRUM");
3408 avalue->DrawCopy();
3409 canvasMISC->cd(7);
3410 Style(gPad, "PEARSON");
3411 rm->DrawCopy("colz");
3412 canvasMISC->cd(8);
3413 Style(gPad, "GRID");
3414 eff->DrawCopy();
3415 } else if(rm && eff) {
3416 Style(rm);
3417 Style(eff);
3418 canvasMISC->cd(7);
3419 Style(gPad, "PEARSON");
3420 rm->DrawCopy("colz");
3421 canvasMISC->cd(8);
3422 Style(gPad, "GRID");
3423 eff->DrawCopy();
3424 }
3425 }
3426 TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
3427 TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
3428 TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
3429 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3430 if(defaultUnfoldedJetSpectrumOut) {
3431 Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
3432 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirName.Data())));
3433 temp->Divide(unfoldedSpectrum);
3434 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3435 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3436 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3437 canvasMasterOut->cd(j);
3438 temp->GetYaxis()->SetRangeUser(0., 2.);
3439 temp->DrawCopy();
3440 }
3441 TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
3442 canvasSpectraOut->cd(j);
3443 Style(gPad);
3444 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3445 unfoldedSpectrum->DrawCopy();
3446 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3447 inputSpectrum->DrawCopy("same");
3448 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3449 refoldedSpectrum->DrawCopy("same");
3450 TLegend* l(AddLegend(gPad));
3451 Style(l);
3452 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3453 Float_t chi(fitStatus->GetBinContent(1));
3454 Float_t pen(fitStatus->GetBinContent(2));
3455 Int_t dof((int)fitStatus->GetBinContent(3));
3456 Float_t beta(fitStatus->GetBinContent(4));
3457 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3458 } else if (fitStatus) { // only available in SVD fit
3459 Int_t reg((int)fitStatus->GetBinContent(1));
3460 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3461 }
3462 }
3463 }
3464 if(canvasRatio && canvasV2) {
3465 canvasRatio->cd(j);
3466 TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
3467 Double_t _tempx(0), _tempy(0);
3468 if(ratioYield) {
3469 Style(ratioYield);
3470 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3471 ratioYield->GetPoint(b, _tempx, _tempy);
3472 ratioProfile->Fill(_tempx, _tempy);
3473 }
3474 ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
3475 ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3476 ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
3477 ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3478 ratioYield->SetFillColor(kRed);
3479 ratioYield->Draw("ap");
3480 }
3481 canvasV2->cd(j);
3482 TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
3483 if(ratioV2) {
3484 Style(ratioV2);
3485 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3486 ratioV2->GetPoint(b, _tempx, _tempy);
3487 v2Profile->Fill(_tempx, _tempy);
3488
3489 }
3490 v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
3491 v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3492 ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
3493 ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3494 ratioV2->SetFillColor(kRed);
3495 ratioV2->Draw("ap");
3496 }
3497 }
3498 }
3499 TFile output(out.Data(), "RECREATE");
3500 canvasProfiles->cd(1);
3501 Style(ratioProfile);
3502 ratioProfile->DrawCopy();
3503 canvasProfiles->cd(2);
3504 Style(v2Profile);
3505 v2Profile->DrawCopy();
3506 SavePadToPDF(canvasProfiles);
3507 canvasProfiles->Write();
3508 SavePadToPDF(canvasIn);
3509 canvasIn->Write();
3510 if(canvasOut) {
3511 SavePadToPDF(canvasOut);
3512 canvasOut->Write();
3513 }
3514 SavePadToPDF(canvasRatioMeasuredRefoldedIn);
3515 canvasRatioMeasuredRefoldedIn->Write();
3516 if(canvasRatioMeasuredRefoldedOut) {
3517 SavePadToPDF(canvasRatioMeasuredRefoldedOut);
3518 canvasRatioMeasuredRefoldedOut->Write();
3519 }
3520 SavePadToPDF(canvasSpectraIn);
3521 canvasSpectraIn->Write();
3522 if(canvasSpectraOut) {
3523 SavePadToPDF(canvasSpectraOut);
3524 canvasSpectraOut->Write();
3525 SavePadToPDF(canvasRatio);
3526 canvasRatio->Write();
3527 SavePadToPDF(canvasV2);
3528 canvasV2->Write();
3529 }
3530 SavePadToPDF(canvasMasterIn);
3531 canvasMasterIn->Write();
3532 if(canvasMasterOut) {
3533 SavePadToPDF(canvasMasterOut);
3534 canvasMasterOut->Write();
3535 }
3536 SavePadToPDF(canvasMISC);
3537 canvasMISC->Write();
3538 output.Write();
3539 output.Close();
3540}
3541//_____________________________________________________________________________
3542void AliJetFlowTools::BootstrapSpectra(TString def, TString in, TString out) const
3543{
3544 // function to interpret results of bootstrapping routine
3545 // TString def should hold the true emperical distribution
3546 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
3547 TFile readMe(in.Data(), "READ"); // open file read-only
3548 if(readMe.IsZombie()) {
3549 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
3550 return;
3551 }
3552 printf("\n\n\n\t\t BootstrapSpectra \n > Recovered the following file structure : \n <");
3553 readMe.ls();
3554 TList* listOfKeys((TList*)readMe.GetListOfKeys());
3555 if(!listOfKeys) {
3556 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
3557 return;
3558 }
3559 // get an estimate of the number of outputs and find the default set
3560 TDirectoryFile* defDir(0x0);
3561 for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
3562 if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
3563 if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3564 }
3565 }
3566
3567 // extract the default output this is the 'emperical' distribution
3568 // w.r.t. which the other distributions will be evaluated
3569 TH1D* defaultUnfoldedJetSpectrumIn(0x0);
3570 TH1D* defaultUnfoldedJetSpectrumOut(0x0);
3571 TH1D* defaultInputSpectrumIn(0x0);
3572 TH1D* defaultInputSpectrumOut(0x0);
3573 TGraphErrors* v2Emperical(0x0);
3574 if(defDir) {
3575 TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
3576 TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
3577 if(defDirIn) {
3578 defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
3579 defaultInputSpectrumIn = (TH1D*)defDirIn->Get(Form("InputSpectrum_in_%s", def.Data()));
3580 }
3581 if(defDirOut) {
3582 defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
3583 defaultInputSpectrumOut = (TH1D*)defDirOut->Get(Form("InputSpectrum_out_%s", def.Data()));
3584 }
3585 }
3586 if((!defaultUnfoldedJetSpectrumIn && defaultUnfoldedJetSpectrumOut)) {
3587 printf(" BootstrapSpectra: couldn't extract default spectra, aborting! \n");
3588 return;
3589 }
3590 else v2Emperical = GetV2(defaultUnfoldedJetSpectrumIn, defaultUnfoldedJetSpectrumOut, fEventPlaneRes);
3591
3592 // now that we know for sure that the input is in place, reserve the bookkeeping histograms
3593 TH1F* delta0[fBinsTrue->GetSize()-1];
3594 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3595 delta0[i] = new TH1F(Form("delta%i_%.2f-%.2f_gev", i, fBinsTrue->At(i), fBinsTrue->At(i+1)), Form("#Delta_{0, %i} p_{T} %.2f-%.2f GeV/c", i, fBinsTrue->At(i), fBinsTrue->At(i+1)), 30, -1., 1.);//.15, .15);
3596 delta0[i]->Sumw2();
3597 }
3598 // and a canvas for illustration purposes only
3599 TCanvas* spectraIn(new TCanvas("spectraIn", "spectraIn"));
3600 TCanvas* spectraOut(new TCanvas("spectraOut", "spectraOut"));
3601 // common reference (in this case the generated v2)
3602 TF1* commonReference = new TF1("v2_strong_bkg_comb", "(x<=3)*.07+(x>3&&x<5)*(.07-(x-3)*.035)+(x>30&&x<40)*(x-30)*.005+(x>40)*.05", 0, 200);
3603
3604 // loop through the directories, only plot the graphs if the deconvolution converged
3605 TDirectoryFile* tempDir(0x0);
3606 TDirectoryFile* tempIn(0x0);
3607 TDirectoryFile* tempOut(0x0);
3608 TH1D* unfoldedSpectrumIn(0x0);
3609 TH1D* unfoldedSpectrumOut(0x0);
3610 TH1D* measuredSpectrumIn(0x0);
3611 TH1D* measuredSpectrumOut(0x0);
3612 for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
3613 // read the directory by its name
3614 tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3615 if(!tempDir) continue;
3616 TString dirName(tempDir->GetName());
3617 // read only bootstrapped outputs
3618 if(!dirName.Contains("bootstrap")) continue;
3619 // try to read the in- and out of plane subdirs
3620 tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
3621 tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
3622 j++;
3623 // extract the unfolded spectra only if both in- and out-of-plane converted (i.e. pearson coefficients were saved)
3624 if(tempIn) {
3625 if(!(TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data()))) continue;
3626 unfoldedSpectrumIn = (TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data()));
3627 measuredSpectrumIn = (TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data()));
3628 spectraIn->cd();
3629 (j==1) ? measuredSpectrumIn->DrawCopy() : measuredSpectrumIn->DrawCopy("same");
3630 }
3631 if(tempOut) {
3632 if(!(TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data()))) continue;
3633 unfoldedSpectrumOut = (TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data()));
3634 measuredSpectrumOut = (TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data()));
3635 spectraOut->cd();
3636 (j==1) ? measuredSpectrumOut->DrawCopy() : measuredSpectrumOut->DrawCopy("same");
3637 }
3638 // get v2 with statistical uncertainties from the extracted spectra
3639 TGraphErrors* v2Bootstrapped(GetV2(unfoldedSpectrumIn, unfoldedSpectrumOut, fEventPlaneRes));
3640 // and loop over all bins to fill the bookkeeping histograms
3641 Double_t yBoot(0), yEmp(0), xDud(0);
3642 // optional: common reference (in this case the sampled v2 value)
3643 for(Int_t k(0); k < fBinsTrue->GetSize()-1; k++) {
3644 // read values point by point (passed by reference)
3645 v2Bootstrapped->GetPoint(k, xDud, yBoot);
3646 v2Emperical->GetPoint(k, xDud, yEmp);
3647 if(commonReference) {
3648 if(!commonReference->Eval(xDud)<1e-10) {
3649 yEmp/=commonReference->Eval(xDud);
3650 yBoot/=commonReference->Eval(xDud);
3651 } else { // if reference equals zero, take emperical distribution as reference
3652 yEmp/=yEmp; // 1
3653 yBoot/=yEmp;
3654 }
3655 }
3656 cout << " yEmp " << yEmp << " yBoot " << yBoot << endl;
3657 // save relative difference per pt bin
3658 if(TMath::Abs(yBoot)>1e-10) delta0[k]->Fill(yEmp - yBoot);
3659 }
3660 }
3661 // extracting final results now, as first estimate just a gaus fit to the distributions
3662 // (should be changed perhaps to proper rms eventually)
3663 // attach relevant data to current buffer in the same loop
3664 TFile output(out.Data(), "RECREATE");
3665 defaultInputSpectrumIn->SetLineColor(kRed);
3666 spectraIn->cd();
3667 defaultInputSpectrumIn->DrawCopy("same");
3668 defaultInputSpectrumOut->SetLineColor(kRed);
3669 spectraOut->cd();
3670 defaultInputSpectrumOut->DrawCopy("same");
3671 spectraIn->Write();
3672 spectraOut->Write();
3673
3674 TH1F* delta0spread(new TH1F("delta0spread", "#sigma(#Delta_{0})", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3675 TH1F* unfoldingError(new TH1F("unfoldingError", "error from unfolding algorithm", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3676 TF1* gaus(new TF1("gaus", "gaus"/*[0]*exp(-0.5*((x-[1])/[2])**2)"*/, -1., 1.));
3677 Double_t xDud(0), yDud(0);
3678 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3679 delta0[i]->Fit(gaus);
3680 delta0[i]->GetYaxis()->SetTitle("counts");
3681 delta0[i]->GetXaxis()->SetTitle("(v_{2, jet}^{measured} - v_{2, jet}^{generated}) / input v_{2}");
3682 delta0spread->SetBinContent(i+1, gaus->GetParameter(1)); // mean of gaus
3683 delta0spread->SetBinError(i+1, gaus->GetParameter(2)); // sigma of gaus
3684 cout << " mean " << gaus->GetParameter(1) << endl;
3685 cout << " sigm " << gaus->GetParameter(2) << endl;
3686 delta0[i]->Write();
3687 v2Emperical->GetPoint(i, xDud, yDud);
3688 unfoldingError->SetBinContent(i+1, 1e-10/*gaus->GetParameter(1)*/);
3689 if(commonReference && !commonReference->Eval(xDud)<1e-10) unfoldingError->SetBinError(i+1, v2Emperical->GetErrorY(i)/(commonReference->Eval(xDud)));
3690 else if(yDud>10e-10) unfoldingError->SetBinError(i+1, v2Emperical->GetErrorY(i)/yDud);
3691 else unfoldingError->SetBinError(i+1, 0.);
3692 }
3693 delta0spread->GetXaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
3694 delta0spread->GetYaxis()->SetTitle("(mean v_{2, jet}^{measured} - v_{2, jet}^{generated}) / input v_{2}");
3695 delta0spread->Write();
3696 unfoldingError->GetXaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
3697 unfoldingError->GetYaxis()->SetTitle("(mean v_{2, jet}^{measured} - v_{2, jet}^{generated}) / input v_{2}");
3698 unfoldingError->Write();
3699 // write the buffer and close the file
3700 output.Write();
3701 output.Close();
3702}
3703//_____________________________________________________________________________
3704Bool_t AliJetFlowTools::SetRawInput (
3705 TH2D* detectorResponse, // detector response matrix
3706 TH1D* jetPtIn, // in plane jet spectrum
3707 TH1D* jetPtOut, // out of plane jet spectrum
3708 TH1D* dptIn, // in plane delta pt distribution
3709 TH1D* dptOut, // out of plane delta pt distribution
3710 Int_t eventCount) {
3711 // set input histograms manually
3712 fDetectorResponse = detectorResponse;
3713 fSpectrumIn = jetPtIn;
3714 fSpectrumOut = jetPtOut;
3715 fDptInDist = dptIn;
3716 fDptOutDist = dptOut;
3717 fRawInputProvided = kTRUE;
3718 // check if all data is provided
3719 if(!fDetectorResponse) {
3720 printf(" fDetectorResponse not found \n ");
3721 return kFALSE;
3722 }
3723 // check if the pt bin for true and rec have been set
3724 if(!fBinsTrue || !fBinsRec) {
3725 printf(" No true or rec bins set, please set binning ! \n");
3726 return kFALSE;
3727 }
3728 if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
3729 // procedures, these profiles will be nonsensical, user is responsible
3730 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3731 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3732 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3733 }
3734 // normalize spectra to event count if requested
3735 if(fNormalizeSpectra) {
3736 fEventCount = eventCount;
3737 if(fEventCount > 0) {
3738 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
3739 fSpectrumOut->Sumw2();
3740 fSpectrumIn->Scale(1./((double)fEventCount));
3741 fSpectrumOut->Scale(1./((double)fEventCount));
3742 }
3743 }
3744 if(!fNormalizeSpectra && fEventCount > 0) {
3745 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
3746 fSpectrumOut->Sumw2();
3747 fSpectrumIn->Scale(1./((double)fEventCount));
3748 fSpectrumOut->Scale(1./((double)fEventCount));
3749 }
3750 fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
3751 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
3752 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
3753 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
3754 fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
3755 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
3756 fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
3757 fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
3758
3759 return kTRUE;
3760}
3761//_____________________________________________________________________________
3762TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
3763{
3764 // return ratio of h1 / h2
3765 // histograms can have different binning. errors are propagated as uncorrelated
3766 if(!(h1 && h2) ) {
3767 printf(" GetRatio called with NULL argument(s) \n ");
3768 return 0x0;
3769 }
3770 Int_t j(0);
3771 TGraphErrors *gr = new TGraphErrors();
3772 gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3773 Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
3774 TH1* dud((TH1*)h1->Clone("dud"));
3775 dud->Sumw2();
3776 h1->Sumw2();
3777 h2->Sumw2();
3778 if(!dud->Divide(h1, h2)) {
3779 printf(" > ROOT failed to divide two histogams, dividing manually \n < ");
3780 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3781 binCent = h1->GetXaxis()->GetBinCenter(i);
3782 if(xmax > 0. && binCent > xmax) continue;
3783 j = h2->FindBin(binCent);
3784 binWidth = h1->GetXaxis()->GetBinWidth(i);
3785 if(h2->GetBinContent(j) > 0.) {
3786 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
3787 Double_t A = h1->GetBinError(i)/h1->GetBinContent(i);
3788 Double_t B = h2->GetBinError(i)/h2->GetBinContent(i);
3789 error2 = ratio*ratio*A*A+ratio*ratio*B*B;
3790 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
3791 gr->SetPoint(i-1, binCent, ratio);
3792 gr->SetPointError(i-1, 0.5*binWidth, error2);
3793 }
3794 }
3795 } else {
3796 printf( " > using ROOT to divide two histograms < \n");
3797 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3798 binCent = dud->GetXaxis()->GetBinCenter(i);
3799 if(xmax > 0. && binCent > xmax) continue;
3800 binWidth = dud->GetXaxis()->GetBinWidth(i);
3801 gr->SetPoint(i-1,binCent,dud->GetBinContent(i));
3802 gr->SetPointError(i-1, 0.5*binWidth,dud->GetBinError(i));
3803 }
3804 }
3805
3806 if(appendFit) {
3807 TF1* fit(new TF1("lin", "pol0", 10, 100));
3808 gr->Fit(fit);
3809 }
3810 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
3811 if(dud) delete dud;
3812 return gr;
3813}
3814//_____________________________________________________________________________
3815TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
3816{
3817 // get v2 from difference of in plane, out of plane yield
3818 // h1 must hold the in-plane yield, h2 holds the out of plane yield
3819 // r is the event plane resolution for the chosen centrality
3820 if(!(h1 && h2) ) {
3821 printf(" GetV2 called with NULL argument(s) \n ");
3822 return 0x0;
3823 }
3824 TGraphErrors *gr = new TGraphErrors();
3825 gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3826 Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
3827 Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
3828
3829 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3830 binCent = h1->GetXaxis()->GetBinCenter(i);
3831 binWidth = h1->GetXaxis()->GetBinWidth(i);
3832 if(h2->GetBinContent(i) > 0.) {
3833 in = h1->GetBinContent(i);
3834 ein = h1->GetBinError(i);
3835 out = h2->GetBinContent(i);
3836 eout = h2->GetBinError(i);
3837 ratio = pre*((in-out)/(in+out));
3838 error2 = (4.*out*out/(TMath::Power(in+out, 4)))*ein*ein+(4.*in*in/(TMath::Power(in+out, 4)))*eout*eout;
3839 error2 = error2*pre*pre;
3840 if(error2 > 0) error2 = TMath::Sqrt(error2);
3841 gr->SetPoint(i-1,binCent,ratio);
3842 gr->SetPointError(i-1,0.5*binWidth,error2);
3843 }
3844 }
3845 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
3846 return gr;
3847}
3848//_____________________________________________________________________________
3849TGraphAsymmErrors* AliJetFlowTools::GetV2WithSystematicErrors(
3850 TH1* h1, TH1* h2, Double_t r, TString name,
3851 TH1* relativeErrorInUp,
3852 TH1* relativeErrorInLow,
3853 TH1* relativeErrorOutUp,
3854 TH1* relativeErrorOutLow,
3855 Float_t rho) const
3856{
3857 // get v2 with asymmetric systematic error
3858 // note that this is ONLY the systematic error, no statistical error!
3859 // rho is the pearson correlation coefficient
3860 TGraphErrors* tempV2(GetV2(h1, h2, r, name));
3861 Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
3862 Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
3863 Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
3864 Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
3865 Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
3866 Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
3867 Double_t in(0.), out(0.), einUp(0.), einLow(0.), eoutUp(0.), eoutLow(0.), error2Up(0.), error2Low(0.);
3868 // loop through the bins and do error propagation
3869 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3870 // extract the absolute errors
3871 in = h1->GetBinContent(i+1);
3872 einUp = TMath::Abs(in*relativeErrorInUp->GetBinContent(i+1));
3873 einLow = TMath::Abs(in*relativeErrorInLow->GetBinContent(1+i));
3874 out = h2->GetBinContent(i+1);
3875 eoutUp = TMath::Abs(out*relativeErrorOutUp->GetBinContent(1+i));
3876 eoutLow = TMath::Abs(out*relativeErrorOutLow->GetBinContent(1+i));
3877 // get the error squared
3878 if(rho <= 0) {
3879 error2Up = TMath::Power(((r*4.)/(TMath::Pi())),-2.)*((4.*out*out/(TMath::Power(in+out, 4)))*einUp*einUp+(4.*in*in/(TMath::Power(in+out, 4)))*eoutLow*eoutLow);
3880 error2Low =TMath::Power(((r*4.)/(TMath::Pi())),-2.)*((4.*out*out/(TMath::Power(in+out, 4)))*einLow*einLow+(4.*in*in/(TMath::Power(in+out, 4)))*eoutUp*eoutUp);
3881 } else {
3882 error2Up = TMath::Power(((r*4.)/(TMath::Pi())),-2.)*((4.*out*out/(TMath::Power(in+out, 4)))*einUp*einUp+(4.*in*in/(TMath::Power(in+out, 4)))*eoutUp*eoutUp-((8.*out*in)/(TMath::Power(in+out, 4)))*rho*einUp*eoutUp);
3883 error2Low =TMath::Power(((r*4.)/(TMath::Pi())),-2.)*((4.*out*out/(TMath::Power(in+out, 4)))*einLow*einLow+(4.*in*in/(TMath::Power(in+out, 4)))*eoutLow*eoutLow-((8.*out*in)/(TMath::Power(in+out, 4)))*rho*einLow*eoutLow);
3884 }
3885 if(error2Up > 0) error2Up = TMath::Sqrt(error2Up);
3886 if(error2Low > 0) error2Low = TMath::Sqrt(error2Low);
3887 // set the errors
3888 ayh[i] = error2Up;
3889 ayl[i] = error2Low;
3890 // get the bin width (which is the 'error' on x)
3891 Double_t binWidth(h1->GetBinWidth(i+1));
3892 axl[i] = binWidth/2.;
3893 axh[i] = binWidth/2.;
3894 // now get the coordinate for the poin
3895 tempV2->GetPoint(i, ax[i], ay[i]);
3896 }
3897 // save the nominal ratio
3898 TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
3899 nominalError->SetName("v_{2} with shape uncertainty");
3900 // do some memory management
3901 delete tempV2;
3902 delete[] ax;
3903 delete[] ay;
3904 delete[] axh;
3905 delete[] axl;
3906 delete[] ayh;
3907 delete[] ayl;
3908
3909 return nominalError;
3910}
3911//_____________________________________________________________________________
3912void AliJetFlowTools::WriteObject(TObject* object, TString suffix, Bool_t kill) {
3913 // write object, if a unique identifier is given the object is cloned
3914 // and the clone is saved. setting kill to true will delete the original obect from the heap
3915 if(!object) {
3916 printf(" > WriteObject:: called with NULL arguments \n ");
3917 return;
3918 } else if(!strcmp("", suffix.Data())) object->Write();
3919 else {
3920 TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
3921 newObject->Write();
3922 }
3923 if(kill) delete object;
3924}
3925//_____________________________________________________________________________
3926TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
3927 // construt a delta pt response matrix from supplied dpt distribution
3928 // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
3929 // do a weighted rebinning to a (coarser) dpt distribution
3930 // be careful with the binning of the dpt response: it should be equal to that
3931 // of the response matrix, otherwise dpt and response matrices cannot be multiplied
3932 //
3933 // the response matrix will be square and have the same binning
3934 // (min, max and granularity) of the input histogram
3935 Int_t bins(dpt->GetXaxis()->GetNbins()); // number of bins, will also be no of rows, columns
3936 Double_t _bins[bins+1]; // prepare array with bin borders
3937 for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
3938 _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1); // get upper edge
3939 TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
3940 for(Int_t j(0); j < bins+1; j++) { // loop on pt true slices j
3941 Bool_t skip = kFALSE;
3942 for(Int_t k(0); k < bins+1; k++) { // loop on pt gen slices k
3943 (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
3944 if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
3945 }
3946 }
3947 return res;
3948}
3949//_____________________________________________________________________________
3950TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
3951 if(!binsTrue || !binsRec) {
3952 printf(" > GetUnityResponse:: function called with NULL arguments < \n");
3953 return 0x0;
3954 }
3955 TString name(Form("unityResponse_%s", suffix.Data()));
3956 TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
3957 for(Int_t i(0); i < binsTrue->GetSize(); i++) {
3958 for(Int_t j(0); j < binsRec->GetSize(); j++) {
3959 if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
3960 }
3961 }
3962 return unity;
3963}
3964//_____________________________________________________________________________
3965void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
3966 // save configuration parameters to histogram
3967 TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
3968 summary->SetBinContent(1, fBetaIn);
3969 summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
3970 summary->SetBinContent(2, fBetaOut);
3971 summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
3972 summary->SetBinContent(3, fCentralityArray->At(0));
3973 summary->GetXaxis()->SetBinLabel(3, "fCentralityArray[0]");
3974 summary->SetBinContent(4, (int)convergedIn);
3975 summary->GetXaxis()->SetBinLabel(4, "convergedIn");
3976 summary->SetBinContent(5, (int)convergedOut);
3977 summary->GetXaxis()->SetBinLabel(5, "convergedOut");
3978 summary->SetBinContent(6, (int)fAvoidRoundingError);
3979 summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
3980 summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
3981 summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
3982 summary->SetBinContent(8, (int)fPrior);
3983 summary->GetXaxis()->SetBinLabel(8, "fPrior");
3984 summary->SetBinContent(9, fSVDRegIn);
3985 summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
3986 summary->SetBinContent(10, fSVDRegOut);
3987 summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
3988 summary->SetBinContent(11, (int)fSVDToy);
3989 summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
3990 summary->SetBinContent(12, fJetRadius);
3991 summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
3992 summary->SetBinContent(13, (int)fNormalizeSpectra);
3993 summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
3994 summary->SetBinContent(14, (int)fSmoothenPrior);
3995 summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
3996 summary->SetBinContent(15, (int)fTestMode);
3997 summary->GetXaxis()->SetBinLabel(15, "fTestMode");
3998 summary->SetBinContent(16, (int)fUseDetectorResponse);
3999 summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
4000 summary->SetBinContent(17, fBayesianIterIn);
4001 summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
4002 summary->SetBinContent(18, fBayesianIterOut);
4003 summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
4004 summary->SetBinContent(19, fBayesianSmoothIn);
4005 summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
4006 summary->SetBinContent(20, fBayesianSmoothOut);
4007 summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
4008}
4009//_____________________________________________________________________________
4010void AliJetFlowTools::ResetAliUnfolding() {
4011 // ugly function: reset all unfolding parameters
4012 TVirtualFitter* fitter(TVirtualFitter::GetFitter());
4013 if(fitter) {
4014 printf(" > Found fitter, will delete it < \n");
4015 delete fitter;
4016 }
4017 if(gMinuit) {
4018 printf(" > Found gMinuit, will re-create it < \n");
4019 delete gMinuit;
4020 gMinuit = new TMinuit;
4021 }
4022 AliUnfolding::fgCorrelationMatrix = 0;
4023 AliUnfolding::fgCorrelationMatrixSquared = 0;
4024 AliUnfolding::fgCorrelationCovarianceMatrix = 0;
4025 AliUnfolding::fgCurrentESDVector = 0;
4026 AliUnfolding::fgEntropyAPriori = 0;
4027 AliUnfolding::fgEfficiency = 0;
4028 AliUnfolding::fgUnfoldedAxis = 0;
4029 AliUnfolding::fgMeasuredAxis = 0;
4030 AliUnfolding::fgFitFunction = 0;
4031 AliUnfolding::fgMaxInput = -1;
4032 AliUnfolding::fgMaxParams = -1;
4033 AliUnfolding::fgOverflowBinLimit = -1;
4034 AliUnfolding::fgRegularizationWeight = 10000;
4035 AliUnfolding::fgSkipBinsBegin = 0;
4036 AliUnfolding::fgMinuitStepSize = 0.1;
4037 AliUnfolding::fgMinuitPrecision = 1e-6;
4038 AliUnfolding::fgMinuitMaxIterations = 1000000;
4039 AliUnfolding::fgMinuitStrategy = 1.;
4040 AliUnfolding::fgMinimumInitialValue = kFALSE;
4041 AliUnfolding::fgMinimumInitialValueFix = -1;
4042 AliUnfolding::fgNormalizeInput = kFALSE;
4043 AliUnfolding::fgNotFoundEvents = 0;
4044 AliUnfolding::fgSkipBin0InChi2 = kFALSE;
4045 AliUnfolding::fgBayesianSmoothing = 1;
4046 AliUnfolding::fgBayesianIterations = 10;
4047 AliUnfolding::fgDebug = kFALSE;
4048 AliUnfolding::fgCallCount = 0;
4049 AliUnfolding::fgPowern = 5;
4050 AliUnfolding::fChi2FromFit = 0.;
4051 AliUnfolding::fPenaltyVal = 0.;
4052 AliUnfolding::fAvgResidual = 0.;
4053 AliUnfolding::fgPrintChi2Details = 0;
4054 AliUnfolding::fgCanvas = 0;
4055 AliUnfolding::fghUnfolded = 0;
4056 AliUnfolding::fghCorrelation = 0;
4057 AliUnfolding::fghEfficiency = 0;
4058 AliUnfolding::fghMeasured = 0;
4059 AliUnfolding::SetMinuitStepSize(1.);
4060 AliUnfolding::SetMinuitPrecision(1e-6);
4061 AliUnfolding::SetMinuitMaxIterations(100000);
4062 AliUnfolding::SetMinuitStrategy(2.);
4063 AliUnfolding::SetDebug(1);
4064}
4065//_____________________________________________________________________________
4066TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) const {
4067 // protect heap by adding unique qualifier to name
4068 if(!protect) return 0x0;
4069 TH1D* p = (TH1D*)protect->Clone();
4070 TString tempString(fActiveString);
4071 tempString+=suffix;
4072 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
4073 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4074 if(kill) delete protect;
4075 return p;
4076}
4077//_____________________________________________________________________________
4078TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) const {
4079 // protect heap by adding unique qualifier to name
4080 if(!protect) return 0x0;
4081 TH2D* p = (TH2D*)protect->Clone();
4082 TString tempString(fActiveString);
4083 tempString+=suffix;
4084 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
4085 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4086 if(kill) delete protect;
4087 return p;
4088}
4089//_____________________________________________________________________________
4090TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) const {
4091 // protect heap by adding unique qualifier to name
4092 if(!protect) return 0x0;
4093 TGraphErrors* p = (TGraphErrors*)protect->Clone();
4094 TString tempString(fActiveString);
4095 tempString+=suffix;
4096 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
4097 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4098 if(kill) delete protect;
4099 return p;
4100}
4101//_____________________________________________________________________________
4102void AliJetFlowTools::MakeAU() {
4103 // === azimuthal unfolding ===
4104 //
4105 // unfolds the spectrum in delta phi bins, extracts the yield per bin, and does a fit
4106 // in transverse momentum and azimuthal correlation space to extract v2 params
4107 // settings are equal to the ones used for 'Make()'
4108 //
4109 // basic steps that are followed:
4110 // 1) rebin the raw output of the jet task to the desired binnings
4111 // 2) calls the unfolding routine
4112 // 3) writes output to file
4113 // can be repeated multiple times with different configurations
4114
4115 Int_t low[] = {1, 6, 11, 16, 21, 26, 31, 36};
4116 Int_t up[] = {5, 10, 15, 20, 25, 30, 35, 40};
4117 TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
4118 const Int_t ptBins(fBinsTrue->GetSize()-1);
4119 const Int_t dPhiBins(8);
4120 TH1D* dPtdPhi[fBinsTrue->GetSize()];
4121 for(Int_t i(0); i < ptBins; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), dPhiBins, 0, TMath::Pi());
4122
4123 // for the output initialize a canvas
4124 TCanvas* v2Fits(new TCanvas("v2 fits", "v2 fits"));
4125 v2Fits->Divide(4, TMath::Floor((1+ptBins)/(float)4)+(((1+ptBins)%4)>0));
4126
4127 for(Int_t i(0); i < dPhiBins; i++) {
4128 // 1) manipulation of input histograms
4129 // check if the input variables are present
4130 if(!PrepareForUnfolding(low[i], up[i])) return;
4131 // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
4132 // parts of the spectrum can end up in over or underflow bins
4133 TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, Form("resized_%s", stringArray[i].Data()), kFALSE);
4134
4135 // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
4136 // the template will be used as a prior for the chi2 unfolding
4137 TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, stringArray[i], kFALSE);
4138
4139 // get the full response matrix from the dpt and the detector response
4140 fDetectorResponse = NormalizeTH2D(fDetectorResponse);
4141 // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
4142 // so that unfolding should return the initial spectrum
4143 if(fUseDptResponse && fUseDetectorResponse) fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
4144 else if (fUseDptResponse && !fUseDetectorResponse) fFullResponseIn = fDptIn;
4145 else if (!fUseDptResponse && fUseDetectorResponse) fFullResponseIn = fDetectorResponse;
4146 else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) return;
4147 // normalize each slide of the response to one
4148 NormalizeTH2D(fFullResponseIn);
4149 // resize to desired binning scheme
4150 TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, stringArray[i]);
4151 // get the kinematic efficiency
4152 TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
4153 kinematicEfficiencyIn->SetNameTitle(Form("kin_eff_%s", stringArray[i].Data()), Form("kin_eff_%s", stringArray[i].Data()));
4154 // suppress the errors
4155 for(Int_t j(0); j < kinematicEfficiencyIn->GetXaxis()->GetNbins(); j++) kinematicEfficiencyIn->SetBinError(1+j, 0.);
4156 TH1D* jetFindingEfficiency(0x0);
4157 if(fJetFindingEff) {
4158 jetFindingEfficiency = ProtectHeap(fJetFindingEff);
4159 jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
4160 jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
4161 }
4162 // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
4163 TH1D* unfoldedJetSpectrumIn(0x0);
4164 fActiveDir->cd(); // select active dir
4165 TDirectoryFile* dirIn = new TDirectoryFile(Form("%s___%s", stringArray[i].Data(), fActiveString.Data()), Form("%s___%s", stringArray[i].Data(), fActiveString.Data()));
4166 dirIn->cd(); // select inplane subdir
4167 // select the unfolding method
4168 unfoldedJetSpectrumIn = UnfoldWrapper(
4169 measuredJetSpectrumIn,
4170 resizedResponseIn,
4171 kinematicEfficiencyIn,
4172 measuredJetSpectrumTrueBinsIn,
4173 TString("dPtdPhiUnfolding"),
4174 jetFindingEfficiency);
4175 // arbitrarily save one of the full outputs (same for all dphi bins, avoid duplicates)
4176 if(i+1 == ptBins) {
4177 resizedResponseIn->SetNameTitle(Form("ResponseMatrix_%s", stringArray[i].Data()), Form("response matrix %s", stringArray[i].Data()));
4178 resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
4179 resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
4180 resizedResponseIn = ProtectHeap(resizedResponseIn);
4181 resizedResponseIn->Write();
4182 kinematicEfficiencyIn->SetNameTitle(Form("KinematicEfficiency_%s", stringArray[i].Data()), Form("Kinematic efficiency, %s", stringArray[i].Data()));
4183 kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
4184 kinematicEfficiencyIn->Write();
4185 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
4186 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
4187 fDetectorResponse->Write();
4188 // optional histograms
4189 if(fSaveFull) {
4190 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", Form("[INPUT] Jet spectrum, %s", stringArray[i].Data()));
4191 fSpectrumIn->Write();
4192 fDptInDist->SetNameTitle("[ORIG]DeltaPt", Form("#delta p_{T} distribution, %s", stringArray[i].Data()));
4193 fDptInDist->Write();
4194 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix", Form("#delta p_{T} matrix, %s", stringArray[i].Data()));
4195 fDptIn->Write();
4196 fFullResponseIn->SetNameTitle("ResponseMatrix", Form("Response matrix, %s", stringArray[i].Data()));
4197 fFullResponseIn->Write();
4198 }
4199 }
4200 fActiveDir->cd();
4201 fDeltaPtDeltaPhi->Write();
4202 fJetPtDeltaPhi->Write();
4203
4204 TH1D* dud(ProtectHeap(unfoldedJetSpectrumIn, kTRUE, stringArray[i]));;
4205 Double_t integralError(0);
4206 // at this point in the code, the spectrum has been unfolded in a certain region of dPhi space
4207 // next step is splitting it in pt space as well to estimate the yield differentially in pt
4208 for(Int_t j(0); j < ptBins; j++) {
4209 // get the integrated
4210 Double_t integral(dud->IntegralAndError(j+1, j+2, integralError));
4211 dPtdPhi[j]->SetBinContent(i+1, integral);
4212 dPtdPhi[j]->SetBinError(i+1, integralError);
4213 }
4214 dud->Write();
4215 // save the current state of the unfolding object
4216 SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, kFALSE);
4217 }
4218 TF1* fourier = new TF1("fourier", "[0]*(1.+0.5*[1]*(TMath::Cos(2.*x)))", 0, TMath::Pi());
4219 TH1D* v2(new TH1D("v2FromFit", "v2FromFit", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
4220 for(Int_t i(0); i < ptBins; i++) {
4221 v2Fits->cd(i+1);
4222 dPtdPhi[i]->Fit(fourier, "VI");
4223 Style(gPad, "PEARSON");
4224 Style(dPtdPhi[i], kBlue, kDeltaPhi);
4225 dPtdPhi[i]->DrawCopy();
4226 AliJetFlowTools::AddText(
4227 TString(Form("%.2f #LT p_{T} #LT %.2f", fBinsTrue->At(i), fBinsTrue->At(i+1))),
4228 kBlack,
4229 .38,
4230 .56,
4231 .62,
4232 .65
4233 );
4234 v2->SetBinContent(1+i, fourier->GetParameter(1));
4235 v2->SetBinError(1+i, fourier->GetParError(1));
4236 }
4237 v2Fits->cd(1+ptBins);
4238 Style(gPad, "PEARSON");
4239 Style(v2, kBlack, kV2, kTRUE);
4240 v2->DrawCopy();
4241 v2Fits->Write();
4242}
4243//_____________________________________________________________________________
4244void AliJetFlowTools::ReplaceBins(TArrayI* array, TGraphErrors* graph) {
4245 // replace bins
4246 Double_t x(0), y(0);
4247 graph->GetPoint(0, x, y);
4248 graph->SetPoint(array->At(0)-1, fBinsTrue->At(array->At(0)), y);
4249 graph->SetPointError(array->At(0)-1, 10, graph->GetErrorY(0));
4250 graph->SetPoint(array->At(1)-1, -5, -5);
4251}
4252//_____________________________________________________________________________
4253void AliJetFlowTools::ReplaceBins(TArrayI* array, TGraphAsymmErrors* graph) {
4254 // replace bins
4255 Double_t x(0), y(0);
4256 graph->GetPoint(0, x, y);
4257 graph->SetPoint(array->At(0)-1, fBinsTrue->At(array->At(0)), y);
4258 Double_t yl = graph->GetErrorYlow(0);
4259 Double_t yh = graph->GetErrorYhigh(0);
4260 graph->SetPointError(array->At(0)-1, 10, 10, yl, yh);
4261 graph->SetPoint(array->At(1)-1, -5, -5);
4262}
4263//_____________________________________________________________________________
4264void AliJetFlowTools::GetSignificance(
4265 TGraphErrors* n, // points with stat error
4266 TGraphAsymmErrors* shape, // points with shape error
4267 TGraphAsymmErrors* corr, // points with stat error
4268 Int_t low, // lower bin (tgraph starts at 0)
4269 Int_t up // upper bin
4270 )
4271{
4272 // main use of this function is filling the static buffers
4273 Double_t statE(0), shapeE(0), corrE(0), y(0), x(0), chi2(0);
4274
4275 // print some stuff
4276 printf(" double v2[] = {\n");
4277 Int_t iterator(0);
4278 for(Int_t i(low); i < up+1; i++) {
4279 n->GetPoint(i, x, y);
4280 if(i==up) printf("%.4f}; \n\n", y);
4281 else printf("%.4f, \n", y);
4282 gV2->SetAt(y, iterator);
4283 iterator++;
4284 }
4285 iterator = 0;
4286 printf(" double stat[] = {\n");
4287 for(Int_t i(low); i < up+1; i++) {
4288 y = n->GetErrorYlow(i);
4289 if(i==up) printf("%.4f}; \n\n", y);
4290 else printf("%.4f, \n", y);
4291 gStat->SetAt(y, iterator);
4292 iterator++;
4293 }
4294 iterator = 0;
4295 printf(" double shape[] = {\n");
4296 for(Int_t i(low); i < up+1; i++) {
4297 y = shape->GetErrorYhigh(i);
4298 y*=gReductionFactor;
4299 shape->SetPointEYhigh(i, y);
4300 y = shape->GetErrorYlow(i);
4301 y*=gReductionFactor;
4302 shape->SetPointEYlow(i, y);
4303 if(i==up) printf("%.4f}; \n\n", y);
4304 else printf("%.4f, \n", y);
4305 gShape->SetAt(y, iterator);
4306 iterator++;
4307 }
4308 iterator = 0;
4309 printf(" double corr[] = {\n");
4310 for(Int_t i(low); i < up+1; i++) {
4311 y = corr->GetErrorYlow(i);
4312 if(i==up) printf("%.4f}; \n\n", y);
4313 else printf("%.4f, \n", y);
4314 gCorr->SetAt(y, iterator);
4315 iterator++;
4316 }
4317
4318 // to plot the average error as function of number of events
4319 Float_t ctr(0);
4320 for(Int_t i(low); i < up+1; i++) {
4321 ctr = ctr + 1.;
4322 // set some flags to 0
4323 x = 0.;
4324 y = 0.;
4325 // get the nominal point
4326 n->GetPoint(i, x, y);
4327 statE += n->GetErrorY(i);
4328 shapeE += shape->GetErrorY(i);
4329 corrE += corr->GetErrorY(i);
4330 // combine the errors
4331 }
4332 printf(" ======================================\n");
4333 printf(" > between %i and %i GeV/c \n", low, up);
4334 cout << " AVERAGE_SHAPE " << shapeE/ctr << endl;
4335 cout << " AVERAGE_CORR " << corrE/ctr << endl;
4336 cout << " AVERAGE_STAT " << statE/ctr << endl;
4337}
4338//_____________________________________________________________________________
4339void AliJetFlowTools::MinimizeChi2nd()
4340{
4341 // Choose method upon creation between:
4342 // kMigrad, kSimplex, kCombined,
4343 // kScan, kFumili
4344 ROOT::Minuit2::Minuit2Minimizer min ( ROOT::Minuit2::kMigrad );
4345 min.SetMaxFunctionCalls(1000000);
4346 min.SetMaxIterations(100000);
4347 min.SetTolerance(0.001);
4348
4349 ROOT::Math::Functor f(&PhenixChi2nd,2);
4350 double step[] = {0.0000001, 0.0000001};
4351 double variable[] = {-1., -1.};
4352
4353 min.SetFunction(f);
4354 // Set the free variables to be minimized!
4355 min.SetVariable(0,"epsilon_c",variable[0], step[0]);
4356 min.SetVariable(1,"epsilon_b",variable[1], step[1]);
4357
4358
4359 min.Minimize();
4360}
4361//_____________________________________________________________________________
4362Double_t AliJetFlowTools::PhenixChi2nd(const Double_t *xx )
4363{
4364 // define arrays with results and errors here, see example at PhenixChi2()
4365
4366 // return the function value at certain epsilon
4367 const Double_t epsc = xx[0];
4368 const Double_t epsb = xx[1];
4369 Double_t chi2(0);
4370 Int_t counts(gV2->GetSize() + gOffsetStop);
4371
4372 // altered implemtation of eq 3 of arXiv:0801.1665v2
4373 // see analysis note and QM2014 poster for validation
4374 for(Int_t i(gOffsetStart); i < counts; i++) {
4375
4376 // quadratic sum of statistical and uncorrelated systematic error
4377 Double_t e = gStat->At(i);;
4378
4379 // sum of v2 plus epsilon times correlated error minus hypothesis (0)
4380 // also the numerator of equation 3 of phenix paper
4381 Double_t numerator = TMath::Power(gV2->At(i)+epsc*gCorr->At(i)+epsb, 2);
4382
4383 // modified denominator of equation 3 of phenix paper
4384 Double_t denominator = e*e;
4385
4386 // add to the sum
4387 chi2 += numerator/denominator;
4388 }
4389 // add the square of epsilon to the total chi2 as penalty
4390
4391 Double_t sumEpsb(0);
4392 for(Int_t j(gOffsetStart); j < counts; j++) sumEpsb += (epsb*epsb)/(gShape->At(j)*gShape->At(j));
4393 chi2 += epsc*epsc + sumEpsb/((float)counts);
4394
4395 return chi2;
4396}
4397//_____________________________________________________________________________
4398Double_t AliJetFlowTools::ConstructFunctionnd(Double_t *x, Double_t *par)
4399{
4400 // internal use only: evaluate the function at a given point
4401 return AliJetFlowTools::PhenixChi2nd(x);
4402}
4403//_____________________________________________________________________________
4404TF2* AliJetFlowTools::ReturnFunctionnd(Double_t &p)
4405{
4406 // return the fitting function, pass the p-value w.r.t. 0 by reference
4407 const Int_t DOF(4);
4408 TF2 *f1 = new TF2("ndhist",AliJetFlowTools::ConstructFunctionnd, -100, 100, -100, 100, 0);
4409 printf(" > locating minima < \n");
4410 Double_t x(0), y(0);
4411 f1->GetMinimumXY(x, y);
4412 f1->GetXaxis()->SetTitle("#epsilon{b}");
4413 f1->GetXaxis()->SetTitle("#epsilon_{c}");
4414 f1->GetZaxis()->SetTitle("#chi^{2}");
4415
4416 printf(" ===============================================================================\n");
4417 printf(" > minimal chi2 f(%.8f, %.8f) = %.8f (i should be ok ... ) \n", x, y, f1->Eval(x, y));
4418 cout << " so the probability of finding data at least as imcompatible with 0 as the actually" << endl;
4419 cout << " observed data is " << TMath::Prob(f1->Eval(x, y), DOF) << endl;
4420 printf(" ===============================================================================\n");
4421
4422 // pass the p-value by reference and return the function
4423 p = TMath::Prob(f1->Eval(x, y), DOF);
4424 return f1;
4425}
4426//_____________________________________________________________________________