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