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