bugfix in function which selects unfolding routine
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliJetFlowTools.cxx
CommitLineData
dc1455ee 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)
53547ff2 18//
dc1455ee 19// Tools class for Jet Flow Analysis, replaces 'extractJetFlow.C' macro
4292ca60 20//
21// The task uses input from two analysis tasks:
0396b3c6 22// - $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetV2.cxx
4292ca60 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
5e11c41c 28// libraries must be present on the system
29// ( see http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html ).
4292ca60 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
dc1455ee 38
39// root includes
e20bfd76 40#include "TH1.h"
e4b6d21a 41#include "TF2.h"
dc1455ee 42#include "TH2D.h"
53547ff2 43#include "TGraph.h"
4292ca60 44#include "TGraphErrors.h"
a39e4b2b 45#include "TGraphAsymmErrors.h"
18698978 46#include "TLine.h"
d7ec324f 47#include "TCanvas.h"
48#include "TLegend.h"
dc1455ee 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"
dc1455ee 56#include "TLine.h"
ad04a83c 57#include "TMath.h"
4292ca60 58#include "TVirtualFitter.h"
ef12d5a5 59#include "TFitResultPtr.h"
1e01bfd6 60#include "Minuit2/Minuit2Minimizer.h"
61#include "Math/Functor.h"
4292ca60 62// aliroot includes
dc1455ee 63#include "AliUnfolding.h"
64#include "AliAnaChargedJetResponseMaker.h"
4292ca60 65// class includes
dc1455ee 66#include "AliJetFlowTools.h"
51e6bc5a 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"
549b5f40 71#include "RooUnfoldBayes.h"
51e6bc5a 72#include "TSVDUnfold.h"
dc1455ee 73
74using namespace std;
dc1455ee 75//_____________________________________________________________________________
76AliJetFlowTools::AliJetFlowTools() :
4292ca60 77 fResponseMaker (new AliAnaChargedJetResponseMaker()),
3e698d27 78 fRMS (kTRUE),
79 fSymmRMS (kTRUE),
6c3fa11d 80 fRho0 (kFALSE),
8bc7d79b 81 fBootstrap (kFALSE),
ef12d5a5 82 fPower (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
9f892925 83 fSaveFull (kTRUE),
dc1455ee 84 fActiveString (""),
ad04a83c 85 fActiveDir (0x0),
dc1455ee 86 fInputList (0x0),
4292ca60 87 fRefreshInput (kTRUE),
dc1455ee 88 fOutputFileName ("UnfoldedSpectra.root"),
ad04a83c 89 fOutputFile (0x0),
67d11165 90 fCentralityArray (0x0),
1e01bfd6 91 fMergeBinsArray (0x0),
c03f7598 92 fCentralityWeights (0x0),
dc1455ee 93 fDetectorResponse (0x0),
53547ff2 94 fJetFindingEff (0x0),
4292ca60 95 fBetaIn (.1),
96 fBetaOut (.1),
549b5f40
RAB
97 fBayesianIterIn (4),
98 fBayesianIterOut (4),
99 fBayesianSmoothIn (0.),
100 fBayesianSmoothOut (0.),
ef12d5a5 101 fAvoidRoundingError (kFALSE),
51e6bc5a 102 fUnfoldingAlgorithm (kChi2),
103 fPrior (kPriorMeasured),
18698978 104 fPriorUser (0x0),
dc1455ee 105 fBinsTrue (0x0),
106 fBinsRec (0x0),
ef12d5a5 107 fBinsTruePrior (0x0),
108 fBinsRecPrior (0x0),
4292ca60 109 fSVDRegIn (5),
110 fSVDRegOut (5),
51e6bc5a 111 fSVDToy (kTRUE),
112 fJetRadius (0.3),
20abfcc4 113 fEventCount (-1),
549b5f40
RAB
114 fNormalizeSpectra (kFALSE),
115 fSmoothenPrior (kFALSE),
4292ca60 116 fFitMin (60.),
549b5f40 117 fFitMax (300.),
4292ca60 118 fFitStart (75.),
549b5f40 119 fSmoothenCounts (kTRUE),
ef12d5a5 120 fTestMode (kFALSE),
4292ca60 121 fRawInputProvided (kFALSE),
ef12d5a5 122 fEventPlaneRes (.63),
123 fUseDetectorResponse(kTRUE),
549b5f40 124 fUseDptResponse (kTRUE),
ef12d5a5 125 fTrainPower (kTRUE),
486fb24e 126 fDphiUnfolding (kTRUE),
127 fDphiDptUnfolding (kFALSE),
128 fExLJDpt (kTRUE),
18698978 129 fTitleFontSize (-999.),
4292ca60 130 fRMSSpectrumIn (0x0),
131 fRMSSpectrumOut (0x0),
132 fRMSRatio (0x0),
ef12d5a5 133 fRMSV2 (0x0),
ad04a83c 134 fDeltaPtDeltaPhi (0x0),
135 fJetPtDeltaPhi (0x0),
dc1455ee 136 fSpectrumIn (0x0),
137 fSpectrumOut (0x0),
138 fDptInDist (0x0),
139 fDptOutDist (0x0),
140 fDptIn (0x0),
141 fDptOut (0x0),
142 fFullResponseIn (0x0),
549b5f40 143 fFullResponseOut (0x0) { // class constructor
486fb24e 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.);
4292ca60 147}
dc1455ee 148//_____________________________________________________________________________
e20bfd76 149void AliJetFlowTools::Make(TH1* customIn, TH1* customOut) {
ad04a83c 150 // core function of the class
486fb24e 151 if(fDphiDptUnfolding) {
152 // to extract the yield as function of Dphi, Dpt - experimental
153 MakeAU();
154 return;
155 }
4292ca60 156 // 1) rebin the raw output of the jet task to the desired binnings
ad04a83c 157 // 2) calls the unfolding routine
158 // 3) writes output to file
4292ca60 159 // can be repeated multiple times with different configurations
160
ad04a83c 161 // 1) manipulation of input histograms
dc1455ee 162 // check if the input variables are present
4292ca60 163 if(fRefreshInput) {
e20bfd76 164 if(!PrepareForUnfolding(customIn, customOut)) {
4292ca60 165 printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
166 return;
167 }
dc1455ee 168 }
4292ca60 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
8bc7d79b 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 ... )
a98c2619 175/*
8bc7d79b 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 }
a98c2619 181*/
549b5f40 182 TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
f3ba6c8e 183 TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec, TString("resized_out_"), kFALSE);
a98c2619 184
8bc7d79b 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)
a98c2619 190
8bc7d79b 191
f3ba6c8e 192
549b5f40 193 // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
4292ca60 194 // the template will be used as a prior for the chi2 unfolding
549b5f40
RAB
195 TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
196 TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
4292ca60 197 // get the full response matrix from the dpt and the detector response
dc1455ee 198 fDetectorResponse = NormalizeTH2D(fDetectorResponse);
ef12d5a5 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) {
549b5f40
RAB
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 }
ef12d5a5 215 } else {
216 fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
217 fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
218 }
4292ca60 219 // normalize each slide of the response to one
dc1455ee 220 NormalizeTH2D(fFullResponseIn);
221 NormalizeTH2D(fFullResponseOut);
4292ca60 222 // resize to desired binning scheme
549b5f40
RAB
223 TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
224 TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
4292ca60 225 // get the kinematic efficiency
549b5f40 226 TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
4292ca60 227 kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
549b5f40 228 TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
4292ca60 229 kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
230 // suppress the errors
51e6bc5a 231 for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
232 kinematicEfficiencyIn->SetBinError(1+i, 0.);
4292ca60 233 kinematicEfficiencyOut->SetBinError(1+i, 0.);
51e6bc5a 234 }
53547ff2
RAB
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 }
ad04a83c 241 // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
549b5f40
RAB
242 TH1D* unfoldedJetSpectrumIn(0x0);
243 TH1D* unfoldedJetSpectrumOut(0x0);
ad04a83c 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
486fb24e 247 // do the inplane unfolding
248 unfoldedJetSpectrumIn = UnfoldWrapper(
249 measuredJetSpectrumIn,
250 resizedResponseIn,
251 kinematicEfficiencyIn,
252 measuredJetSpectrumTrueBinsIn,
253 TString("in"),
254 jetFindingEfficiency);
549b5f40 255 resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
20bd6574 256 resizedResponseIn->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
257 resizedResponseIn->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
549b5f40
RAB
258 resizedResponseIn = ProtectHeap(resizedResponseIn);
259 resizedResponseIn->Write();
4292ca60 260 kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
261 kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
ad04a83c 262 kinematicEfficiencyIn->Write();
263 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
4292ca60 264 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
ad04a83c 265 fDetectorResponse->Write();
4292ca60 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 }
ad04a83c 277 fActiveDir->cd();
486fb24e 278 if(fDphiUnfolding) {
5e11c41c 279 TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
280 dirOut->cd();
486fb24e 281 // do the out of plane unfolding
282 unfoldedJetSpectrumOut = UnfoldWrapper(
5e11c41c 283 measuredJetSpectrumOut,
284 resizedResponseOut,
285 kinematicEfficiencyOut,
286 measuredJetSpectrumTrueBinsOut,
287 TString("out"),
288 jetFindingEfficiency);
5e11c41c 289 resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
20bd6574 290 resizedResponseOut->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
291 resizedResponseOut->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
5e11c41c 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();
ef12d5a5 311 }
5e11c41c 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");
f3ba6c8e 319 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
5e11c41c 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 }
35c03ef1 331 TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
5e11c41c 332 if(v2) {
333 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
f3ba6c8e 334 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
5e11c41c 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");
f3ba6c8e 343 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
5e11c41c 344 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
345 ratio = ProtectHeap(ratio);
346 ratio->Write();
347 }
35c03ef1 348 TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
5e11c41c 349 if(v2) {
350 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
f3ba6c8e 351 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
5e11c41c 352 v2->GetYaxis()->SetTitle("v_{2}");
353 v2 = ProtectHeap(v2);
354 v2->Write();
355 }
4292ca60 356 }
486fb24e 357 } // end of if(fDphiUnfolding)
ad04a83c 358 fDeltaPtDeltaPhi->Write();
9f892925 359 unfoldedJetSpectrumIn->Sumw2();
360 ProtectHeap(unfoldedJetSpectrumIn, kFALSE);
361 unfoldedJetSpectrumIn->Write();
362 unfoldedJetSpectrumOut->Sumw2();
363 ProtectHeap(unfoldedJetSpectrumOut, kFALSE);
364 unfoldedJetSpectrumOut->Write();
ad04a83c 365 fJetPtDeltaPhi->Write();
549b5f40
RAB
366 // save the current state of the unfolding object
367 SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
f3ba6c8e 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
dc1455ee 374}
375//_____________________________________________________________________________
486fb24e 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
87233f72 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;
fb0fb175 391 else if(fUnfoldingAlgorithm == kFold) myFunction = &AliJetFlowTools::FoldSpectrum;
87233f72 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()));
c03f7598 395 return clone;//RebinTH1D(clone, fBinsTrue, clone->GetName(), kFALSE);
486fb24e 396 }
87233f72 397 else return 0x0;
486fb24e 398 // do the actual unfolding with the selected function
399 return (this->*myFunction)( measuredJetSpectrum, resizedResponse, kinematicEfficiency, measuredJetSpectrumTrueBins, suffix, jetFindingEfficiency);
400}
401//_____________________________________________________________________________
549b5f40
RAB
402TH1D* AliJetFlowTools::UnfoldSpectrumChi2(
403 const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
486fb24e 404 const TH2D* resizedResponse, // response matrix
549b5f40
RAB
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)
dc1455ee 409{
51e6bc5a 410 // unfold the spectrum using chi2 minimization
411
4292ca60 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);
dc1455ee 421
549b5f40
RAB
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'
4292ca60 425
549b5f40
RAB
426 // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
427 TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
4292ca60 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
dc1455ee 431 // full response matrix and kinematic efficiency
549b5f40 432 TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
51e6bc5a 433 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
d7ec324f 434
4292ca60 435 // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
549b5f40
RAB
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);
ef12d5a5 439
4292ca60 440 // step 2) start the unfolding
51e6bc5a 441 Int_t status(-1), i(0);
442 while(status < 0 && i < 100) {
4292ca60 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)
549b5f40 449 measuredJetSpectrumLocal, // measured spectrum
4292ca60 450 priorLocal, // initial conditions (set NULL to use measured spectrum)
451 unfoldedLocal); // results
452 // status holds the minuit fit status (where 0 means convergence)
51e6bc5a 453 i++;
454 }
4292ca60 455 // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
1e01bfd6 456 TH2D* hPearson(0x0);
4292ca60 457 if(status == 0 && gMinuit->fISW[1] == 3) {
458 // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
ad04a83c 459 TVirtualFitter *fitter(TVirtualFitter::GetFitter());
4292ca60 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();
1e01bfd6 464 hPearson = new TH2D(*pearson);
4292ca60 465 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
466 hPearson = ProtectHeap(hPearson);
ad04a83c 467 hPearson->Write();
1e01bfd6 468 if(fMergeBinsArray) unfoldedLocal = MergeSpectrumBins(fMergeBinsArray, unfoldedLocal, hPearson);
4292ca60 469 } else status = -1;
d7ec324f 470
4292ca60 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()));
549b5f40 475 TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
4292ca60 476 if(ratio) {
53547ff2
RAB
477 ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
478 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
4292ca60 479 ratio = ProtectHeap(ratio);
480 ratio->Write();
481 }
d7ec324f 482
483 // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
484 // objects are cloned using 'ProtectHeap()'
549b5f40
RAB
485 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
486 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
487 measuredJetSpectrumLocal->Write();
d7ec324f 488
4292ca60 489 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
490 resizedResponseLocal->Write();
d7ec324f 491
4292ca60 492 unfoldedLocal = ProtectHeap(unfoldedLocal);
53547ff2 493 if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
51e6bc5a 494 unfoldedLocal->Write();
d7ec324f 495
4292ca60 496 foldedLocal = ProtectHeap(foldedLocal);
497 foldedLocal->Write();
d7ec324f 498
ef12d5a5 499 priorLocal = ProtectHeap(priorLocal);
500 priorLocal->Write();
d7ec324f 501
502 // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
53547ff2 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));
d7ec324f 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");
53547ff2
RAB
510 fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
511 fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
d7ec324f 512 fitStatus->Write();
1e01bfd6 513
549b5f40 514 return unfoldedLocal;
dc1455ee 515}
516//_____________________________________________________________________________
549b5f40
RAB
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)
51e6bc5a 524{
549b5f40
RAB
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
20abfcc4 539 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
540
4292ca60 541 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
549b5f40
RAB
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);
ef12d5a5 550 // kinematic efficiency
51e6bc5a 551 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
ef12d5a5 552 // place holder histos
d7ec324f 553 TH1D *unfoldedLocalSVD(0x0);
554 TH1D *foldedLocalSVD(0x0);
4292ca60 555 cout << " 2) setup necessary input " << endl;
51e6bc5a 556 // 3) configure routine
557 RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
d7ec324f 558 cout << " step 3) configured routine " << endl;
559
4292ca60 560 // 4) get transpose matrices
549b5f40
RAB
561 // a) get the transpose of the full response matrix
562 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
4292ca60 563 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
549b5f40
RAB
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;
4292ca60 568
569 // 5) get response for SVD unfolding
570 RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
4292ca60 571 cout << " 5) retrieved roo unfold response object " << endl;
549b5f40 572
4292ca60 573 // 6) actualy unfolding loop
549b5f40 574 RooUnfoldSvd unfoldSVD(&responseSVD, measuredJetSpectrumLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
20abfcc4 575 unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
549b5f40
RAB
576 // correct the spectrum for the kinematic efficiency
577 unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
578
579 // get the pearson coefficients from the covariance matrix
20abfcc4 580 TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
581 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
1e01bfd6 582 TH2D* hPearson(0x0);
4292ca60 583 if(pearson) {
1e01bfd6 584 hPearson = new TH2D(*pearson);
4292ca60 585 pearson->Print();
d7ec324f 586 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
587 hPearson = ProtectHeap(hPearson);
4292ca60 588 hPearson->Write();
1e01bfd6 589 if(fMergeBinsArray) unfoldedLocalSVD = MergeSpectrumBins(fMergeBinsArray, unfoldedLocalSVD, hPearson);
549b5f40 590 } else return 0x0; // return if unfolding didn't converge
4292ca60 591
592 // plot singular values and d_i vector
20abfcc4 593 TSVDUnfold* svdUnfold(unfoldSVD.Impl());
51e6bc5a 594 TH1* hSVal(svdUnfold->GetSV());
595 TH1D* hdi(svdUnfold->GetD());
4292ca60 596 hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
51e6bc5a 597 hSVal->SetXTitle("singular values");
598 hSVal->Write();
4292ca60 599 hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
51e6bc5a 600 hdi->SetXTitle("|d_{i}^{kreg}|");
601 hdi->Write();
4292ca60 602 cout << " plotted singular values and d_i vector " << endl;
603
604 // 7) refold the unfolded spectrum
549b5f40
RAB
605 foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, resizedResponseLocalNorm, kinematicEfficiencyLocal);
606 TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalSVD, "ratio measured / re-folded", kTRUE));
53547ff2 607 ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
4292ca60 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
549b5f40
RAB
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
d7ec324f 618 unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
619 unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
53547ff2 620 if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
4292ca60 621 unfoldedLocalSVD->Write(); // unfolded spectrum
d7ec324f 622 foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
623 foldedLocalSVD = ProtectHeap(foldedLocalSVD);
4292ca60 624 foldedLocalSVD->Write(); // re-folded spectrum
d7ec324f 625
549b5f40
RAB
626 // save more general bookkeeeping histograms to the output directory
627 responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
f3ba6c8e 628 responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
629 responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
4292ca60 630 responseMatrixLocalTransposePrior->Write();
549b5f40
RAB
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();
53547ff2
RAB
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
549b5f40 644 return unfoldedLocalSVD;
51e6bc5a 645}
646//_____________________________________________________________________________
549b5f40
RAB
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)
d7ec324f 654{
549b5f40
RAB
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);
4e4f12b6
RAB
666 AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
667
549b5f40
RAB
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
1e01bfd6 708 TH2D* hPearson(0x0);
549b5f40
RAB
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();
1e01bfd6 716 hPearson= new TH2D(*pearson);
549b5f40
RAB
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;
1e01bfd6 721 if(fMergeBinsArray) unfoldedLocal = MergeSpectrumBins(fMergeBinsArray, unfoldedLocal, hPearson);
549b5f40
RAB
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();
d7ec324f 733 }
549b5f40
RAB
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();
1e01bfd6 765
549b5f40
RAB
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");
d7ec324f 791 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
d7ec324f 792
793 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
549b5f40
RAB
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);
d7ec324f 802 // kinematic efficiency
803 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
804 // place holder histos
549b5f40
RAB
805 TH1D *unfoldedLocalBayes(0x0);
806 TH1D *foldedLocalBayes(0x0);
d7ec324f 807 cout << " 2) setup necessary input " << endl;
549b5f40
RAB
808 // 4) get transpose matrices
809 // a) get the transpose of the full response matrix
810 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
d7ec324f 811 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
549b5f40
RAB
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()));
d7ec324f 817
549b5f40
RAB
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);
d7ec324f 826 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
1e01bfd6 827 TH2D* hPearson(0x0);
d7ec324f 828 if(pearson) {
1e01bfd6 829 hPearson = new TH2D(*pearson);
d7ec324f 830 pearson->Print();
831 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
832 hPearson = ProtectHeap(hPearson);
833 hPearson->Write();
1e01bfd6 834 if(fMergeBinsArray) unfoldedLocalBayes = MergeSpectrumBins(fMergeBinsArray, unfoldedLocalBayes, hPearson);
549b5f40 835 } else return 0x0; // return if unfolding didn't converge
d7ec324f 836
549b5f40
RAB
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()));
d7ec324f 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
549b5f40
RAB
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");
f3ba6c8e 861 responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
862 responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
d7ec324f 863 responseMatrixLocalTransposePrior->Write();
549b5f40
RAB
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();
53547ff2
RAB
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));
549b5f40
RAB
873 fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fBayesianIterIn : fBayesianIterOut);
874 fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fBayesianIterIn" : "fBayesianIterOut");
53547ff2
RAB
875 fitStatus->Write();
876
549b5f40 877 return unfoldedLocalBayes;
d7ec324f 878}
879//_____________________________________________________________________________
e20bfd76 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)
dc1455ee 925{
926 // prepare for unfolding
4292ca60 927 if(fRawInputProvided) return kTRUE;
dc1455ee 928 if(!fInputList) {
929 printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
930 return kFALSE;
931 }
549b5f40 932 if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
4292ca60 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;
dc1455ee 937 }
486fb24e 938 if(!fRMSSpectrumIn && fDphiUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
4292ca60 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());
dc1455ee 943 }
ef12d5a5 944 if(!fTrainPower) {
549b5f40 945 // clear minuit state to avoid constraining the fit with the results of the previous iteration
ef12d5a5 946 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
947 }
c03f7598 948 // extract the spectra
949 TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
6c3fa11d 950 if(fRho0) spectrumName = Form("fHistJetPsi2PtRho0_%i", fCentralityArray->At(0));
c03f7598 951 if(!fInputList->FindObject(spectrumName.Data())) {
dc1455ee 952 printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
953 return kFALSE;
954 }
67d11165 955
c03f7598 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
4292ca60 960 fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
c03f7598 961 fJetPtDeltaPhi->Scale(fCentralityWeights->At(0));
24005d85 962 printf("Extracted %s wight weight %.2f \n", spectrumName.Data(), fCentralityWeights->At(0));
c03f7598 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
dc1455ee 970 // in plane spectrum
486fb24e 971 if(!fDphiUnfolding) {
5e11c41c 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");
ef12d5a5 974 } else {
5e11c41c 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"));
ef12d5a5 977 fSpectrumIn = ProtectHeap(fSpectrumIn);
978 // out of plane spectrum
5e11c41c 979 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30, "e");
ef12d5a5 980 fSpectrumOut = ProtectHeap(fSpectrumOut);
981 }
e20bfd76 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
20abfcc4 986 // normalize spectra to event count if requested
987 if(fNormalizeSpectra) {
c03f7598 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));
67d11165 992 }
ef12d5a5 993 if(!rho) return 0x0;
994 Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
995 if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
20abfcc4 996 if(fEventCount > 0) {
4292ca60 997 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
998 fSpectrumOut->Sumw2();
ef12d5a5 999 Double_t pt(0);
d7ec324f 1000 Double_t error(0); // lots of issues with the errors here ...
ef12d5a5 1001 for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
4292ca60 1002 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount; // normalized count
d7ec324f 1003 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
4292ca60 1004 fSpectrumIn->SetBinContent(1+i, pt);
ef12d5a5 1005 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
d7ec324f 1006 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
ef12d5a5 1007 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
4292ca60 1008 }
ef12d5a5 1009 for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
1010 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount; // normalized count
d7ec324f 1011 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
ef12d5a5 1012 fSpectrumOut->SetBinContent(1+i, pt);
1013 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
d7ec324f 1014 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
ef12d5a5 1015 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
4292ca60 1016 }
20abfcc4 1017 }
ef12d5a5 1018 if(normalizeToFullSpectrum) fEventCount = -1;
20abfcc4 1019 }
dc1455ee 1020 // extract the delta pt matrices
486fb24e 1021 TString deltaptName("");
6c3fa11d 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 }
ad04a83c 1027 fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1028 if(!fDeltaPtDeltaPhi) {
dc1455ee 1029 printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
549b5f40
RAB
1030 printf(" > may be ok, depending no what you want to do < \n");
1031 fRefreshInput = kTRUE;
1032 return kTRUE;
dc1455ee 1033 }
67d11165 1034
c03f7598 1035 // clone the distribution on the heap and if requested merge with other centralities
4292ca60 1036 fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
c03f7598 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
4292ca60 1045 // in plane delta pt distribution
486fb24e 1046 if(!fDphiUnfolding) {
5e11c41c 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");
ef12d5a5 1049 } else {
5e11c41c 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"));
ef12d5a5 1052 // out of plane delta pt distribution
5e11c41c 1053 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30, "e");
ef12d5a5 1054 fDptInDist = ProtectHeap(fDptInDist);
1055 fDptOutDist = ProtectHeap(fDptOutDist);
1056 // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
1057 }
4292ca60 1058
dc1455ee 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
ad04a83c 1062 Bool_t skip = kFALSE;
dc1455ee 1063 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
ad04a83c 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;
dc1455ee 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
ad04a83c 1070 Bool_t skip = kFALSE;
dc1455ee 1071 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
ad04a83c 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;
dc1455ee 1074 }
1075 }
1076 fDptIn = new TH2D(*rfIn);
c03f7598 1077 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
f3ba6c8e 1078 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1079 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
4292ca60 1080 fDptIn = ProtectHeap(fDptIn);
dc1455ee 1081 fDptOut = new TH2D(*rfOut);
c03f7598 1082 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
f3ba6c8e 1083 fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1084 fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
4292ca60 1085 fDptOut = ProtectHeap(fDptOut);
1086
1087 fRefreshInput = kTRUE; // force cloning of the input
dc1455ee 1088 return kTRUE;
1089}
1090//_____________________________________________________________________________
486fb24e 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 }
6c3fa11d 1108 // extract the spectra
c03f7598 1109 TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
6c3fa11d 1110 if(fRho0) spectrumName = Form("fHistJetPsi2PtRho0_%i", fCentralityArray->At(0));
486fb24e 1111 fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
1112 if(!fJetPtDeltaPhi) {
1113 printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
1114 return kFALSE;
1115 }
67d11165 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 }
486fb24e 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("");
6c3fa11d 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 }
486fb24e 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 }
67d11165 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
486fb24e 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);
c03f7598 1159 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
f3ba6c8e 1160 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1161 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
486fb24e 1162 fDptIn = ProtectHeap(fDptIn);
1163
1164 return kTRUE;
1165}
1166//_____________________________________________________________________________
549b5f40
RAB
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
d31812e1 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;
18698978 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;
549b5f40 1196 case kPriorChi2 : {
1e01bfd6 1197 TArrayI* placeHolder(0x0);
1198 if(fMergeBinsArray) {
1199 placeHolder = fMergeBinsArray;
1200 fMergeBinsArray = 0x0;
1201 }
549b5f40
RAB
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;
1e01bfd6 1207 // for the prior, do not re-bin the output
549b5f40
RAB
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");
1e01bfd6 1223 if(placeHolder) fMergeBinsArray = placeHolder;
549b5f40
RAB
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");
1e01bfd6 1240 if(placeHolder) fMergeBinsArray = placeHolder;
549b5f40
RAB
1241 return 0x0;
1242 }
1e01bfd6 1243 if(placeHolder) fMergeBinsArray = placeHolder;
549b5f40
RAB
1244 }
1245 break;
1e01bfd6 1246
549b5f40
RAB
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//_____________________________________________________________________________
dc1455ee 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//_____________________________________________________________________________
512ced40 1306TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo, Bool_t noError) {
dc1455ee 1307 // general method to normalize all vertical slices of a th2 to unity
1308 // i.e. get a probability matrix
1309 if(!histo) {
d0ef8ed2 1310 printf(" > NormalizeTH2D:: NULL pointer passed, returning NULL < \n");
dc1455ee 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);
512ced40
RAB
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);
dc1455ee 1327 }
1328 }
1329 return histo;
1330}
1331//_____________________________________________________________________________
8bc7d79b 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()))));
a98c2619 1342 bootstrapped->Reset();
1343
1344 /* OLD method - slightly underestimates fluctuations
8bc7d79b 1345 // reset the content
1346 bootstrapped->Reset();
1347 // resample the input histogram
a98c2619 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
8bc7d79b 1370 // if requested kill input histo
1371 if(kill) delete hist;
1372 // return resampled histogram
1373 return bootstrapped;
1374}
1375//_____________________________________________________________________________
53547ff2 1376TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
dc1455ee 1377 // return a TH1D with the supplied histogram rebinned to the supplied bins
53547ff2 1378 // the returned histogram is new, the original is deleted from the heap if kill is true
dc1455ee 1379 if(!histo || !bins) {
53547ff2 1380 printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
dc1455ee 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++) {
4292ca60 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));
dc1455ee 1391 }
53547ff2 1392 if(kill) delete histo;
dc1455ee 1393 return rebinned;
1394}
1395//_____________________________________________________________________________
4292ca60 1396TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
53547ff2 1397 // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
8bc7d79b 1398 // not static as it is just a wrapper for the response maker object
4292ca60 1399 if(!fResponseMaker || !binsTrue || !binsRec) {
1400 printf(" > RebinTH2D:: function called with NULL arguments < \n");
1401 return 0x0;
dc1455ee 1402 }
4292ca60 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);
dc1455ee 1405}
1406//_____________________________________________________________________________
d7ec324f 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);
512ced40 1420 c->SetBinError(x2, y1, 0.);
dc1455ee 1421 }
1422 }
d7ec324f 1423 if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1424 return c;
dc1455ee 1425}
1426//_____________________________________________________________________________
d7ec324f 1427TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale)
1428{
dc1455ee 1429 // normalize a th1d to a certain scale
4292ca60 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;
dc1455ee 1436}
1437//_____________________________________________________________________________
1e01bfd6 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//_____________________________________________________________________________
51e6bc5a 1467TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
dc1455ee 1468{
1469 // Calculate pearson coefficients from covariance matrix
51e6bc5a 1470 TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1471 Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
ad04a83c 1472 Double_t pearson(0.);
dc1455ee 1473 if(nrows==0 && ncols==0) return 0x0;
ef12d5a5 1474 for(Int_t row = 0; row < nrows; row++) {
1475 for(Int_t col = 0; col<ncols; col++) {
51e6bc5a 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;
dc1455ee 1478 }
1479 }
51e6bc5a 1480 return pearsonCoefficients;
dc1455ee 1481}
1482//_____________________________________________________________________________
549b5f40 1483TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
d7ec324f 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)
549b5f40
RAB
1489 if(!spectrum || !function) return 0x0;
1490 if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
d7ec324f 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
87233f72 1493 TFitResultPtr r = temp->Fit(function, "", "", min, max);
d7ec324f 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
486fb24e 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));
d7ec324f 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//_____________________________________________________________________________
20bd6574 1506void AliJetFlowTools::Style(Bool_t legacy)
18698978 1507{
1508 // set global style for your current aliroot session
1509 if(!gStyle) return;
20bd6574 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 }
18698978 1559}
1560//_____________________________________________________________________________
d7ec324f 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//_____________________________________________________________________________
20bd6574 1582void AliJetFlowTools::Style(TVirtualPad* c, TString style, Bool_t legacy)
d7ec324f 1583{
20bd6574 1584 // set a default style for a canva
1585
1586 if(legacy) {
1587 c->SetLeftMargin(.25);
1588 c->SetBottomMargin(.25);
1589 }
1590 else Style();
d7ec324f 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;
18698978 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();
d7ec324f 1612 } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1613}
1614//_____________________________________________________________________________
53547ff2
RAB
1615void AliJetFlowTools::Style(TLegend* l)
1616{
1617 // set a default style for a legend
53547ff2 1618 l->SetFillColor(0);
18698978 1619 l->SetBorderSize(0);
20bd6574 1620 if(gStyle) l->SetTextSize(gStyle->GetTextSize()*.08);
53547ff2
RAB
1621}
1622//_____________________________________________________________________________
20bd6574 1623void AliJetFlowTools::Style(TH1* h, EColor col, histoType type, Bool_t legacy)
53547ff2
RAB
1624{
1625 // style a histo
e4b6d21a 1626 h->GetYaxis()->SetNdivisions(505);
53547ff2
RAB
1627 h->SetLineColor(col);
1628 h->SetMarkerColor(col);
20bd6574 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();
53547ff2
RAB
1640 switch (type) {
1641 case kInPlaneSpectrum : {
1642 h->SetTitle("IN PLANE");
20bd6574 1643 h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
18698978 1644 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
53547ff2
RAB
1645 } break;
1646 case kOutPlaneSpectrum : {
1647 h->SetTitle("OUT OF PLANE");
18698978 1648 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
e4b6d21a 1649 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
53547ff2
RAB
1650 } break;
1651 case kUnfoldedSpectrum : {
1652 h->SetTitle("UNFOLDED");
18698978 1653 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
e4b6d21a 1654 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
53547ff2
RAB
1655 } break;
1656 case kFoldedSpectrum : {
1657 h->SetTitle("FOLDED");
18698978 1658 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
e4b6d21a 1659 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
53547ff2
RAB
1660 } break;
1661 case kMeasuredSpectrum : {
1662 h->SetTitle("MEASURED");
18698978 1663 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
e4b6d21a 1664 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{cht} (GeV/#it{c})");
53547ff2 1665 } break;
18698978 1666 case kBar : {
1667 h->SetFillColor(col);
1668 h->SetBarWidth(.6);
e4b6d21a 1669 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
18698978 1670 h->SetBarOffset(0.2);
1671 }
1e01bfd6 1672 case kRatio : {
1673 h->SetMarkerStyle(8);
1674 h->SetMarkerSize(1);
1675 } break;
6c3fa11d 1676 case kDeltaPhi : {
1677 h->GetYaxis()->SetTitle("[counts]");
1678 h->GetXaxis()->SetTitle("#Delta #phi");
1679 }
53547ff2
RAB
1680 default : break;
1681 }
1682}
1683//_____________________________________________________________________________
20bd6574 1684void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type, Bool_t legacy)
53547ff2 1685{
20bd6574 1686 // style a tgraph
e4b6d21a 1687 h->GetYaxis()->SetNdivisions(505);
53547ff2
RAB
1688 h->SetLineColor(col);
1689 h->SetMarkerColor(col);
20bd6574 1690 h->SetLineWidth(2);
1691 h->SetMarkerSize(1);
18698978 1692 h->SetTitle("");
d06dbffe 1693 h->SetFillColor(kCyan);
20bd6574 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();
e4b6d21a 1702 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
53547ff2
RAB
1703 switch (type) {
1704 case kInPlaneSpectrum : {
1705 h->SetTitle("IN PLANE");
18698978 1706 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
53547ff2
RAB
1707 } break;
1708 case kOutPlaneSpectrum : {
1709 h->SetTitle("OUT OF PLANE");
18698978 1710 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
53547ff2
RAB
1711 } break;
1712 case kUnfoldedSpectrum : {
1713 h->SetTitle("UNFOLDED");
18698978 1714 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
53547ff2
RAB
1715 } break;
1716 case kFoldedSpectrum : {
1717 h->SetTitle("FOLDED");
18698978 1718 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
53547ff2
RAB
1719 } break;
1720 case kMeasuredSpectrum : {
1721 h->SetTitle("MEASURED");
18698978 1722 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
53547ff2 1723 } break;
a39e4b2b 1724 case kRatio : {
20bd6574 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))");
a39e4b2b 1727 } break;
1728 case kV2 : {
20bd6574 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}}}");
3fbc711e 1730 h->GetYaxis()->SetTitle("#it{v}_{2}^{ch, jet} \{EP, |#Delta#eta|>0.9 \} ");
a39e4b2b 1731 h->GetYaxis()->SetRangeUser(-.5, 1.);
1e01bfd6 1732 h->SetMarkerStyle(8);
1733 h->SetMarkerSize(1);
a39e4b2b 1734 } break;
53547ff2
RAB
1735 default : break;
1736 }
1737}
1738//_____________________________________________________________________________
a39e4b2b 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
a39e4b2b 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),
24005d85 1775 fBinsTrue->At(fBinsTrue->GetSize()-1),
a39e4b2b 1776 readMe,
1777 "nominal_values");
35c03ef1 1778 v2 = GetV2(nominalIn, nominalOut, fEventPlaneRes, "nominal v_{2}");
a39e4b2b 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();
a39e4b2b 1783}
1784//_____________________________________________________________________________
18698978 1785void AliJetFlowTools::GetCorrelatedUncertainty(
d31812e1 1786 TGraphAsymmErrors*& corrRatio, // correlated uncertainty pointer reference
1787 TGraphAsymmErrors*& corrV2, // correlated uncertainty pointer reference
1788 TArrayI* variationsIn, // variantions in plane
a39e4b2b 1789 TArrayI* variationsOut, // variantions out of plane
7b88bd32 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
18698978 1794 TString type, // type of uncertaitny
24005d85 1795 TString type2,
18698978 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
5159ed9e 1799 Float_t corr, // correlation strength
18698978 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
f3ba6c8e 1803{
18698978 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);
98d5b614 1820 TH1D* relativeError2ndVariationInLow(0x0);
1821 TH1D* relativeError2ndVariationInUp(0x0);
1822 TH1D* relativeError2ndVariationOutLow(0x0);
1823 TH1D* relativeError2ndVariationOutUp(0x0);
18698978 1824 TH1D* relativeStatisticalErrorIn(0x0);
1825 TH1D* relativeStatisticalErrorOut(0x0);
a39e4b2b 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()));
18698978 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
a39e4b2b 1840 relativeStatisticalErrorIn, // pointer reference
1841 relativeStatisticalErrorOut, // pointer reference
1842 nominal,
1843 nominalIn,
1844 nominalOut,
18698978 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();
14ce5e08 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 }
18698978 1891 }
1892 }
98d5b614 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,
24005d85 1911 type2);
98d5b614 1912 if(relativeError2ndVariationInUp) {
1913 // canvas with the error from variation strength
24005d85 1914 TCanvas* relativeError2ndVariation(new TCanvas(Form("relativeError2nd_%s", type2.Data()), Form("relativeError2nd_%s", type2.Data())));
98d5b614 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 }
18698978 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
a39e4b2b 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);
98d5b614 1952 if(relativeError2ndVariationInUp) bInUp = relativeError2ndVariationInUp->GetBinContent(b+1);
1953 if(relativeError2ndVariationOutUp) bInLow = relativeError2ndVariationOutUp->GetBinContent(b+1);
a39e4b2b 1954 dInUp = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp;
7b88bd32 1955 // for a symmetric first variation
1956 if(sym) dInUp += aInLow*aInLow;
a39e4b2b 1957 if(dInUp > 0) relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(dInUp));
1958 dOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp;
d06dbffe 1959 if(sym) dOutUp += aOutLow*aOutLow;
a39e4b2b 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);
98d5b614 1964 if(relativeError2ndVariationInLow) bInLow = relativeError2ndVariationInLow->GetBinContent(b+1);
1965 if(relativeError2ndVariationOutLow) bOutLow = relativeError2ndVariationOutLow->GetBinContent(b+1);
a39e4b2b 1966 dInLow = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow;
7b88bd32 1967 if(sym) dInLow += aInUp*aInUp;
5159ed9e 1968 if(dInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1*TMath::Sqrt(dInLow));
a39e4b2b 1969 dOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow;
d06dbffe 1970 if(sym) dOutLow += aOutUp*aOutUp;
a39e4b2b 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
b3b03df7 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 }
5159ed9e 1995 ayl[i] = TMath::Sqrt(TMath::Abs(lowErr))*nominal->GetBinContent(i+1);
1996 ayh[i] = TMath::Sqrt(TMath::Abs(upErr))*nominal->GetBinContent(i+1);
a39e4b2b 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);
d06dbffe 2013 Style(nominalError, kCyan, kRatio);
2014 nominalError->SetLineColor(kCyan);
2015 nominalError->SetMarkerColor(kCyan);
a39e4b2b 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,
35c03ef1 2027 fEventPlaneRes,
5159ed9e 2028 "v_{2} with shape uncertainty",
a39e4b2b 2029 relativeErrorInUp,
2030 relativeErrorInLow,
24005d85 2031 relativeErrorOutUp,
5159ed9e 2032 relativeErrorOutLow,
2033 corr));
a39e4b2b 2034 // pass the nominal values to the pointer references
2035 corrV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
35c03ef1 2036 TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
a39e4b2b 2037 nominalCanvas->cd(2);
2038 Style(nominalV2, kBlack);
d06dbffe 2039 Style(nominalV2Error, kCyan, kV2);
2040 nominalV2Error->SetLineColor(kCyan);
2041 nominalV2Error->SetMarkerColor(kCyan);
a39e4b2b 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
18698978 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();
a39e4b2b 2079// output->Close();
18698978 2080}
2081//_____________________________________________________________________________
2082void AliJetFlowTools::GetShapeUncertainty(
a39e4b2b 2083 TGraphAsymmErrors*& shapeRatio, // pointer reference to final shape uncertainty of ratio
2084 TGraphAsymmErrors*& shapeV2, // pointer reference to final shape uncertainty of v2
18698978 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
b2150106 2091 TArrayI* methodIn, // method in
2092 TArrayI* methodOut, // method out
18698978 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
b9bd4252 2096 Float_t corr, // correlation strength
18698978 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);
b2150106 2119 TH1D* relativeErrorMethodInLow(0x0);
2120 TH1D* relativeErrorMethodInUp(0x0);
18698978 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);
b2150106 2129 TH1D* relativeErrorMethodOutLow(0x0);
2130 TH1D* relativeErrorMethodOutUp(0x0);
a39e4b2b 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()));
18698978 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
a39e4b2b 2147 nominal,
2148 nominalIn,
2149 nominalOut,
18698978 2150 columns,
2151 rangeLow,
2152 rangeUp,
2153 readMe,
d06dbffe 2154 "regularization",
3e698d27 2155 fRMS);
18698978 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,
a39e4b2b 2184 nominal,
2185 nominalIn,
2186 nominalOut,
18698978 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,
24005d85 2213 relativeErrorRecBinInUp, // pointer reference
2214 relativeErrorRecBinInLow, // pointer reference
2215 relativeErrorRecBinOutUp, // pointer reference
2216 relativeErrorRecBinOutLow, // pointer reference
18698978 2217 relativeStatisticalErrorIn,
2218 relativeStatisticalErrorOut,
a39e4b2b 2219 nominal,
2220 nominalIn,
2221 nominalOut,
18698978 2222 columns,
2223 rangeLow,
2224 rangeUp,
2225 readMe,
14ce5e08 2226 "recBin",
2227 fRMS);
18698978 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 }
b2150106 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,
db02017e 2263 "method"
2264 );
b2150106 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 }
18698978 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
b2150106 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.);
18698978 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);
b2150106 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));
18698978 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);
b2150106 2320 if(relativeErrorMethodInLow) dInLow = relativeErrorMethodInLow->GetBinContent(b+1);
2321 if(relativeErrorMethodOutLow) dOutLow = relativeErrorMethodOutLow->GetBinContent(b+1);
db02017e 2322 if(fSymmRMS) { // take first category as symmetric
d0ef8ed2 2323 aInLow = aInUp;
2324 aOutLow = aOutUp;
14ce5e08 2325 cInLow = cInUp;
2326 cOutLow = cOutUp; // other sources
db02017e 2327 if(dInLow < dInUp) dInLow = dInUp;
2328 if(dOutLow < dOutUp) dOutLow = dOutUp;
2329 }
2330
b2150106 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));
18698978 2335 }
a39e4b2b 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
24005d85 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);
a39e4b2b 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);
d06dbffe 2370 Style(nominalError, kCyan, kRatio);
2371 nominalError->SetLineColor(kCyan);
2372 nominalError->SetMarkerColor(kCyan);
a39e4b2b 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,
35c03ef1 2384 fEventPlaneRes,
a39e4b2b 2385 "v_{2} with shape uncertainty",
2386 relativeErrorInUp,
2387 relativeErrorInLow,
2388 relativeErrorOutUp,
b9bd4252 2389 relativeErrorOutLow,
2390 corr));
a39e4b2b 2391 shapeV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
35c03ef1 2392 TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
a39e4b2b 2393 nominalCanvas->cd(2);
2394 Style(nominalV2, kBlack);
d06dbffe 2395 Style(nominalV2Error, kCyan, kV2);
2396 nominalV2Error->SetLineColor(kCyan);
2397 nominalV2Error->SetMarkerColor(kCyan);
a39e4b2b 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
18698978 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();
a39e4b2b 2435// output->Close();
18698978 2436}
2437//_____________________________________________________________________________
a39e4b2b 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)
d06dbffe 2454 TString source, // source of the variation
2455 Bool_t RMS // return RMS of distribution of variations as error
a39e4b2b 2456 ) const
18698978 2457{
2458 // intermediate systematic check function. first index of supplied array is nominal value
18698978 2459 TList* listOfKeys((TList*)readMe->GetListOfKeys());
f3ba6c8e 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 }
18698978 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())));
f3ba6c8e 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());
18698978 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
d06dbffe 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
1e01bfd6 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());
18698978 2495 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
3e698d27 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 }
18698978 2515 }
3e698d27 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
18698978 2521 // define an output histogram with the systematic error from this systematic constribution
2522 if(!relativeStatisticalErrorIn && !relativeStatisticalErrorOut) {
1e01bfd6 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());
18698978 2525 }
2526
f3ba6c8e 2527 // prepare necessary canvasses
18698978 2528 TCanvas* canvasIn(new TCanvas(Form("SYST_%s_PearsonIn", source.Data()), Form("SYST_%s_PearsonIn", source.Data())));
f3ba6c8e 2529 TCanvas* canvasOut(0x0);
18698978 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())));
f3ba6c8e 2532 TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
18698978 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())));
f3ba6c8e 2535 TCanvas* canvasSpectraOut(0x0);
18698978 2536 if(fDphiUnfolding) canvasSpectraOut = new TCanvas(Form("SYST_%s_SpectraOut", source.Data()), Form("SYST_%s_SpectraOut", source.Data()));
f3ba6c8e 2537 TCanvas* canvasRatio(0x0);
18698978 2538 if(fDphiUnfolding) canvasRatio = new TCanvas(Form("SYST_%s_Ratio", source.Data()), Form("SYST_%s_Ratio", source.Data()));
f3ba6c8e 2539 TCanvas* canvasV2(0x0);
18698978 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())));
f3ba6c8e 2543 TCanvas* canvasMasterOut(0x0);
18698978 2544 if(fDphiUnfolding) canvasMasterOut = new TCanvas(Form("SYST_%s_defaultOut", source.Data()), Form("SYST_%s_defaultOut", source.Data()));
f3ba6c8e 2545 (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
a39e4b2b 2546
18698978 2547 TCanvas* canvasProfiles(new TCanvas(Form("SYST_%s_canvasProfiles", source.Data()), Form("SYST_%s_canvasProfiles", source.Data())));
f3ba6c8e 2548 canvasProfiles->Divide(2);
18698978 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()));
f3ba6c8e 2551 // get an estimate of the number of outputs and find the default set
5c636e0d 2552
2553 Int_t rows = 1;
2554 columns = variationsIn->GetSize()-1;
2555 (TMath::Floor(variationsIn->GetSize()/(float)columns)+((variationsIn->GetSize()%columns)>0));
f3ba6c8e 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
18698978 2568 TCanvas* canvasNominalIn(new TCanvas(Form("Nominal_%s_PearsonIn", source.Data()), Form("Nominal_%s_PearsonIn", source.Data())));
f3ba6c8e 2569 TCanvas* canvasNominalOut(0x0);
18698978 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())));
f3ba6c8e 2572 TCanvas* canvasNominalRatioMeasuredRefoldedOut(0x0);
18698978 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())));
f3ba6c8e 2575 TCanvas* canvasNominalSpectraOut(0x0);
18698978 2576 if(fDphiUnfolding) canvasNominalSpectraOut = new TCanvas(Form("Nominal_%s_SpectraOut", source.Data()), Form("Nominal_%s_SpectraOut", source.Data()));
f3ba6c8e 2577 TCanvas* canvasNominalRatio(0x0);
18698978 2578 if(fDphiUnfolding) canvasNominalRatio = new TCanvas(Form("Nominal_%s_Ratio", source.Data()), Form("Nominal_%s_Ratio", source.Data()));
f3ba6c8e 2579 TCanvas* canvasNominalV2(0x0);
18698978 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())));
f3ba6c8e 2583 TCanvas* canvasNominalMasterOut(0x0);
18698978 2584 if(fDphiUnfolding) canvasNominalMasterOut = new TCanvas(Form("Nominal_%s_defaultOut", source.Data()), Form("Nominal_%s_defaultOut", source.Data()));
f3ba6c8e 2585 (fDphiUnfolding) ? canvasNominalMISC->Divide(4, 2) : canvasNominalMISC->Divide(4, 1);
18698978 2586
2587 canvasNominalSpectraIn->Divide(2);
2588 if(canvasNominalSpectraOut) canvasNominalSpectraOut->Divide(2);
2589
2590 canvasNominalMasterIn->Divide(2);
2591 if(canvasNominalMasterOut) canvasNominalMasterOut->Divide(2);
f3ba6c8e 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");
a39e4b2b 2601
f3ba6c8e 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++) {
18698978 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())));
f3ba6c8e 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);
18698978 2637 Style(gPad, "GRID");
f3ba6c8e 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);
2835b296 2664 Style(gPad, "GRID");
f3ba6c8e 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);
2835b296 2675 Style(gPad, "GRID");
f3ba6c8e 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);
18698978 2688 // get the absolute relative error
2689 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
d06dbffe 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 }
3e698d27 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) {
d06dbffe 2704 relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
3e698d27 2705 relativeErrorInUpN[b]++;
d06dbffe 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));
3e698d27 2710 relativeErrorInLowN[b]++;
d06dbffe 2711 }
3e698d27 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]++;
18698978 2716 }
2717 if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorIn->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2718 }
f3ba6c8e 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);
f3ba6c8e 2723 temp->GetYaxis()->SetRangeUser(0., 2);
18698978 2724 Style(gPad, "GRID");
f3ba6c8e 2725 temp->DrawCopy();
18698978 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 }
f3ba6c8e 2735 }
2736 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirNameIn.Data())));
2737 canvasSpectraIn->cd(j);
18698978 2738 if(i==0) canvasNominalSpectraIn->cd(1);
f3ba6c8e 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 }
18698978 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");
f3ba6c8e 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);
18698978 2781 Style(gPad, "GRID");
f3ba6c8e 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);
2835b296 2808 Style(gPad, "GRID");
f3ba6c8e 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);
2835b296 2819 Style(gPad, "GRID");
f3ba6c8e 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);
18698978 2832 // get the absolute relative error
2833 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
d06dbffe 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 }
3e698d27 2845 } else if (RMS && !fSymmRMS) {
2846 printf(" OOps \n ");
d06dbffe 2847 if(temp->GetBinContent(b+1) < 1) {
2848 relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
3e698d27 2849 relativeErrorOutUpN[b]++;
d06dbffe 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));
3e698d27 2853 relativeErrorOutLowN[b]++;
d06dbffe 2854 }
3e698d27 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]++;
18698978 2859 }
2860 if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorOut->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
d06dbffe 2861 }
18698978 2862 temp->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
f3ba6c8e 2863 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2864 temp->GetYaxis()->SetTitle("ratio");
2865 canvasMasterOut->cd(j);
18698978 2866 temp->GetYaxis()->SetRangeUser(0., 2);
2867 Style(gPad, "GRID");
f3ba6c8e 2868 temp->DrawCopy();
18698978 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 }
f3ba6c8e 2878 }
2879 TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirNameOut.Data())));
2880 canvasSpectraOut->cd(j);
18698978 2881 if(i==0) canvasNominalSpectraOut->cd(1);
f3ba6c8e 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));