]> 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
Merge branch 'master' into dev
[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");
1492 v2 = GetV2(nominalIn, nominalOut, .5666666, "nominal v_{2}");
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,
1707 .56,
1708 "v_{2} with correlated uncertainty",
1709 relativeErrorInUp,
1710 relativeErrorInLow,
1711 relativeErrorOutUp,
1712 relativeErrorOutLow));
1713 // pass the nominal values to the pointer references
1714 corrV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
1715 TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, .56, "v_{2}"));
1716 nominalCanvas->cd(2);
1717 Style(nominalV2, kBlack);
1718 Style(nominalV2Error, kYellow, kV2);
1719 nominalV2Error->SetLineColor(kYellow);
1720 nominalV2Error->SetMarkerColor(kYellow);
1721 nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1722 nominalV2Error->DrawClone("a2");
1723 nominalV2->DrawClone("same E1");
1724 Style(AddLegend(gPad));
1725 Style(gPad, "GRID");
1726 Style(nominalCanvas);
1727 SavePadToPDF(nominalCanvas);
1728 nominalCanvas->Write();
1729 }
1730
18698978 1731 TCanvas* relativeError(new TCanvas("relativeCorrelatedError"," relativeCorrelatedError"));
1732 relativeError->Divide(2);
1733 relativeError->cd(1);
1734 Style(gPad, "GRID");
1735 relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
1736 Style(relativeErrorInUp, kBlue, kBar);
1737 Style(relativeErrorInLow, kGreen, kBar);
1738 relativeErrorInUp->DrawCopy("b");
1739 relativeErrorInLow->DrawCopy("same b");
1740 Style(relativeStatisticalErrorIn, kRed);
1741 relativeStatisticalErrorIn->DrawCopy("same");
1742 Style(AddLegend(gPad));
1743 relativeError->cd(2);
1744 Style(gPad, "GRID");
1745 relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
1746 Style(relativeErrorOutUp, kBlue, kBar);
1747 Style(relativeErrorOutLow, kGreen, kBar);
1748 relativeErrorOutUp->DrawCopy("b");
1749 relativeErrorOutLow->DrawCopy("same b");
1750 Style(relativeStatisticalErrorOut, kRed);
1751 relativeStatisticalErrorOut->DrawCopy("same");
1752 Style(AddLegend(gPad));
1753
1754 // write the buffered file to disk and close the file
1755 SavePadToPDF(relativeError);
1756 relativeError->Write();
1757 output->Write();
a39e4b2b 1758// output->Close();
18698978 1759}
1760//_____________________________________________________________________________
1761void AliJetFlowTools::GetShapeUncertainty(
a39e4b2b 1762 TGraphAsymmErrors*& shapeRatio, // pointer reference to final shape uncertainty of ratio
1763 TGraphAsymmErrors*& shapeV2, // pointer reference to final shape uncertainty of v2
18698978 1764 TArrayI* regularizationIn, // regularization strength systematics, in plane
1765 TArrayI* regularizationOut, // regularization strenght systematics, out of plane
1766 TArrayI* trueBinIn, // true bin ranges, in plane
1767 TArrayI* trueBinOut, // true bin ranges, out of plane
1768 TArrayI* recBinIn, // rec bin ranges, in plane
1769 TArrayI* recBinOut, // rec bin ranges, out of plane
1770 Int_t columns, // divide the output canvasses in this many columns
1771 Float_t rangeLow, // lower pt range
1772 Float_t rangeUp, // upper pt range
1773 TString in, // input file name (created by this unfolding class)
1774 TString out // output file name (which will hold results of the systematic test)
1775 ) const
1776{
1777 // do full systematics
1778 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close(); // if for some weird reason the unfolding output is still mutable
1779 TFile* readMe(new TFile(in.Data(), "READ")); // open unfolding output read-only
1780 if(readMe->IsZombie()) {
1781 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1782 return;
1783 }
1784 printf("\n\n\n\t\t DOSYSTEMATICS \n > Recovered the following file structure : \n <");
1785 readMe->ls();
1786 TFile* output(new TFile(out.Data(), "RECREATE")); // create new output file
1787
1788 // create some null placeholder pointers
1789 TH1D* relativeErrorRegularizationInLow(0x0);
1790 TH1D* relativeErrorRegularizationInUp(0x0);
1791 TH1D* relativeErrorTrueBinInLow(0x0);
1792 TH1D* relativeErrorTrueBinInUp(0x0);
1793 TH1D* relativeErrorRecBinInLow(0x0);
1794 TH1D* relativeErrorRecBinInUp(0x0);
1795 TH1D* relativeErrorRegularizationOutLow(0x0);
1796 TH1D* relativeErrorRegularizationOutUp(0x0);
1797 TH1D* relativeErrorTrueBinOutLow(0x0);
1798 TH1D* relativeErrorTrueBinOutUp(0x0);
1799 TH1D* relativeErrorRecBinOutLow(0x0);
1800 TH1D* relativeErrorRecBinOutUp(0x0);
1801 TH1D* relativeStatisticalErrorIn(0x0);
1802 TH1D* relativeStatisticalErrorOut(0x0);
a39e4b2b 1803 // histo for the nominal ratio in / out
1804 TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1805 TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1806 TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
18698978 1807
1808 // call the functions
1809 if(regularizationIn && regularizationOut) {
1810 DoIntermediateSystematics(
1811 regularizationIn,
1812 regularizationOut,
1813 relativeErrorRegularizationInUp, // pointer reference
1814 relativeErrorRegularizationInLow, // pointer reference
1815 relativeErrorRegularizationOutUp, // pointer reference
1816 relativeErrorRegularizationOutLow, // pointer reference
1817 relativeStatisticalErrorIn, // pointer reference
1818 relativeStatisticalErrorOut, // pointer reference
a39e4b2b 1819 nominal,
1820 nominalIn,
1821 nominalOut,
18698978 1822 columns,
1823 rangeLow,
1824 rangeUp,
1825 readMe,
1826 "regularization");
1827 if(relativeErrorRegularizationInUp) {
1828 // canvas with the error from regularization strength
1829 TCanvas* relativeErrorRegularization(new TCanvas("relativeErrorRegularization", "relativeErrorRegularization"));
1830 relativeErrorRegularization->Divide(2);
1831 relativeErrorRegularization->cd(1);
1832 Style(gPad, "GRID");
1833 relativeErrorRegularizationInUp->DrawCopy("b");
1834 relativeErrorRegularizationInLow->DrawCopy("same b");
1835 Style(AddLegend(gPad));
1836 relativeErrorRegularization->cd(2);
1837 Style(gPad, "GRID");
1838 relativeErrorRegularizationOutUp->DrawCopy("b");
1839 relativeErrorRegularizationOutLow->DrawCopy("same b");
1840 SavePadToPDF(relativeErrorRegularization);
1841 Style(AddLegend(gPad));
1842 relativeErrorRegularization->Write();
1843 }
1844 }
1845 if(trueBinIn && trueBinOut) {
1846 DoIntermediateSystematics(
1847 trueBinIn,
1848 trueBinOut,
1849 relativeErrorTrueBinInUp, // pointer reference
1850 relativeErrorTrueBinInLow, // pointer reference
1851 relativeErrorTrueBinOutUp, // pointer reference
1852 relativeErrorTrueBinOutLow, // pointer reference
1853 relativeStatisticalErrorIn,
1854 relativeStatisticalErrorOut,
a39e4b2b 1855 nominal,
1856 nominalIn,
1857 nominalOut,
18698978 1858 columns,
1859 rangeLow,
1860 rangeUp,
1861 readMe,
1862 "trueBin");
1863 if(relativeErrorTrueBinInUp) {
1864 TCanvas* relativeErrorTrueBin(new TCanvas("relativeErrorTrueBin", "relativeErrorTrueBin"));
1865 relativeErrorTrueBin->Divide(2);
1866 relativeErrorTrueBin->cd(1);
1867 Style(gPad, "GRID");
1868 relativeErrorTrueBinInUp->DrawCopy("b");
1869 relativeErrorTrueBinInLow->DrawCopy("same b");
1870 Style(AddLegend(gPad));
1871 relativeErrorTrueBin->cd(2);
1872 Style(gPad, "GRID");
1873 relativeErrorTrueBinOutUp->DrawCopy("b");
1874 relativeErrorTrueBinOutLow->DrawCopy("same b");
1875 SavePadToPDF(relativeErrorTrueBin);
1876 Style(AddLegend(gPad));
1877 relativeErrorTrueBin->Write();
1878 }
1879 }
1880 if(recBinIn && recBinOut) {
1881 DoIntermediateSystematics(
1882 recBinIn,
1883 recBinOut,
1884 relativeErrorRecBinInLow, // pointer reference
1885 relativeErrorRecBinInUp, // pointer reference
1886 relativeErrorRecBinOutLow, // pointer reference
1887 relativeErrorRecBinOutUp, // pointer reference
1888 relativeStatisticalErrorIn,
1889 relativeStatisticalErrorOut,
a39e4b2b 1890 nominal,
1891 nominalIn,
1892 nominalOut,
18698978 1893 columns,
1894 rangeLow,
1895 rangeUp,
1896 readMe,
1897 "recBin");
1898 if(relativeErrorRecBinOutUp) {
1899 // canvas with the error from regularization strength
1900 TCanvas* relativeErrorRecBin(new TCanvas("relativeErrorRecBin"," relativeErrorRecBin"));
1901 relativeErrorRecBin->Divide(2);
1902 relativeErrorRecBin->cd(1);
1903 Style(gPad, "GRID");
1904 relativeErrorRecBinInUp->DrawCopy("b");
1905 relativeErrorRecBinInLow->DrawCopy("same b");
1906 Style(AddLegend(gPad));
1907 relativeErrorRecBin->cd(2);
1908 Style(gPad, "GRID");
1909 relativeErrorRecBinOutUp->DrawCopy("b");
1910 relativeErrorRecBinOutLow->DrawCopy("same b");
1911 SavePadToPDF(relativeErrorRecBin);
1912 Style(AddLegend(gPad));
1913 relativeErrorRecBin->Write();
1914 }
1915 }
1916
1917 // and the placeholder for the final systematic
1918 TH1D* relativeErrorInUp(new TH1D("max shape uncertainty in plane", "max shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1919 TH1D* relativeErrorOutUp(new TH1D("max shape uncertainty out of plane", "max shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1920 TH1D* relativeErrorInLow(new TH1D("min shape uncertainty in plane", "min shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1921 TH1D* relativeErrorOutLow(new TH1D("min shape uncertainty out of plane", "min shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1922 relativeErrorInUp->SetYTitle("relative uncertainty");
1923 relativeErrorOutUp->SetYTitle("relative uncertainty");
1924 relativeErrorInLow->SetYTitle("relative uncertainty");
1925 relativeErrorOutLow->SetYTitle("relative uncertainty");
1926
1927 // sum of squares for the total systematic error
1928 Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.);
1929 Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.);
1930 Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.);
1931 Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.);
1932
1933 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
1934 // for the upper bound
1935 if(relativeErrorRegularizationInUp) aInUp = relativeErrorRegularizationInUp->GetBinContent(b+1);
1936 if(relativeErrorRegularizationOutUp) aOutUp = relativeErrorRegularizationOutUp->GetBinContent(b+1);
1937 if(relativeErrorTrueBinInUp) bInUp = relativeErrorTrueBinInUp->GetBinContent(b+1);
1938 if(relativeErrorTrueBinOutUp) bOutUp = relativeErrorTrueBinOutUp->GetBinContent(b+1);
1939 if(relativeErrorRecBinInUp) cInUp = relativeErrorRecBinInUp->GetBinContent(b+1);
1940 if(relativeErrorRecBinOutUp) cOutUp = relativeErrorRecBinOutUp->GetBinContent(b+1);
1941 dInUp = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp;
1942 if(dInUp > 0) relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(dInUp));
1943 dOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp;
1944 if(dOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(dOutUp));
1945 // for the lower bound
1946 if(relativeErrorRegularizationInLow) aInLow = relativeErrorRegularizationInLow->GetBinContent(b+1);
1947 if(relativeErrorRegularizationOutLow) aOutLow = relativeErrorRegularizationOutLow->GetBinContent(b+1);
1948 if(relativeErrorTrueBinInLow) bInLow = relativeErrorTrueBinInLow->GetBinContent(b+1);
1949 if(relativeErrorTrueBinOutLow) bOutLow = relativeErrorTrueBinOutLow->GetBinContent(b+1);
1950 if(relativeErrorRecBinInLow) cInLow = relativeErrorRecBinInLow->GetBinContent(b+1);
1951 if(relativeErrorRecBinOutLow) cOutLow = relativeErrorRecBinOutLow->GetBinContent(b+1);
1952 dInLow = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow;
1953 if(dInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(dInLow));
1954 dOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow;
1955 if(dOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(dOutLow));
1956 }
a39e4b2b 1957 // project the estimated errors on the nominal ratio
1958 if(nominal) {
1959 Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
1960 Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
1961 Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
1962 Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
1963 Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
1964 Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
1965 Double_t lowErr(0.), upErr(0.);
1966 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
1967 // add the in and out of plane errors in quadrature
1968 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(1+i)*relativeErrorOutLow->GetBinContent(i+1);
1969 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1);
1970 // set the errors
1971 ayl[i] = TMath::Sqrt(lowErr)*nominal->GetBinContent(i+1);
1972 ayh[i] = TMath::Sqrt(upErr)*nominal->GetBinContent(i+1);
1973 // get the bin width (which is the 'error' on x
1974 Double_t binWidth(nominal->GetBinWidth(i+1));
1975 axl[i] = binWidth/2.;
1976 axh[i] = binWidth/2.;
1977 // now get the coordinate for the point
1978 ax[i] = nominal->GetBinCenter(i+1);
1979 ay[i] = nominal->GetBinContent(i+1);
1980 }
1981 // save the nominal ratio
1982 TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
1983 shapeRatio = (TGraphAsymmErrors*)nominalError->Clone();
1984 nominalError->SetName("shape uncertainty");
1985 TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
1986 nominalCanvas->Divide(2);
1987 nominalCanvas->cd(1);
1988 Style(nominal, kBlack);
1989 Style(nominalError, kYellow, kRatio);
1990 nominalError->SetLineColor(kYellow);
1991 nominalError->SetMarkerColor(kYellow);
1992 nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1993 nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
1994 nominalError->DrawClone("a2");
1995 nominal->DrawCopy("same E1");
1996 Style(AddLegend(gPad));
1997 Style(gPad, "GRID");
1998 Style(nominalCanvas);
1999 // save nominal v2 and systematic errors
2000 TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
2001 nominalIn,
2002 nominalOut,
2003 .56,
2004 "v_{2} with shape uncertainty",
2005 relativeErrorInUp,
2006 relativeErrorInLow,
2007 relativeErrorOutUp,
2008 relativeErrorOutLow));
2009 shapeV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
2010 TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, .56, "v_{2}"));
2011 nominalCanvas->cd(2);
2012 Style(nominalV2, kBlack);
2013 Style(nominalV2Error, kYellow, kV2);
2014 nominalV2Error->SetLineColor(kYellow);
2015 nominalV2Error->SetMarkerColor(kYellow);
2016 nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2017 nominalV2Error->DrawClone("a2");
2018 nominalV2->DrawClone("same E1");
2019 Style(AddLegend(gPad));
2020 Style(gPad, "GRID");
2021 Style(nominalCanvas);
2022 SavePadToPDF(nominalCanvas);
2023 nominalCanvas->Write();
2024 }
2025
18698978 2026 TCanvas* relativeError(new TCanvas("relativeError"," relativeError"));
2027 relativeError->Divide(2);
2028 relativeError->cd(1);
2029 Style(gPad, "GRID");
2030 relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2031 Style(relativeErrorInUp, kBlue, kBar);
2032 Style(relativeErrorInLow, kGreen, kBar);
2033 relativeErrorInUp->DrawCopy("b");
2034 relativeErrorInLow->DrawCopy("same b");
2035 Style(relativeStatisticalErrorIn, kRed);
2036 relativeStatisticalErrorIn->DrawCopy("same");
2037 Style(AddLegend(gPad));
2038 relativeError->cd(2);
2039 Style(gPad, "GRID");
2040 relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2041 Style(relativeErrorOutUp, kBlue, kBar);
2042 Style(relativeErrorOutLow, kGreen, kBar);
2043 relativeErrorOutUp->DrawCopy("b");
2044 relativeErrorOutLow->DrawCopy("same b");
2045 Style(relativeStatisticalErrorOut, kRed);
2046 relativeStatisticalErrorOut->DrawCopy("same");
2047 Style(AddLegend(gPad));
2048
2049 // write the buffered file to disk and close the file
2050 SavePadToPDF(relativeError);
2051 relativeError->Write();
2052 output->Write();
a39e4b2b 2053// output->Close();
18698978 2054}
2055//_____________________________________________________________________________
a39e4b2b 2056 void AliJetFlowTools::DoIntermediateSystematics(
2057 TArrayI* variationsIn, // variantions in plane
2058 TArrayI* variationsOut, // variantions out of plane
2059 TH1D*& relativeErrorInUp, // pointer reference to minimum relative error histogram in plane
2060 TH1D*& relativeErrorInLow, // pointer reference to maximum relative error histogram in plane
2061 TH1D*& relativeErrorOutUp, // pointer reference to minimum relative error histogram out of plane
2062 TH1D*& relativeErrorOutLow, // pointer reference to maximum relative error histogram out of plane
2063 TH1D*& relativeStatisticalErrorIn, // relative systematic error on ratio
2064 TH1D*& relativeStatisticalErrorOut, // relative systematic error on ratio
2065 TH1D*& nominal, // clone of the nominal ratio in / out of plane
2066 TH1D*& nominalIn, // clone of the nominal in plane yield
2067 TH1D*& nominalOut, // clone of the nominal out of plane yield
2068 Int_t columns, // divide the output canvasses in this many columns
2069 Float_t rangeLow, // lower pt range
2070 Float_t rangeUp, // upper pt range
2071 TFile* readMe, // input file name (created by this unfolding class)
2072 TString source // source of the variation
2073 ) const
18698978 2074{
2075 // intermediate systematic check function. first index of supplied array is nominal value
f3ba6c8e 2076 //
18698978 2077 TList* listOfKeys((TList*)readMe->GetListOfKeys());
f3ba6c8e 2078 if(!listOfKeys) {
2079 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2080 return;
2081 }
2082 // check input params
2083 if(variationsIn->GetSize() != variationsOut->GetSize()) {
2084 printf(" > DoSystematics: fatal error, input arrays have different sizes ! < \n ");
2085 return;
2086 }
18698978 2087 TDirectoryFile* defRootDirIn(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(0))->GetName())));
2088 TDirectoryFile* defRootDirOut(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(0))->GetName())));
f3ba6c8e 2089 if(!(defRootDirIn && defRootDirOut)) {
2090 printf(" > DoSystematics: fatal error, couldn't retrieve nominal values ! < \n ");
2091 return;
2092 }
2093 TString defIn(defRootDirIn->GetName());
2094 TString defOut(defRootDirOut->GetName());
18698978 2095
2096 // define lines to make the output prettier
2097 TLine* lineLow(new TLine(rangeLow, 0., rangeLow, 2.));
2098 TLine* lineUp(new TLine(rangeUp, 0., rangeUp, 2.));
2099 lineLow->SetLineColor(11);
2100 lineUp->SetLineColor(11);
2101 lineLow->SetLineWidth(3);
2102 lineUp->SetLineWidth(3);
2103
2104 // define an output histogram with the maximum relative error from this systematic constribution
2105 relativeErrorInUp = new TH1D(Form("max #sigma/|x| from %s", source.Data()), Form("max #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2106 relativeErrorInLow = new TH1D(Form("min #sigma/|x| from %s", source.Data()), Form("min #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2107 relativeErrorOutUp = new TH1D(Form("max #sigma/|x| from %s", source.Data()), Form("max #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2108 relativeErrorOutLow = new TH1D(Form("min #sigma/|x| from %s", source.Data()), Form("min #sigma/|x| from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2109 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2110 relativeErrorInUp->SetBinContent(b+1, 1.);
2111 relativeErrorInUp->SetBinError(b+1, 0.);
2112 relativeErrorOutUp->SetBinContent(b+1, 1.);
2113 relativeErrorOutUp->SetBinError(b+1, .0);
2114 relativeErrorInLow->SetBinContent(b+1, 1.);
2115 relativeErrorInLow->SetBinError(b+1, 0.);
2116 relativeErrorOutLow->SetBinContent(b+1, 1.);
2117 relativeErrorOutLow->SetBinError(b+1, .0);
2118 }
2119 // define an output histogram with the systematic error from this systematic constribution
2120 if(!relativeStatisticalErrorIn && !relativeStatisticalErrorOut) {
2121 relativeStatisticalErrorIn = new TH1D("relative statistical error, in plane", "#sigma/|x|, statistical, in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2122 relativeStatisticalErrorOut = new TH1D("relative statistical error, out of plane", "#sigma/|x|, statistical, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2123 }
2124
f3ba6c8e 2125 // prepare necessary canvasses
18698978 2126 TCanvas* canvasIn(new TCanvas(Form("SYST_%s_PearsonIn", source.Data()), Form("SYST_%s_PearsonIn", source.Data())));
f3ba6c8e 2127 TCanvas* canvasOut(0x0);
18698978 2128 if(fDphiUnfolding) canvasOut = new TCanvas(Form("SYST_%s_PearsonOut", source.Data()), Form("SYST_%s_PearsonOut", source.Data()));
2129 TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas(Form("SYST_%s_RefoldedIn", source.Data()), Form("SYST_%s_RefoldedIn", source.Data())));
f3ba6c8e 2130 TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
18698978 2131 if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas(Form("SYST_%s_RefoldedOut", source.Data()), Form("SYST_%s_RefoldedOut", source.Data()));
2132 TCanvas* canvasSpectraIn(new TCanvas(Form("SYST_%s_SpectraIn", source.Data()), Form("SYST_%s_SpectraIn", source.Data())));
f3ba6c8e 2133 TCanvas* canvasSpectraOut(0x0);
18698978 2134 if(fDphiUnfolding) canvasSpectraOut = new TCanvas(Form("SYST_%s_SpectraOut", source.Data()), Form("SYST_%s_SpectraOut", source.Data()));
f3ba6c8e 2135 TCanvas* canvasRatio(0x0);
18698978 2136 if(fDphiUnfolding) canvasRatio = new TCanvas(Form("SYST_%s_Ratio", source.Data()), Form("SYST_%s_Ratio", source.Data()));
f3ba6c8e 2137 TCanvas* canvasV2(0x0);
18698978 2138 if(fDphiUnfolding) canvasV2 = new TCanvas(Form("SYST_%s_V2", source.Data()), Form("SYST_%s_V2", source.Data()));
2139 TCanvas* canvasMISC(new TCanvas(Form("SYST_%s_MISC", source.Data()), Form("SYST_%s_MISC", source.Data())));
2140 TCanvas* canvasMasterIn(new TCanvas(Form("SYST_%s_defaultIn", source.Data()), Form("SYST_%s_defaultIn", source.Data())));
f3ba6c8e 2141 TCanvas* canvasMasterOut(0x0);
18698978 2142 if(fDphiUnfolding) canvasMasterOut = new TCanvas(Form("SYST_%s_defaultOut", source.Data()), Form("SYST_%s_defaultOut", source.Data()));
f3ba6c8e 2143 (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
a39e4b2b 2144
18698978 2145 TCanvas* canvasProfiles(new TCanvas(Form("SYST_%s_canvasProfiles", source.Data()), Form("SYST_%s_canvasProfiles", source.Data())));
f3ba6c8e 2146 canvasProfiles->Divide(2);
18698978 2147 TProfile* ratioProfile(new TProfile(Form("SYST_%s_ratioProfile", source.Data()), Form("SYST_%s_ratioProfile", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2148 TProfile* v2Profile(new TProfile(Form("SYST_%s_v2Profile", source.Data()), Form("SYST_%s_v2Profile", source.Data()),fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
f3ba6c8e 2149 // get an estimate of the number of outputs and find the default set
5c636e0d 2150
2151 Int_t rows = 1;
2152 columns = variationsIn->GetSize()-1;
2153 (TMath::Floor(variationsIn->GetSize()/(float)columns)+((variationsIn->GetSize()%columns)>0));
f3ba6c8e 2154 canvasIn->Divide(columns, rows);
2155 if(canvasOut) canvasOut->Divide(columns, rows);
2156 canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
2157 if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
2158 canvasSpectraIn->Divide(columns, rows);
2159 if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
2160 if(canvasRatio) canvasRatio->Divide(columns, rows);
2161 if(canvasV2) canvasV2->Divide(columns, rows);
2162 canvasMasterIn->Divide(columns, rows);
2163 if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
2164
2165 // prepare a separate set of canvases to hold the nominal points
18698978 2166 TCanvas* canvasNominalIn(new TCanvas(Form("Nominal_%s_PearsonIn", source.Data()), Form("Nominal_%s_PearsonIn", source.Data())));
f3ba6c8e 2167 TCanvas* canvasNominalOut(0x0);
18698978 2168 if(fDphiUnfolding) canvasNominalOut = new TCanvas(Form("Nominal_%s_PearsonOut", source.Data()), Form("Nominal_%s_PearsonOut", source.Data()));
2169 TCanvas* canvasNominalRatioMeasuredRefoldedIn(new TCanvas(Form("Nominal_%s_RefoldedIn", source.Data()), Form("Nominal_%s_RefoldedIn", source.Data())));
f3ba6c8e 2170 TCanvas* canvasNominalRatioMeasuredRefoldedOut(0x0);
18698978 2171 if(fDphiUnfolding) canvasNominalRatioMeasuredRefoldedOut = new TCanvas(Form("Nominal_%s_RefoldedOut", source.Data()), Form("Nominal_%s_RefoldedOut", source.Data()));
2172 TCanvas* canvasNominalSpectraIn(new TCanvas(Form("Nominal_%s_SpectraIn", source.Data()), Form("Nominal_%s_SpectraIn", source.Data())));
f3ba6c8e 2173 TCanvas* canvasNominalSpectraOut(0x0);
18698978 2174 if(fDphiUnfolding) canvasNominalSpectraOut = new TCanvas(Form("Nominal_%s_SpectraOut", source.Data()), Form("Nominal_%s_SpectraOut", source.Data()));
f3ba6c8e 2175 TCanvas* canvasNominalRatio(0x0);
18698978 2176 if(fDphiUnfolding) canvasNominalRatio = new TCanvas(Form("Nominal_%s_Ratio", source.Data()), Form("Nominal_%s_Ratio", source.Data()));
f3ba6c8e 2177 TCanvas* canvasNominalV2(0x0);
18698978 2178 if(fDphiUnfolding) canvasNominalV2 = new TCanvas(Form("Nominal_%s_V2", source.Data()), Form("Nominal_%s_V2", source.Data()));
2179 TCanvas* canvasNominalMISC(new TCanvas(Form("Nominal_%s_MISC", source.Data()), Form("Nominal_%s_MISC", source.Data())));
2180 TCanvas* canvasNominalMasterIn(new TCanvas(Form("Nominal_%s_defaultIn", source.Data()), Form("Nominal_%s_defaultIn", source.Data())));
f3ba6c8e 2181 TCanvas* canvasNominalMasterOut(0x0);
18698978 2182 if(fDphiUnfolding) canvasNominalMasterOut = new TCanvas(Form("Nominal_%s_defaultOut", source.Data()), Form("Nominal_%s_defaultOut", source.Data()));
f3ba6c8e 2183 (fDphiUnfolding) ? canvasNominalMISC->Divide(4, 2) : canvasNominalMISC->Divide(4, 1);
18698978 2184
2185 canvasNominalSpectraIn->Divide(2);
2186 if(canvasNominalSpectraOut) canvasNominalSpectraOut->Divide(2);
2187
2188 canvasNominalMasterIn->Divide(2);
2189 if(canvasNominalMasterOut) canvasNominalMasterOut->Divide(2);
f3ba6c8e 2190
2191 // extract the default output
2192 TH1D* defaultUnfoldedJetSpectrumIn(0x0);
2193 TH1D* defaultUnfoldedJetSpectrumOut(0x0);
2194 TDirectoryFile* defDirIn = (TDirectoryFile*)defRootDirIn->Get(Form("InPlane___%s", defIn.Data()));
2195 TDirectoryFile* defDirOut = (TDirectoryFile*)defRootDirOut->Get(Form("OutOfPlane___%s", defOut.Data()));
2196 if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", defIn.Data()));
2197 if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", defOut.Data()));
2198 printf(" > succesfully extracted default results < \n");
a39e4b2b 2199
f3ba6c8e 2200 // loop through the directories, only plot the graphs if the deconvolution converged
2201 TDirectoryFile* tempDirIn(0x0);
2202 TDirectoryFile* tempDirOut(0x0);
2203 TDirectoryFile* tempIn(0x0);
2204 TDirectoryFile* tempOut(0x0);
2205 TH1D* unfoldedSpectrumInForRatio(0x0);
2206 TH1D* unfoldedSpectrumOutForRatio(0x0);
2207 for(Int_t i(0), j(-1); i < variationsIn->GetSize(); i++) {
18698978 2208 tempDirIn = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(i))->GetName())));
2209 tempDirOut = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(i))->GetName())));
f3ba6c8e 2210 if(!(tempDirIn && tempDirOut)) {
2211 printf(" > DoSystematics: couldn't get a set of variations < \n");
2212 continue;
2213 }
2214 TString dirNameIn(tempDirIn->GetName());
2215 TString dirNameOut(tempDirOut->GetName());
2216 // try to read the in- and out of plane subdirs
2217 tempIn = (TDirectoryFile*)tempDirIn->Get(Form("InPlane___%s", dirNameIn.Data()));
2218 tempOut = (TDirectoryFile*)tempDirOut->Get(Form("OutOfPlane___%s", dirNameOut.Data()));
2219 j++;
2220 if(tempIn) {
2221 // to see if the unfolding converged try to extract the pearson coefficients
2222 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirNameIn.Data())));
2223 if(pIn) {
2224 printf(" - %s in plane converged \n", dirNameIn.Data());
2225 canvasIn->cd(j);
2226 if(i==0) canvasNominalIn->cd(j);
2227 Style(gPad, "PEARSON");
2228 pIn->DrawCopy("colz");
2229 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirNameIn.Data())));
2230 if(rIn) {
2231 Style(rIn);
2232 printf(" > found RatioRefoldedMeasured < \n");
2233 canvasRatioMeasuredRefoldedIn->cd(j);
2234 if(i==0) canvasNominalRatioMeasuredRefoldedIn->cd(j);
18698978 2235 Style(gPad, "GRID");
f3ba6c8e 2236 rIn->SetFillColor(kRed);
2237 rIn->Draw("ap");
2238 }
2239 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2240 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2241 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirNameIn.Data())));
2242 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirNameIn.Data())));
2243 if(dvector && avalue && rm && eff) {
2244 Style(dvector);
2245 Style(avalue);
2246 Style(rm);
2247 Style(eff);
2248 canvasMISC->cd(1);
2249 if(i==0) canvasNominalMISC->cd(1);
2250 Style(gPad, "SPECTRUM");
2251 dvector->DrawCopy();
2252 canvasMISC->cd(2);
2253 if(i==0) canvasNominalMISC->cd(2);
2254 Style(gPad, "SPECTRUM");
2255 avalue->DrawCopy();
2256 canvasMISC->cd(3);
2257 if(i==0) canvasNominalMISC->cd(3);
2258 Style(gPad, "PEARSON");
2259 rm->DrawCopy("colz");
2260 canvasMISC->cd(4);
2261 if(i==0) canvasNominalMISC->cd(4);
2262 eff->DrawCopy();
2263 } else if(rm && eff) {
2264 Style(rm);
2265 Style(eff);
2266 canvasMISC->cd(3);
2267 if(i==0) canvasNominalMISC->cd(3);
2268 Style(gPad, "PEARSON");
2269 rm->DrawCopy("colz");
2270 canvasMISC->cd(4);
2271 if(i==0) canvasNominalMISC->cd(4);
2272 eff->DrawCopy();
2273 }
2274 }
2275 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirNameIn.Data())));
2276 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirNameIn.Data())));
2277 unfoldedSpectrumInForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2278 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirNameIn.Data())));
2279 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2280 if(defaultUnfoldedJetSpectrumIn) {
2281 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2282 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirNameIn.Data())));
2283 temp->Divide(unfoldedSpectrum);
18698978 2284 // get the absolute relative error
2285 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2286 // check if the error is larger than the current maximum
2287 if( temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorInUp->GetBinContent(b+1)) {
2288 relativeErrorInUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2289 relativeErrorInUp->SetBinError(b+1, 0.);
2290 }
2291 // check if the error is smaller than the current minimum
2292 else if(temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorInLow->GetBinContent(b+1)) {
2293 relativeErrorInLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2294 relativeErrorInLow->SetBinError(b+1, 0.);
2295 }
2296 if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorIn->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2297 }
f3ba6c8e 2298 temp->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2299 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2300 temp->GetYaxis()->SetTitle("ratio");
2301 canvasMasterIn->cd(j);
f3ba6c8e 2302 temp->GetYaxis()->SetRangeUser(0., 2);
18698978 2303 Style(gPad, "GRID");
f3ba6c8e 2304 temp->DrawCopy();
18698978 2305 canvasNominalMasterIn->cd(1);
2306 Style(gPad, "GRID");
2307 if(i > 0 ) {
2308 TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2309 tempSyst->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2310 Style(tempSyst, (EColor)(i+2));
2311 if(i==1) tempSyst->DrawCopy();
2312 else tempSyst->DrawCopy("same");
2313 }
f3ba6c8e 2314 }
2315 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirNameIn.Data())));
2316 canvasSpectraIn->cd(j);
18698978 2317 if(i==0) canvasNominalSpectraIn->cd(1);
f3ba6c8e 2318 Style(gPad);
2319 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2320 unfoldedSpectrum->DrawCopy();
2321 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2322 inputSpectrum->DrawCopy("same");
2323 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2324 refoldedSpectrum->DrawCopy("same");
2325 TLegend* l(AddLegend(gPad));
2326 Style(l);
2327 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2328 Float_t chi(fitStatus->GetBinContent(1));
2329 Float_t pen(fitStatus->GetBinContent(2));
2330 Int_t dof((int)fitStatus->GetBinContent(3));
2331 Float_t beta(fitStatus->GetBinContent(4));
2332 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2333 } else if (fitStatus) { // only available in SVD fit
2334 Int_t reg((int)fitStatus->GetBinContent(1));
2335 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2336 }
18698978 2337 canvasNominalSpectraIn->cd(2);
2338 TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2339 tempSyst->SetTitle(Form("[%s]", dirNameIn.Data()));
2340 Style(tempSyst, (EColor)(i+2));
2341 Style(gPad, "SPECTRUM");
2342 if(i==0) tempSyst->DrawCopy();
2343 else tempSyst->DrawCopy("same");
f3ba6c8e 2344 }
2345 }
2346 if(tempOut) {
2347 TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirNameOut.Data())));
2348 if(pOut) {
2349 printf(" - %s out of plane converged \n", dirNameOut.Data());
2350 canvasOut->cd(j);
2351 if(i==0) canvasNominalOut->cd(j);
2352 Style(gPad, "PEARSON");
2353 pOut->DrawCopy("colz");
2354 TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirNameOut.Data())));
2355 if(rOut) {
2356 Style(rOut);
2357 printf(" > found RatioRefoldedMeasured < \n");
2358 canvasRatioMeasuredRefoldedOut->cd(j);
2359 if(i==0) canvasNominalRatioMeasuredRefoldedOut->cd(j);
18698978 2360 Style(gPad, "GRID");
f3ba6c8e 2361 rOut->SetFillColor(kRed);
2362 rOut->Draw("ap");
2363 }
2364 TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2365 TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2366 TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirNameOut.Data())));
2367 TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirNameOut.Data())));
2368 if(dvector && avalue && rm && eff) {
2369 Style(dvector);
2370 Style(avalue);
2371 Style(rm);
2372 Style(eff);
2373 canvasMISC->cd(5);
2374 if(i==0) canvasNominalMISC->cd(5);
2375 Style(gPad, "SPECTRUM");
2376 dvector->DrawCopy();
2377 canvasMISC->cd(6);
2378 if(i==0) canvasNominalMISC->cd(6);
2379 Style(gPad, "SPECTRUM");
2380 avalue->DrawCopy();
2381 canvasMISC->cd(7);
2382 if(i==0) canvasNominalMISC->cd(7);
2383 Style(gPad, "PEARSON");
2384 rm->DrawCopy("colz");
2385 canvasMISC->cd(8);
2386 if(i==0) canvasNominalMISC->cd(8);
2387 eff->DrawCopy();
2388 } else if(rm && eff) {
2389 Style(rm);
2390 Style(eff);
2391 canvasMISC->cd(7);
2392 if(i==0) canvasNominalMISC->cd(7);
2393 Style(gPad, "PEARSON");
2394 rm->DrawCopy("colz");
2395 canvasMISC->cd(8);
2396 if(i==0) canvasNominalMISC->cd(8);
2397 eff->DrawCopy();
2398 }
2399 }
2400 TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirNameOut.Data())));
2401 TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirNameOut.Data())));
2402 unfoldedSpectrumOutForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2403 TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirNameOut.Data())));
2404 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2405 if(defaultUnfoldedJetSpectrumOut) {
2406 Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
2407 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirNameOut.Data())));
2408 temp->Divide(unfoldedSpectrum);
18698978 2409 // get the absolute relative error
2410 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2411 // check if the error is larger than the current maximum
2412 if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorOutUp->GetBinContent(b+1)) {
2413 relativeErrorOutUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2414 relativeErrorOutUp->SetBinError(b+1, 0.);
2415 }
2416 // check if the error is smaller than the current minimum
2417 else if(temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorOutLow->GetBinContent(b+1)) {
2418 relativeErrorOutLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2419 relativeErrorOutLow->SetBinError(b+1, 0.);
2420 }
2421 if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorOut->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
a39e4b2b 2422 }
18698978 2423 temp->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
f3ba6c8e 2424 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2425 temp->GetYaxis()->SetTitle("ratio");
2426 canvasMasterOut->cd(j);
18698978 2427 temp->GetYaxis()->SetRangeUser(0., 2);
2428 Style(gPad, "GRID");
f3ba6c8e 2429 temp->DrawCopy();
18698978 2430 canvasNominalMasterOut->cd(1);
2431 Style(gPad, "GRID");
2432 if(i > 0 ) {
2433 TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2434 tempSyst->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2435 Style(tempSyst, (EColor)(i+2));
2436 if(i==1) tempSyst->DrawCopy();
2437 else tempSyst->DrawCopy("same");
2438 }
f3ba6c8e 2439 }
2440 TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirNameOut.Data())));
2441 canvasSpectraOut->cd(j);
18698978 2442 if(i==0) canvasNominalSpectraOut->cd(1);
f3ba6c8e 2443 Style(gPad);
2444 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2445 unfoldedSpectrum->DrawCopy();
2446 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2447 inputSpectrum->DrawCopy("same");
2448 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2449 refoldedSpectrum->DrawCopy("same");
2450 TLegend* l(AddLegend(gPad));
2451 Style(l);
2452 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2453 Float_t chi(fitStatus->GetBinContent(1));
2454 Float_t pen(fitStatus->GetBinContent(2));
2455 Int_t dof((int)fitStatus->GetBinContent(3));
2456 Float_t beta(fitStatus->GetBinContent(4));
2457 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2458 } else if (fitStatus) { // only available in SVD fit
2459 Int_t reg((int)fitStatus->GetBinContent(1));
2460 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2461 }
18698978 2462 canvasNominalSpectraOut->cd(2);
2463 TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2464 tempSyst->SetTitle(Form("[%s]", dirNameOut.Data()));
2465 Style(tempSyst, (EColor)(i+2));
2466 Style(gPad, "SPECTRUM");
2467 if(i==0) tempSyst->DrawCopy();
2468 else tempSyst->DrawCopy("same");
f3ba6c8e 2469 }
2470 }
2471 if(canvasRatio && canvasV2) {
2472 canvasRatio->cd(j);
a39e4b2b 2473 if(i==0) {
2474 canvasNominalRatio->cd(j);
2475 if(nominal && nominalIn && nominalOut) {
2476 // if a nominal ratio is requested, delete the placeholder and update the nominal point
2477 delete nominal;
2478 delete nominalIn;
2479 delete nominalOut;
2480 nominalIn = (TH1D*)unfoldedSpectrumInForRatio->Clone("in plane jet yield");
2481 nominalOut = (TH1D*)unfoldedSpectrumOutForRatio->Clone("out of plane jet yield");
2482 nominal = (TH1D*)unfoldedSpectrumInForRatio->Clone("ratio in plane out of plane");
2483 nominal->Divide(nominal, unfoldedSpectrumOutForRatio); // standard root divide for uncorrelated histos
2484 }
2485 }
f3ba6c8e 2486 TGraphErrors* ratioYield(GetRatio(unfoldedSpectrumInForRatio, unfoldedSpectrumOutForRatio, TString(Form("ratio [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2487 Double_t _tempx(0), _tempy(0);
2488 if(ratioYield) {
2489 Style(ratioYield);
2490 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
a39e4b2b 2491 ratioYield->GetPoint(b, _tempx, _tempy);
f3ba6c8e 2492 ratioProfile->Fill(_tempx, _tempy);
2493 }
2494 ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
2495 ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2496 ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
2497 ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2498 ratioYield->SetFillColor(kRed);
2499 ratioYield->Draw("ap");
2500 }
2501 canvasV2->cd(j);
2502 if(i==0) canvasNominalV2->cd(j);
a39e4b2b 2503 TGraphErrors* ratioV2(GetV2(unfoldedSpectrumInForRatio,unfoldedSpectrumOutForRatio, .56, TString(Form("v_{2} [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
f3ba6c8e 2504 if(ratioV2) {
2505 Style(ratioV2);
2506 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
a39e4b2b 2507 ratioV2->GetPoint(b, _tempx, _tempy);
f3ba6c8e 2508 v2Profile->Fill(_tempx, _tempy);
2509
2510 }
2511 v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
2512 v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2513 ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
2514 ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2515 ratioV2->SetFillColor(kRed);
2516 ratioV2->Draw("ap");
2517 }
2518 }
2519 delete unfoldedSpectrumInForRatio;
2520 delete unfoldedSpectrumOutForRatio;
2521 }
f3ba6c8e 2522 // save the canvasses
2523 canvasProfiles->cd(1);
2524 Style(ratioProfile);
2525 ratioProfile->DrawCopy();
2526 canvasProfiles->cd(2);
2527 Style(v2Profile);
2528 v2Profile->DrawCopy();
2529 SavePadToPDF(canvasProfiles);
2530 canvasProfiles->Write();
2531 SavePadToPDF(canvasIn);
2532 canvasIn->Write();
2533 if(canvasOut) {
2534 SavePadToPDF(canvasOut);
2535 canvasOut->Write();
2536 }
2537 SavePadToPDF(canvasRatioMeasuredRefoldedIn);
2538 canvasRatioMeasuredRefoldedIn->Write();
2539 if(canvasRatioMeasuredRefoldedOut) {
2540 SavePadToPDF(canvasRatioMeasuredRefoldedOut);
2541 canvasRatioMeasuredRefoldedOut->Write();
2542 }
2543 SavePadToPDF(canvasSpectraIn);
2544 canvasSpectraIn->Write();
2545 if(canvasSpectraOut) {
2546 SavePadToPDF(canvasSpectraOut);
2547 canvasSpectraOut->Write();
2548 SavePadToPDF(canvasRatio);
2549 canvasRatio->Write();
2550 SavePadToPDF(canvasV2);
2551 canvasV2->Write();
2552 }
2553 SavePadToPDF(canvasMasterIn);
2554 canvasMasterIn->Write();
2555 if(canvasMasterOut) {
2556 SavePadToPDF(canvasMasterOut);
2557 canvasMasterOut->Write();
2558 }
2559 SavePadToPDF(canvasMISC);
2560 canvasMISC->Write();
2561 // save the nomial canvasses
2562 SavePadToPDF(canvasNominalIn);
2563 canvasNominalIn->Write();
2564 if(canvasNominalOut) {
2565 SavePadToPDF(canvasNominalOut);
2566 canvasNominalOut->Write();
2567 }
2568 SavePadToPDF(canvasNominalRatioMeasuredRefoldedIn);
2569 canvasNominalRatioMeasuredRefoldedIn->Write();
2570 if(canvasNominalRatioMeasuredRefoldedOut) {
2571 SavePadToPDF(canvasNominalRatioMeasuredRefoldedOut);
2572 canvasNominalRatioMeasuredRefoldedOut->Write();
2573 }
18698978 2574 canvasNominalSpectraIn->cd(2);
2575 Style(AddLegend(gPad));
f3ba6c8e 2576 SavePadToPDF(canvasNominalSpectraIn);
2577 canvasNominalSpectraIn->Write();
2578 if(canvasNominalSpectraOut) {
18698978 2579 canvasNominalSpectraOut->cd(2);
2580 Style(AddLegend(gPad));
f3ba6c8e 2581 SavePadToPDF(canvasNominalSpectraOut);
2582 canvasNominalSpectraOut->Write();
2583 SavePadToPDF(canvasNominalRatio);
2584 canvasNominalRatio->Write();
2585 SavePadToPDF(canvasNominalV2);
2586 canvasNominalV2->Write();
2587 }
18698978 2588 canvasNominalMasterIn->cd(1);
2589 Style(AddLegend(gPad));
2590 lineUp->DrawClone("same");
2591 lineLow->DrawClone("same");
f3ba6c8e 2592 SavePadToPDF(canvasNominalMasterIn);
2593 canvasNominalMasterIn->Write();
2594 if(canvasNominalMasterOut) {
18698978 2595 canvasNominalMasterOut->cd(1);
2596 Style(AddLegend(gPad));
2597 lineUp->DrawClone("same");
2598 lineLow->DrawClone("same");
f3ba6c8e 2599 SavePadToPDF(canvasNominalMasterOut);
2600 canvasNominalMasterOut->Write();
2601 }
2602 SavePadToPDF(canvasNominalMISC);
2603 canvasNominalMISC->Write();
2604
18698978 2605 // save the relative errors
2606 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2607 relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)-1);
2608 relativeErrorInUp->SetBinError(b+1, 0.);
2609 relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)-1);
2610 relativeErrorOutUp->SetBinError(b+1, .0);
2611 relativeErrorInLow->SetBinContent(b+1, relativeErrorInLow->GetBinContent(b+1)-1);
2612 relativeErrorInLow->SetBinError(b+1, 0.);
2613 relativeErrorOutLow->SetBinContent(b+1, relativeErrorOutLow->GetBinContent(b+1)-1);
2614 relativeErrorOutLow->SetBinError(b+1, .0);
2615 }
2616 relativeErrorInUp->SetYTitle("relative uncertainty");
2617 relativeErrorOutUp->SetYTitle("relative uncertainty");
2618 relativeErrorInLow->SetYTitle("relative uncertainty");
2619 relativeErrorOutLow->SetYTitle("relative uncertainty");
2620 relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2621 relativeErrorInLow->GetYaxis()->SetRangeUser(-1.5, 3.);
2622 relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2623 relativeErrorOutLow->GetYaxis()->SetRangeUser(-1.5, 3.);
2624
2625 canvasNominalMasterIn->cd(2);
2626 Style(gPad, "GRID");
2627 Style(relativeErrorInUp, kBlue, kBar);
2628 Style(relativeErrorInLow, kGreen, kBar);
2629 relativeErrorInUp->DrawCopy("b");
2630 relativeErrorInLow->DrawCopy("same b");
2631 Style(AddLegend(gPad));
2632 SavePadToPDF(canvasNominalMasterIn);
2633 canvasNominalMasterIn->Write();
2634 canvasNominalMasterOut->cd(2);
2635 Style(gPad, "GRID");
2636 Style(relativeErrorOutUp, kBlue, kBar);
2637 Style(relativeErrorOutLow, kGreen, kBar);
2638 relativeErrorOutUp->DrawCopy("b");
2639 relativeErrorOutLow->DrawCopy("same b");
2640 Style(AddLegend(gPad));
2641 SavePadToPDF(canvasNominalMasterOut);
2642 canvasNominalMasterOut->Write();
f3ba6c8e 2643}
2644//_____________________________________________________________________________
87233f72 2645void AliJetFlowTools::PostProcess(TString def, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
d7ec324f 2646{
2647 // go through the output file and perform post processing routines
2648 // can either be performed in one go with the unfolding, or at a later stage
53547ff2 2649 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
d7ec324f 2650 TFile readMe(in.Data(), "READ"); // open file read-only
2651 if(readMe.IsZombie()) {
2652 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
2653 return;
2654 }
2655 printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
2656 readMe.ls();
2657 TList* listOfKeys((TList*)readMe.GetListOfKeys());
2658 if(!listOfKeys) {
2659 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2660 return;
2661 }
2662 // prepare necessary canvasses
53547ff2 2663 TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
5e11c41c 2664 TCanvas* canvasOut(0x0);
486fb24e 2665 if(fDphiUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
53547ff2 2666 TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
5e11c41c 2667 TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
486fb24e 2668 if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
53547ff2 2669 TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
5e11c41c 2670 TCanvas* canvasSpectraOut(0x0);
486fb24e 2671 if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
5e11c41c 2672 TCanvas* canvasRatio(0x0);
486fb24e 2673 if(fDphiUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
5e11c41c 2674 TCanvas* canvasV2(0x0);
486fb24e 2675 if(fDphiUnfolding) canvasV2 = new TCanvas("V2", "V2");
53547ff2
RAB		 2676 TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
2677 TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
5e11c41c 2678 TCanvas* canvasMasterOut(0x0);
486fb24e 2679 if(fDphiUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
2680 (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
d7ec324f 2681 TDirectoryFile* defDir(0x0);
87233f72 2682 TCanvas* canvasProfiles(new TCanvas("canvasProfiles", "canvasProfiles"));
2683 canvasProfiles->Divide(2);
2684 TProfile* ratioProfile(new TProfile("ratioProfile", "ratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2685 TProfile* v2Profile(new TProfile("v2Profile", "v2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
d7ec324f 2686 // get an estimate of the number of outputs and find the default set
2687 Int_t cacheMe(0);
2688 for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
2689 if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
2690 if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
2691 cacheMe++;
2692 }
2693 }
f3ba6c8e 2694 Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%columns)>0));
53547ff2 2695 canvasIn->Divide(columns, rows);
5e11c41c 2696 if(canvasOut) canvasOut->Divide(columns, rows);
53547ff2 2697 canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
5e11c41c 2698 if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
53547ff2 2699 canvasSpectraIn->Divide(columns, rows);
5e11c41c 2700 if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
2701 if(canvasRatio) canvasRatio->Divide(columns, rows);
2702 if(canvasV2) canvasV2->Divide(columns, rows);
d7ec324f 2703
53547ff2 2704 canvasMasterIn->Divide(columns, rows);
5e11c41c 2705 if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
d7ec324f 2706 // extract the default output
ab474fb0 2707 TH1D* defaultUnfoldedJetSpectrumIn(0x0);
2708 TH1D* defaultUnfoldedJetSpectrumOut(0x0);
d7ec324f 2709 if(defDir) {
2710 TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
2711 TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
ab474fb0 2712 if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
2713 if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
d7ec324f 2714 printf(" > succesfully extracted default results < \n");
2715 }
a39e4b2b 2716
d7ec324f 2717 // loop through the directories, only plot the graphs if the deconvolution converged
2718 TDirectoryFile* tempDir(0x0);
2719 TDirectoryFile* tempIn(0x0);
2720 TDirectoryFile* tempOut(0x0);
2721 for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
ab474fb0 2722 // read the directory by its name
d7ec324f 2723 tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
2724 if(!tempDir) continue;
2725 TString dirName(tempDir->GetName());
ab474fb0 2726 // try to read the in- and out of plane subdirs
d7ec324f 2727 tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
2728 tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
2729 j++;
486fb24e 2730 if(!tempIn) { // attempt to get MakeAU output
2731 TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
2732 TCanvas* tempPearson(new TCanvas(Form("pearson_%s", dirName.Data()), Form("pearson_%s", dirName.Data())));
2733 tempPearson->Divide(4,2);
2734 TCanvas* tempRatio(new TCanvas(Form("ratio_%s", dirName.Data()), Form("ratio_%s", dirName.Data())));
2735 tempRatio->Divide(4,2);
2736 TCanvas* tempResult(new TCanvas(Form("result_%s", dirName.Data()), Form("result_%s", dirName.Data())));
2737 tempResult->Divide(4,2);
2738 for(Int_t q(0); q < 8; q++) {
2739 tempIn = (TDirectoryFile*)tempDir->Get(Form("%s___%s", stringArray[q].Data(), dirName.Data()));
2740 if(tempIn) {
a39e4b2b 2741 // to see if the unfolding converged try to extract the pearson coefficients
2742 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
2743 if(pIn) {
486fb24e 2744 printf(" - %s in plane converged \n", dirName.Data());
a39e4b2b 2745 tempPearson->cd(1+q);
486fb24e 2746 Style(gPad, "PEARSON");
a39e4b2b 2747 pIn->DrawCopy("colz");
2748 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2749 if(rIn) {
2750 Style(rIn);
2751 printf(" > found RatioRefoldedMeasured < \n");
2752 tempRatio->cd(q+1);
2753 rIn->SetFillColor(kRed);
2754 rIn->Draw("ap");
2755 }
2756 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2757 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2758 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
2759 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
2760 if(dvector && avalue && rm && eff) {
2761 Style(dvector);
2762 Style(avalue);
2763 Style(rm);
2764 Style(eff);
2765 canvasMISC->cd(1);
2766 Style(gPad, "SPECTRUM");
2767 dvector->DrawCopy();
2768 canvasMISC->cd(2);
2769 Style(gPad, "SPECTRUM");
2770 avalue->DrawCopy();
2771 canvasMISC->cd(3);
2772 Style(gPad, "PEARSON");
2773 rm->DrawCopy("colz");
2774 canvasMISC->cd(4);
2775 eff->DrawCopy();
2776 } else if(rm && eff) {
2777 Style(rm);
2778 Style(eff);
2779 canvasMISC->cd(3);
2780 Style(gPad, "PEARSON");
2781 rm->DrawCopy("colz");
2782 canvasMISC->cd(4);
2783 eff->DrawCopy();
2784 }
486fb24e 2785 }
a39e4b2b 2786 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
2787 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
2788 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
2789 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2790 if(defaultUnfoldedJetSpectrumIn) {
2791 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2792 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
2793 temp->Divide(unfoldedSpectrum);
2794 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
2795 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2796 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
2797 canvasMasterIn->cd(j);
2798 temp->GetYaxis()->SetRangeUser(0., 2);
2799 temp->DrawCopy();
2800 }
2801 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
2802 tempResult->cd(q+1);
2803 Style(gPad);
2804 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2805 unfoldedSpectrum->DrawCopy();
2806 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2807 inputSpectrum->DrawCopy("same");
2808 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2809 refoldedSpectrum->DrawCopy("same");
2810 TLegend* l(AddLegend(gPad));
2811 Style(l);
2812 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2813 Float_t chi(fitStatus->GetBinContent(1));
2814 Float_t pen(fitStatus->GetBinContent(2));
2815 Int_t dof((int)fitStatus->GetBinContent(3));
2816 Float_t beta(fitStatus->GetBinContent(4));
2817 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2818 } else if (fitStatus) { // only available in SVD fit
2819 Int_t reg((int)fitStatus->GetBinContent(1));
2820 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2821 }
2822 }
486fb24e 2823 }
2824 }
2825 }
d7ec324f 2826 if(tempIn) {
2827 // to see if the unfolding converged try to extract the pearson coefficients
2828 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
2829 if(pIn) {
2830 printf(" - %s in plane converged \n", dirName.Data());
2831 canvasIn->cd(j);
2832 Style(gPad, "PEARSON");
2833 pIn->DrawCopy("colz");
2834 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2835 if(rIn) {
53547ff2 2836 Style(rIn);
d7ec324f 2837 printf(" > found RatioRefoldedMeasured < \n");
2838 canvasRatioMeasuredRefoldedIn->cd(j);
ab474fb0 2839 rIn->SetFillColor(kRed);
2840 rIn->Draw("ap");
d7ec324f 2841 }
2842 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2843 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2844 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
2845 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
2846 if(dvector && avalue && rm && eff) {
53547ff2
RAB		 2847 Style(dvector);
2848 Style(avalue);
2849 Style(rm);
2850 Style(eff);
d7ec324f 2851 canvasMISC->cd(1);
2852 Style(gPad, "SPECTRUM");
2853 dvector->DrawCopy();
2854 canvasMISC->cd(2);
2855 Style(gPad, "SPECTRUM");
2856 avalue->DrawCopy();
2857 canvasMISC->cd(3);
2858 Style(gPad, "PEARSON");
2859 rm->DrawCopy("colz");
2860 canvasMISC->cd(4);
2861 eff->DrawCopy();
53547ff2
RAB		 2862 } else if(rm && eff) {
2863 Style(rm);
2864 Style(eff);
2865 canvasMISC->cd(3);
2866 Style(gPad, "PEARSON");
2867 rm->DrawCopy("colz");
2868 canvasMISC->cd(4);
2869 eff->DrawCopy();
d7ec324f 2870 }
2871 }
2872 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
2873 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
2874 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
2875 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
ab474fb0 2876 if(defaultUnfoldedJetSpectrumIn) {
2877 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2878 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
d7ec324f 2879 temp->Divide(unfoldedSpectrum);
2880 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
f3ba6c8e 2881 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
d7ec324f 2882 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
2883 canvasMasterIn->cd(j);
ab474fb0 2884 temp->GetYaxis()->SetRangeUser(0., 2);
d7ec324f 2885 temp->DrawCopy();
2886 }
2887 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
2888 canvasSpectraIn->cd(j);
2889 Style(gPad);
53547ff2 2890 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
d7ec324f 2891 unfoldedSpectrum->DrawCopy();
53547ff2 2892 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
d7ec324f 2893 inputSpectrum->DrawCopy("same");
53547ff2 2894 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
d7ec324f 2895 refoldedSpectrum->DrawCopy("same");
2896 TLegend* l(AddLegend(gPad));
53547ff2
RAB		 2897 Style(l);
2898 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2899 Float_t chi(fitStatus->GetBinContent(1));
2900 Float_t pen(fitStatus->GetBinContent(2));
d7ec324f 2901 Int_t dof((int)fitStatus->GetBinContent(3));
53547ff2
RAB		 2902 Float_t beta(fitStatus->GetBinContent(4));
2903 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2904 } else if (fitStatus) { // only available in SVD fit
2905 Int_t reg((int)fitStatus->GetBinContent(1));
2906 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
d7ec324f 2907 }
2908 }
2909 }
2910 if(tempOut) {
2911 TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
2912 if(pOut) {
2913 printf(" - %s out of plane converged \n", dirName.Data());
2914 canvasOut->cd(j);
2915 Style(gPad, "PEARSON");
2916 pOut->DrawCopy("colz");
2917 TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2918 if(rOut) {
53547ff2 2919 Style(rOut);
d7ec324f 2920 printf(" > found RatioRefoldedMeasured < \n");
2921 canvasRatioMeasuredRefoldedOut->cd(j);
ab474fb0 2922 rOut->SetFillColor(kRed);
2923 rOut->Draw("ap");
d7ec324f 2924 }
2925 TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2926 TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2927 TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
2928 TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
2929 if(dvector && avalue && rm && eff) {
53547ff2
RAB		 2930 Style(dvector);
2931 Style(avalue);
2932 Style(rm);
2933 Style(eff);
d7ec324f 2934 canvasMISC->cd(5);
2935 Style(gPad, "SPECTRUM");
2936 dvector->DrawCopy();
2937 canvasMISC->cd(6);
2938 Style(gPad, "SPECTRUM");
2939 avalue->DrawCopy();
2940 canvasMISC->cd(7);
2941 Style(gPad, "PEARSON");
2942 rm->DrawCopy("colz");
2943 canvasMISC->cd(8);
2944 eff->DrawCopy();
53547ff2
RAB		 2945 } else if(rm && eff) {
2946 Style(rm);
2947 Style(eff);
549b5f40 2948 canvasMISC->cd(7);
53547ff2
RAB		 2949 Style(gPad, "PEARSON");
2950 rm->DrawCopy("colz");
549b5f40 2951 canvasMISC->cd(8);
53547ff2 2952 eff->DrawCopy();
d7ec324f 2953 }
2954 }
2955 TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
2956 TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
2957 TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
2958 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
ab474fb0 2959 if(defaultUnfoldedJetSpectrumOut) {
2960 Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
2961 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirName.Data())));
d7ec324f 2962 temp->Divide(unfoldedSpectrum);
2963 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
f3ba6c8e 2964 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
d7ec324f 2965 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
2966 canvasMasterOut->cd(j);
ab474fb0 2967 temp->GetYaxis()->SetRangeUser(0., 2.);
d7ec324f 2968 temp->DrawCopy();
2969 }
2970 TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
2971 canvasSpectraOut->cd(j);
2972 Style(gPad);
53547ff2 2973 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
d7ec324f 2974 unfoldedSpectrum->DrawCopy();
53547ff2 2975 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
d7ec324f 2976 inputSpectrum->DrawCopy("same");
53547ff2 2977 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
d7ec324f 2978 refoldedSpectrum->DrawCopy("same");
2979 TLegend* l(AddLegend(gPad));
53547ff2
RAB		 2980 Style(l);
2981 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2982 Float_t chi(fitStatus->GetBinContent(1));
2983 Float_t pen(fitStatus->GetBinContent(2));
2984 Int_t dof((int)fitStatus->GetBinContent(3));
2985 Float_t beta(fitStatus->GetBinContent(4));
2986 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2987 } else if (fitStatus) { // only available in SVD fit
2988 Int_t reg((int)fitStatus->GetBinContent(1));
2989 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
d7ec324f 2990 }
2991 }
2992 }
5e11c41c 2993 if(canvasRatio && canvasV2) {
2994 canvasRatio->cd(j);
2995 TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
87233f72 2996 Double_t _tempx(0), _tempy(0);
5e11c41c 2997 if(ratioYield) {
2998 Style(ratioYield);
87233f72 2999 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
a39e4b2b 3000 ratioYield->GetPoint(b, _tempx, _tempy);
87233f72 3001 ratioProfile->Fill(_tempx, _tempy);
3002 }
9f892925 3003 ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
87233f72 3004 ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
9f892925 3005 ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
87233f72 3006 ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
ab474fb0 3007 ratioYield->SetFillColor(kRed);
3008 ratioYield->Draw("ap");
5e11c41c 3009 }
3010 canvasV2->cd(j);
3011 TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
3012 if(ratioV2) {
3013 Style(ratioV2);
87233f72 3014 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
a39e4b2b 3015 ratioV2->GetPoint(b, _tempx, _tempy);
87233f72 3016 v2Profile->Fill(_tempx, _tempy);
3017
3018 }
9f892925 3019 v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
87233f72 3020 v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
5e11c41c 3021 ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
87233f72 3022 ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
ab474fb0 3023 ratioV2->SetFillColor(kRed);
3024 ratioV2->Draw("ap");
5e11c41c 3025 }
d7ec324f 3026 }
3027 }
3028 TFile output(out.Data(), "RECREATE");
87233f72 3029 canvasProfiles->cd(1);
3030 Style(ratioProfile);
3031 ratioProfile->DrawCopy();
3032 canvasProfiles->cd(2);
3033 Style(v2Profile);
3034 v2Profile->DrawCopy();
3035 SavePadToPDF(canvasProfiles);
3036 canvasProfiles->Write();
512ced40 3037 SavePadToPDF(canvasIn);
d7ec324f 3038 canvasIn->Write();
5e11c41c 3039 if(canvasOut) {
3040 SavePadToPDF(canvasOut);
3041 canvasOut->Write();
3042 }
512ced40 3043 SavePadToPDF(canvasRatioMeasuredRefoldedIn);
d7ec324f 3044 canvasRatioMeasuredRefoldedIn->Write();
5e11c41c 3045 if(canvasRatioMeasuredRefoldedOut) {
3046 SavePadToPDF(canvasRatioMeasuredRefoldedOut);
3047 canvasRatioMeasuredRefoldedOut->Write();
3048 }
512ced40 3049 SavePadToPDF(canvasSpectraIn);
d7ec324f 3050 canvasSpectraIn->Write();
5e11c41c 3051 if(canvasSpectraOut) {
3052 SavePadToPDF(canvasSpectraOut);
3053 canvasSpectraOut->Write();
3054 SavePadToPDF(canvasRatio);
3055 canvasRatio->Write();
3056 SavePadToPDF(canvasV2);
3057 canvasV2->Write();
3058 }
512ced40 3059 SavePadToPDF(canvasMasterIn);
d7ec324f 3060 canvasMasterIn->Write();
5e11c41c 3061 if(canvasMasterOut) {
3062 SavePadToPDF(canvasMasterOut);
3063 canvasMasterOut->Write();
3064 }
512ced40 3065 SavePadToPDF(canvasMISC);
d7ec324f 3066 canvasMISC->Write();
3067 output.Write();
3068 output.Close();
3069}
3070//_____________________________________________________________________________
4292ca60 3071Bool_t AliJetFlowTools::SetRawInput (
3072 TH2D* detectorResponse, // detector response matrix
3073 TH1D* jetPtIn, // in plane jet spectrum
3074 TH1D* jetPtOut, // out of plane jet spectrum
3075 TH1D* dptIn, // in plane delta pt distribution
3076 TH1D* dptOut, // out of plane delta pt distribution
3077 Int_t eventCount) {
3078 // set input histograms manually
3079 fDetectorResponse = detectorResponse;
3080 fSpectrumIn = jetPtIn;
3081 fSpectrumOut = jetPtOut;
3082 fDptInDist = dptIn;
3083 fDptOutDist = dptOut;
3084 fRawInputProvided = kTRUE;
3085 // check if all data is provided
3086 if(!fDetectorResponse) {
3087 printf(" fDetectorResponse not found \n ");
3088 return kFALSE;
3089 }
3090 // check if the pt bin for true and rec have been set
3091 if(!fBinsTrue || !fBinsRec) {
3092 printf(" No true or rec bins set, please set binning ! \n");
3093 return kFALSE;
3094 }
3095 if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
a39e4b2b 3096 // procedures, these profiles will be nonsensical, user is responsible
4292ca60 3097 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3098 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3099 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3100 }
3101 // normalize spectra to event count if requested
3102 if(fNormalizeSpectra) {
3103 fEventCount = eventCount;
3104 if(fEventCount > 0) {
3105 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
3106 fSpectrumOut->Sumw2();
3107 fSpectrumIn->Scale(1./((double)fEventCount));
3108 fSpectrumOut->Scale(1./((double)fEventCount));
3109 }
3110 }
3111 if(!fNormalizeSpectra && fEventCount > 0) {
3112 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
3113 fSpectrumOut->Sumw2();
3114 fSpectrumIn->Scale(1./((double)fEventCount));
3115 fSpectrumOut->Scale(1./((double)fEventCount));
3116 }
3117 fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
3118 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
f3ba6c8e 3119 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
3120 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
4292ca60 3121 fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
3122 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
f3ba6c8e 3123 fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
3124 fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
4292ca60 3125
3126 return kTRUE;
3127}
3128//_____________________________________________________________________________
3129TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
dc1455ee 3130{
ef12d5a5 3131 // return ratio of h1 / h2
3132 // histograms can have different binning. errors are propagated as uncorrelated
20abfcc4 3133 if(!(h1 && h2) ) {
3134 printf(" GetRatio called with NULL argument(s) \n ");
3135 return 0x0;
3136 }
ad04a83c 3137 Int_t j(0);
4292ca60 3138 TGraphErrors *gr = new TGraphErrors();
f3ba6c8e 3139 gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
ef12d5a5 3140 Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
7bf39895 3141 TH1* dud((TH1*)h1->Clone("dud"));
486fb24e 3142 dud->Sumw2();
3143 h1->Sumw2();
3144 h2->Sumw2();
3145 if(!dud->Divide(h1, h2)) {
3146 printf(" > ROOT failed to divide two histogams, dividing manually \n < ");
7bf39895 3147 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3148 binCent = h1->GetXaxis()->GetBinCenter(i);
3149 if(xmax > 0. && binCent > xmax) continue;
3150 j = h2->FindBin(binCent);
3151 binWidth = h1->GetXaxis()->GetBinWidth(i);
3152 if(h2->GetBinContent(j) > 0.) {
3153 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
486fb24e 3154 /* original propagation of uncertainty changed 08012014
dc1455ee 3155 Double_t A = 1./h2->GetBinContent(j)*h1->GetBinError(i);
4292ca60 3156 Double_t B = 0.;
3157 if(h2->GetBinError(j)>0.) {
3158 B = -1.*h1->GetBinContent(i)/(h2->GetBinContent(j)*h2->GetBinContent(j))*h2->GetBinError(j);
3159 error2 = A*A + B*B;
ab474fb0 3160 } else error2 = A*A; */
7bf39895 3161 Double_t A = h1->GetBinError(i)/h1->GetBinContent(i);
3162 Double_t B = h2->GetBinError(i)/h2->GetBinContent(i);
3163 error2 = ratio*ratio*A*A+ratio*ratio*B*B;
3164 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
486fb24e 3165 gr->SetPoint(gr->GetN(), binCent, ratio);
3166 gr->SetPointError(gr->GetN()-1, 0.5*binWidth, error2);
7bf39895 3167 }
3168 }
3169 } else {
486fb24e 3170 printf( " > using ROOT to divide two histograms < \n");
7bf39895 3171 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3172 binCent = dud->GetXaxis()->GetBinCenter(i);
3173 if(xmax > 0. && binCent > xmax) continue;
3174 binWidth = dud->GetXaxis()->GetBinWidth(i);
3175 gr->SetPoint(gr->GetN(),binCent,dud->GetBinContent(i));
3176 gr->SetPointError(gr->GetN()-1,0.5*binWidth,dud->GetBinError(i));
dc1455ee 3177 }
3178 }
7bf39895 3179
ad04a83c 3180 if(appendFit) {
4292ca60 3181 TF1* fit(new TF1("lin", "pol0", 10, 100));
ad04a83c 3182 gr->Fit(fit);
3183 }
4292ca60 3184 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
7bf39895 3185 if(dud) delete dud;
dc1455ee 3186 return gr;
3187}
3188//_____________________________________________________________________________
ef12d5a5 3189TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
3190{
3191 // get v2 from difference of in plane, out of plane yield
3192 // h1 must hold the in-plane yield, h2 holds the out of plane yield
549b5f40
RAB		 3193 // different binning is allowed but will mean that the error
3194 // propagation is unreliable
ef12d5a5 3195 // r is the event plane resolution for the chosen centrality
3196 if(!(h1 && h2) ) {
3197 printf(" GetV2 called with NULL argument(s) \n ");
3198 return 0x0;
3199 }
3200 Int_t j(0);
3201 TGraphErrors *gr = new TGraphErrors();
f3ba6c8e 3202 gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
ef12d5a5 3203 Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
3204 Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
3205 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3206 binCent = h1->GetXaxis()->GetBinCenter(i);
3207 j = h2->FindBin(binCent);
3208 binWidth = h1->GetXaxis()->GetBinWidth(i);
3209 if(h2->GetBinContent(j) > 0.) {
3210 in = h1->GetBinContent(i);
3211 ein = h1->GetBinError(i);
3212 out = h2->GetBinContent(j);
3213 eout = h2->GetBinError(j);
3214 ratio = pre*((in-out)/(in+out));
549b5f40 3215 error2 =((r*4.)/(TMath::Pi()))*((4.*out*out/(TMath::Power(in+out, 4)))*ein*ein+(4.*in*in/(TMath::Power(in+out, 4)))*eout*eout);
ef12d5a5 3216 if(error2 > 0) error2 = TMath::Sqrt(error2);
3217 gr->SetPoint(gr->GetN(),binCent,ratio);
3218 gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
3219 }
3220 }
3221 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
3222 return gr;
3223}
3224//_____________________________________________________________________________
a39e4b2b 3225TGraphAsymmErrors* AliJetFlowTools::GetV2WithSystematicErrors(
3226 TH1* h1, TH1* h2, Double_t r, TString name,
3227 TH1* relativeErrorInUp,
3228 TH1* relativeErrorInLow,
3229 TH1* relativeErrorOutUp,
3230 TH1* relativeErrorOutLow) const
3231{
3232 // get v2 with asymmetric systematic error
3233 // note that this is ONLY the systematic error, no statistical error!
3234 TGraphErrors* tempV2(GetV2(h1, h2, r, name));
3235 Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
3236 Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
3237 Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
3238 Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
3239 Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
3240 Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
3241 Double_t in(0.), out(0.), einUp(0.), einLow(0.), eoutUp(0.), eoutLow(0.), error2Up(0.), error2Low(0.);
3242 // loop through the bins and do error propagation
3243 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3244 // extract the absolute errors
3245 in = h1->GetBinContent(i+1);
3246 einUp = in*relativeErrorInUp->GetBinContent(i+1);
3247 einLow = in*relativeErrorInLow->GetBinContent(1+i);
3248 out = h2->GetBinContent(i+1);
3249 eoutUp = out*relativeErrorOutUp->GetBinContent(1+i);
3250 eoutLow = out*relativeErrorOutLow->GetBinContent(1+i);
3251 // get the error squared
3252 error2Up =((r*4.)/(TMath::Pi()))*((4.*out*out/(TMath::Power(in+out, 4)))*einUp*einUp+(4.*in*in/(TMath::Power(in+out, 4)))*eoutUp*eoutUp);
3253 error2Low =((r*4.)/(TMath::Pi()))*((4.*out*out/(TMath::Power(in+out, 4)))*einLow*einLow+(4.*in*in/(TMath::Power(in+out, 4)))*eoutLow*eoutLow);
3254 if(error2Up > 0) error2Up = TMath::Sqrt(error2Up);
3255 if(error2Low > 0) error2Low = TMath::Sqrt(error2Low);
3256 // set the errors
3257 ayh[i] = error2Up;
3258 ayl[i] = error2Low;
3259 // get the bin width (which is the 'error' on x)
3260 Double_t binWidth(h1->GetBinWidth(i+1));
3261 axl[i] = binWidth/2.;
3262 axh[i] = binWidth/2.;
3263 // now get the coordinate for the poin
3264 tempV2->GetPoint(i, ax[i], ay[i]);
3265 }
3266 // save the nominal ratio
3267 TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
3268 nominalError->SetName("v_{2} with shape uncertainty");
3269 // do some memory management
3270 delete tempV2;
3271 delete[] ax;
3272 delete[] ay;
3273 delete[] axh;
3274 delete[] axl;
3275 delete[] ayh;
3276 delete[] ayl;
3277
3278 return nominalError;
3279}
3280//_____________________________________________________________________________
549b5f40
RAB		 3281void AliJetFlowTools::WriteObject(TObject* object, TString suffix, Bool_t kill) {
3282 // write object, if a unique identifier is given the object is cloned
3283 // and the clone is saved. setting kill to true will delete the original obect from the heap
4292ca60 3284 if(!object) {
3285 printf(" > WriteObject:: called with NULL arguments \n ");
3286 return;
549b5f40
RAB		 3287 } else if(!strcmp("", suffix.Data())) object->Write();
3288 else {
3289 TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
3290 newObject->Write();
3291 }
3292 if(kill) delete object;
4292ca60 3293}
3294//_____________________________________________________________________________
3295TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
3296 // construt a delta pt response matrix from supplied dpt distribution
3297 // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
3298 // do a weighted rebinning to a (coarser) dpt distribution
3299 // be careful with the binning of the dpt response: it should be equal to that
3300 // of the response matrix, otherwise dpt and response matrices cannot be multiplied
3301 //
3302 // the response matrix will be square and have the same binning
3303 // (min, max and granularity) of the input histogram
3304 Int_t bins(dpt->GetXaxis()->GetNbins()); // number of bins, will also be no of rows, columns
3305 Double_t _bins[bins+1]; // prepare array with bin borders
3306 for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
3307 _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1); // get upper edge
3308 TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
3309 for(Int_t j(0); j < bins+1; j++) { // loop on pt true slices j
3310 Bool_t skip = kFALSE;
3311 for(Int_t k(0); k < bins+1; k++) { // loop on pt gen slices k
3312 (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
3313 if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
3314 }
3315 }
3316 return res;
3317}
ef12d5a5 3318//_____________________________________________________________________________
3319TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
3320 if(!binsTrue || !binsRec) {
3321 printf(" > GetUnityResponse:: function called with NULL arguments < \n");
3322 return 0x0;
3323 }
3324 TString name(Form("unityResponse_%s", suffix.Data()));
3325 TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
3326 for(Int_t i(0); i < binsTrue->GetSize(); i++) {
3327 for(Int_t j(0); j < binsRec->GetSize(); j++) {
3328 if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
3329 }
3330 }
3331 return unity;
3332}
4292ca60 3333//_____________________________________________________________________________
549b5f40 3334void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
4292ca60 3335 // save configuration parameters to histogram
549b5f40 3336 TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
ef12d5a5 3337 summary->SetBinContent(1, fBetaIn);
3338 summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
3339 summary->SetBinContent(2, fBetaOut);
3340 summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
3341 summary->SetBinContent(3, fCentralityBin);
3342 summary->GetXaxis()->SetBinLabel(3, "fCentralityBin");
3343 summary->SetBinContent(4, (int)convergedIn);
3344 summary->GetXaxis()->SetBinLabel(4, "convergedIn");
3345 summary->SetBinContent(5, (int)convergedOut);
3346 summary->GetXaxis()->SetBinLabel(5, "convergedOut");
3347 summary->SetBinContent(6, (int)fAvoidRoundingError);
3348 summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
3349 summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
3350 summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
3351 summary->SetBinContent(8, (int)fPrior);
3352 summary->GetXaxis()->SetBinLabel(8, "fPrior");
3353 summary->SetBinContent(9, fSVDRegIn);
3354 summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
3355 summary->SetBinContent(10, fSVDRegOut);
3356 summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
3357 summary->SetBinContent(11, (int)fSVDToy);
3358 summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
3359 summary->SetBinContent(12, fJetRadius);
3360 summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
3361 summary->SetBinContent(13, (int)fNormalizeSpectra);
3362 summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
549b5f40
RAB		 3363 summary->SetBinContent(14, (int)fSmoothenPrior);
3364 summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
ef12d5a5 3365 summary->SetBinContent(15, (int)fTestMode);
3366 summary->GetXaxis()->SetBinLabel(15, "fTestMode");
3367 summary->SetBinContent(16, (int)fUseDetectorResponse);
3368 summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
549b5f40
RAB		 3369 summary->SetBinContent(17, fBayesianIterIn);
3370 summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
3371 summary->SetBinContent(18, fBayesianIterOut);
3372 summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
3373 summary->SetBinContent(19, fBayesianSmoothIn);
3374 summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
3375 summary->SetBinContent(20, fBayesianSmoothOut);
3376 summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
4292ca60 3377}
3378//_____________________________________________________________________________
3379void AliJetFlowTools::ResetAliUnfolding() {
3380 // ugly function: reset all unfolding parameters
3381 TVirtualFitter* fitter(TVirtualFitter::GetFitter());
3382 if(fitter) {
3383 printf(" > Found fitter, will delete it < \n");
3384 delete fitter;
3385 }
d7ec324f 3386 if(gMinuit) {
3387 printf(" > Found gMinuit, will re-create it < \n");
3388 delete gMinuit;
3389 gMinuit = new TMinuit;
3390 }
4292ca60 3391 AliUnfolding::fgCorrelationMatrix = 0;
3392 AliUnfolding::fgCorrelationMatrixSquared = 0;
3393 AliUnfolding::fgCorrelationCovarianceMatrix = 0;
3394 AliUnfolding::fgCurrentESDVector = 0;
3395 AliUnfolding::fgEntropyAPriori = 0;
3396 AliUnfolding::fgEfficiency = 0;
3397 AliUnfolding::fgUnfoldedAxis = 0;
3398 AliUnfolding::fgMeasuredAxis = 0;
3399 AliUnfolding::fgFitFunction = 0;
3400 AliUnfolding::fgMaxInput = -1;
3401 AliUnfolding::fgMaxParams = -1;
3402 AliUnfolding::fgOverflowBinLimit = -1;
3403 AliUnfolding::fgRegularizationWeight = 10000;
3404 AliUnfolding::fgSkipBinsBegin = 0;
3405 AliUnfolding::fgMinuitStepSize = 0.1;
3406 AliUnfolding::fgMinuitPrecision = 1e-6;
3407 AliUnfolding::fgMinuitMaxIterations = 1000000;
3408 AliUnfolding::fgMinuitStrategy = 1.;
3409 AliUnfolding::fgMinimumInitialValue = kFALSE;
3410 AliUnfolding::fgMinimumInitialValueFix = -1;
3411 AliUnfolding::fgNormalizeInput = kFALSE;
3412 AliUnfolding::fgNotFoundEvents = 0;
3413 AliUnfolding::fgSkipBin0InChi2 = kFALSE;
3414 AliUnfolding::fgBayesianSmoothing = 1;
3415 AliUnfolding::fgBayesianIterations = 10;
3416 AliUnfolding::fgDebug = kFALSE;
3417 AliUnfolding::fgCallCount = 0;
3418 AliUnfolding::fgPowern = 5;
3419 AliUnfolding::fChi2FromFit = 0.;
3420 AliUnfolding::fPenaltyVal = 0.;
3421 AliUnfolding::fAvgResidual = 0.;
3422 AliUnfolding::fgPrintChi2Details = 0;
3423 AliUnfolding::fgCanvas = 0;
3424 AliUnfolding::fghUnfolded = 0;
3425 AliUnfolding::fghCorrelation = 0;
3426 AliUnfolding::fghEfficiency = 0;
3427 AliUnfolding::fghMeasured = 0;
3428 AliUnfolding::SetMinuitStepSize(1.);
3429 AliUnfolding::SetMinuitPrecision(1e-6);
3430 AliUnfolding::SetMinuitMaxIterations(100000);
3431 AliUnfolding::SetMinuitStrategy(2.);
3432 AliUnfolding::SetDebug(1);
3433}
3434//_____________________________________________________________________________
f3ba6c8e 3435TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) const {
4292ca60 3436 // protect heap by adding unique qualifier to name
3437 if(!protect) return 0x0;
3438 TH1D* p = (TH1D*)protect->Clone();
ef12d5a5 3439 TString tempString(fActiveString);
3440 tempString+=suffix;
3441 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
3442 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 3443 if(kill) delete protect;
3444 return p;
3445}
3446//_____________________________________________________________________________
f3ba6c8e 3447TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) const {
4292ca60 3448 // protect heap by adding unique qualifier to name
3449 if(!protect) return 0x0;
3450 TH2D* p = (TH2D*)protect->Clone();
ef12d5a5 3451 TString tempString(fActiveString);
3452 tempString+=suffix;
3453 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
3454 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 3455 if(kill) delete protect;
3456 return p;
3457}
3458//_____________________________________________________________________________
f3ba6c8e 3459TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) const {
4292ca60 3460 // protect heap by adding unique qualifier to name
3461 if(!protect) return 0x0;
3462 TGraphErrors* p = (TGraphErrors*)protect->Clone();
ef12d5a5 3463 TString tempString(fActiveString);
3464 tempString+=suffix;
3465 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
3466 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 3467 if(kill) delete protect;
3468 return p;
3469}
3470//_____________________________________________________________________________
486fb24e 3471void AliJetFlowTools::MakeAU() {
3472 // === azimuthal unfolding ===
3473 //
3474 // unfolds the spectrum in delta phi bins, extracts the yield per bin, and does a fit
3475 // in transverse momentum and azimuthal correlation space to extract v2 params
3476 // settings are equal to the ones used for 'Make()'
3477 //
3478 // basic steps that are followed:
3479 // 1) rebin the raw output of the jet task to the desired binnings
3480 // 2) calls the unfolding routine
3481 // 3) writes output to file
3482 // can be repeated multiple times with different configurations
3483
3484 Int_t low[] = {1, 6, 11, 16, 21, 26, 31, 36};
3485 Int_t up[] = {5, 10, 15, 20, 25, 30, 35, 40};
3486 TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
3487 TH1D* dPtdPhi[6];
3488 for(Int_t i(0); i < 6; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), 8, 0, TMath::Pi());
3489
3490 for(Int_t i(0); i < 8; i++) {
3491 // 1) manipulation of input histograms
3492 // check if the input variables are present
3493 if(!PrepareForUnfolding(low[i], up[i])) return;
3494 // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
3495 // parts of the spectrum can end up in over or underflow bins
3496 TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, Form("resized_%s", stringArray[i].Data()), kFALSE);
3497
3498 // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
3499 // the template will be used as a prior for the chi2 unfolding
3500 TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, stringArray[i], kFALSE);
3501
3502 // get the full response matrix from the dpt and the detector response
3503 fDetectorResponse = NormalizeTH2D(fDetectorResponse);
3504 // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
3505 // so that unfolding should return the initial spectrum
3506 if(fUseDptResponse && fUseDetectorResponse) fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
3507 else if (fUseDptResponse && !fUseDetectorResponse) fFullResponseIn = fDptIn;
3508 else if (!fUseDptResponse && fUseDetectorResponse) fFullResponseIn = fDetectorResponse;
3509 else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) return;
3510 // normalize each slide of the response to one
3511 NormalizeTH2D(fFullResponseIn);
3512 // resize to desired binning scheme
3513 TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, stringArray[i]);
3514 // get the kinematic efficiency
3515 TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
3516 kinematicEfficiencyIn->SetNameTitle(Form("kin_eff_%s", stringArray[i].Data()), Form("kin_eff_%s", stringArray[i].Data()));
3517 // suppress the errors
3518 for(Int_t j(0); j < kinematicEfficiencyIn->GetXaxis()->GetNbins(); j++) kinematicEfficiencyIn->SetBinError(1+j, 0.);
3519 TH1D* jetFindingEfficiency(0x0);
3520 if(fJetFindingEff) {
3521 jetFindingEfficiency = ProtectHeap(fJetFindingEff);
3522 jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
3523 jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
3524 }
3525 // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
3526 TH1D* unfoldedJetSpectrumIn(0x0);
3527 fActiveDir->cd(); // select active dir
3528 TDirectoryFile* dirIn = new TDirectoryFile(Form("%s___%s", stringArray[i].Data(), fActiveString.Data()), Form("%s___%s", stringArray[i].Data(), fActiveString.Data()));
3529 dirIn->cd(); // select inplane subdir
3530 // select the unfolding method
3531 unfoldedJetSpectrumIn = UnfoldWrapper(
3532 measuredJetSpectrumIn,
3533 resizedResponseIn,
3534 kinematicEfficiencyIn,
3535 measuredJetSpectrumTrueBinsIn,
3536 TString("dPtdPhiUnfolding"),
3537 jetFindingEfficiency);
3538 if(i==5) {
3539 resizedResponseIn->SetNameTitle(Form("ResponseMatrix_%s", stringArray[i].Data()), Form("response matrix %s", stringArray[i].Data()));
f3ba6c8e 3540 resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
3541 resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
486fb24e 3542 resizedResponseIn = ProtectHeap(resizedResponseIn);
3543 resizedResponseIn->Write();
3544 kinematicEfficiencyIn->SetNameTitle(Form("KinematicEfficiency_%s", stringArray[i].Data()), Form("Kinematic efficiency, %s", stringArray[i].Data()));
3545 kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
3546 kinematicEfficiencyIn->Write();
3547 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
3548 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
3549 fDetectorResponse->Write();
3550 // optional histograms
3551 if(fSaveFull) {
3552 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", Form("[INPUT] Jet spectrum, %s", stringArray[i].Data()));
3553 fSpectrumIn->Write();
3554 fDptInDist->SetNameTitle("[ORIG]DeltaPt", Form("#delta p_{T} distribution, %s", stringArray[i].Data()));
3555 fDptInDist->Write();
3556 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix", Form("#delta p_{T} matrix, %s", stringArray[i].Data()));
3557 fDptIn->Write();
3558 fFullResponseIn->SetNameTitle("ResponseMatrix", Form("Response matrix, %s", stringArray[i].Data()));
3559 fFullResponseIn->Write();
3560 }
3561 }
3562 fActiveDir->cd();
3563 fDeltaPtDeltaPhi->Write();
3564 fJetPtDeltaPhi->Write();
3565
3566 TH1D* dud(ProtectHeap(unfoldedJetSpectrumIn, kTRUE, stringArray[i]));;
3567 Double_t integralError(0);
3568 for(Int_t j(0); j < 6; j++) {
3569 // get the integrated
a39e4b2b 3570 Double_t integral(dud->IntegralAndError(2*j+1, 2*j+3, integralError));
3571 dPtdPhi[j]->SetBinContent(i+1, integral);
3572 dPtdPhi[j]->SetBinError(i+1, integralError);
486fb24e 3573 }
3574 dud->Write();
3575 // save the current state of the unfolding object
3576 SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, kFALSE);
3577 }
3578 TF1* fourier = new TF1("fourier", "[0]*(1.+0.5*[1]*(TMath::Cos(2.*x)))", 0, TMath::Pi());
3579 TH1D* v2(new TH1D("v2FromFit", "v2FromFit", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3580 for(Int_t i(0); i < 6; i++) {
3581 dPtdPhi[i]->Fit(fourier, "VI");
3582 v2->SetBinContent(1+i, fourier->GetParameter(1));
3583 v2->SetBinError(1+i, fourier->GetParError(1));
3584 dPtdPhi[i]->Write();
3585 }
3586 v2->Write();
3587}
3588//_____________________________________________________________________________