]> 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
trivial fix for ignored function return type
[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"
d7ec324f 45#include "TCanvas.h"
46#include "TLegend.h"
dc1455ee 47#include "TArrayD.h"
48#include "TList.h"
49#include "TMinuit.h"
50#include "TVirtualFitter.h"
51#include "TLegend.h"
52#include "TCanvas.h"
53#include "TStyle.h"
dc1455ee 54#include "TLine.h"
ad04a83c 55#include "TMath.h"
4292ca60 56#include "TVirtualFitter.h"
ef12d5a5 57#include "TFitResultPtr.h"
4292ca60 58// aliroot includes
dc1455ee 59#include "AliUnfolding.h"
60#include "AliAnaChargedJetResponseMaker.h"
4292ca60 61// class includes
dc1455ee 62#include "AliJetFlowTools.h"
51e6bc5a 63// roo unfold includes (make sure you have these available on your system)
64#include "RooUnfold.h"
65#include "RooUnfoldResponse.h"
66#include "RooUnfoldSvd.h"
549b5f40 67#include "RooUnfoldBayes.h"
51e6bc5a 68#include "TSVDUnfold.h"
dc1455ee 69
70using namespace std;
dc1455ee 71//_____________________________________________________________________________
72AliJetFlowTools::AliJetFlowTools() :
4292ca60 73 fResponseMaker (new AliAnaChargedJetResponseMaker()),
ef12d5a5 74 fPower (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
9f892925 75 fSaveFull (kTRUE),
dc1455ee 76 fActiveString (""),
ad04a83c 77 fActiveDir (0x0),
dc1455ee 78 fInputList (0x0),
4292ca60 79 fRefreshInput (kTRUE),
dc1455ee 80 fOutputFileName ("UnfoldedSpectra.root"),
ad04a83c 81 fOutputFile (0x0),
dc1455ee 82 fCentralityBin (0),
dc1455ee 83 fDetectorResponse (0x0),
53547ff2 84 fJetFindingEff (0x0),
4292ca60 85 fBetaIn (.1),
86 fBetaOut (.1),
549b5f40
RAB		 87 fBayesianIterIn (4),
88 fBayesianIterOut (4),
89 fBayesianSmoothIn (0.),
90 fBayesianSmoothOut (0.),
ef12d5a5 91 fAvoidRoundingError (kFALSE),
51e6bc5a 92 fUnfoldingAlgorithm (kChi2),
93 fPrior (kPriorMeasured),
dc1455ee 94 fBinsTrue (0x0),
95 fBinsRec (0x0),
ef12d5a5 96 fBinsTruePrior (0x0),
97 fBinsRecPrior (0x0),
4292ca60 98 fSVDRegIn (5),
99 fSVDRegOut (5),
51e6bc5a 100 fSVDToy (kTRUE),
101 fJetRadius (0.3),
20abfcc4 102 fEventCount (-1),
549b5f40
RAB		 103 fNormalizeSpectra (kFALSE),
104 fSmoothenPrior (kFALSE),
4292ca60 105 fFitMin (60.),
549b5f40 106 fFitMax (300.),
4292ca60 107 fFitStart (75.),
549b5f40 108 fSmoothenCounts (kTRUE),
ef12d5a5 109 fTestMode (kFALSE),
4292ca60 110 fRawInputProvided (kFALSE),
ef12d5a5 111 fEventPlaneRes (.63),
112 fUseDetectorResponse(kTRUE),
549b5f40 113 fUseDptResponse (kTRUE),
ef12d5a5 114 fTrainPower (kTRUE),
486fb24e 115 fDphiUnfolding (kTRUE),
116 fDphiDptUnfolding (kFALSE),
117 fExLJDpt (kTRUE),
4292ca60 118 fRMSSpectrumIn (0x0),
119 fRMSSpectrumOut (0x0),
120 fRMSRatio (0x0),
ef12d5a5 121 fRMSV2 (0x0),
ad04a83c 122 fDeltaPtDeltaPhi (0x0),
123 fJetPtDeltaPhi (0x0),
dc1455ee 124 fSpectrumIn (0x0),
125 fSpectrumOut (0x0),
126 fDptInDist (0x0),
127 fDptOutDist (0x0),
128 fDptIn (0x0),
129 fDptOut (0x0),
130 fFullResponseIn (0x0),
549b5f40 131 fFullResponseOut (0x0) { // class constructor
486fb24e 132 // create response maker weight function (tuned to PYTHIA spectrum)
133 fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0, 200));
134 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
4292ca60 135}
dc1455ee 136//_____________________________________________________________________________
137void AliJetFlowTools::Make() {
ad04a83c 138 // core function of the class
486fb24e 139 if(fDphiDptUnfolding) {
140 // to extract the yield as function of Dphi, Dpt - experimental
141 MakeAU();
142 return;
143 }
4292ca60 144 // 1) rebin the raw output of the jet task to the desired binnings
ad04a83c 145 // 2) calls the unfolding routine
146 // 3) writes output to file
4292ca60 147 // can be repeated multiple times with different configurations
148
ad04a83c 149 // 1) manipulation of input histograms
dc1455ee 150 // check if the input variables are present
4292ca60 151 if(fRefreshInput) {
152 if(!PrepareForUnfolding()) {
153 printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
154 return;
155 }
dc1455ee 156 }
4292ca60 157 // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
158 // parts of the spectrum can end up in over or underflow bins
f3ba6c8e 159/*
160 Double_t binsTrue[] = {-50, -45, -40, -35, -30,-25, -20,-15, -10,-5, 0, 5, 10, 15, 20,25, 30,35, 40,45, 50, 55, 60,65, 70, 75,80,85, 90,95, 100};
161 TArrayD* temparray = new TArrayD(sizeof(binsTrue)/sizeof(binsTrue[0]), binsTrue);
162 TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, temparray, TString("resized_in_"), kFALSE);
163 TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, temparray, TString("resized_out_"), kFALSE);
164 */
549b5f40 165 TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
f3ba6c8e 166 TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec, TString("resized_out_"), kFALSE);
167
549b5f40 168 // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
4292ca60 169 // the template will be used as a prior for the chi2 unfolding
549b5f40
RAB		 170 TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
171 TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
4292ca60 172 // get the full response matrix from the dpt and the detector response
dc1455ee 173 fDetectorResponse = NormalizeTH2D(fDetectorResponse);
ef12d5a5 174 // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
175 // so that unfolding should return the initial spectrum
176 if(!fTestMode) {
549b5f40
RAB		 177 if(fUseDptResponse && fUseDetectorResponse) {
178 fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
179 fFullResponseOut = MatrixMultiplication(fDptOut, fDetectorResponse);
180 } else if (fUseDptResponse && !fUseDetectorResponse) {
181 fFullResponseIn = fDptIn;
182 fFullResponseOut = fDptOut;
183 } else if (!fUseDptResponse && fUseDetectorResponse) {
184 fFullResponseIn = fDetectorResponse;
185 fFullResponseOut = fDetectorResponse;
186 } else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) {
187 printf(" > No response, exiting ! < \n" );
188 return;
189 }
ef12d5a5 190 } else {
191 fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
192 fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
193 }
4292ca60 194 // normalize each slide of the response to one
dc1455ee 195 NormalizeTH2D(fFullResponseIn);
196 NormalizeTH2D(fFullResponseOut);
4292ca60 197 // resize to desired binning scheme
549b5f40
RAB		 198 TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
199 TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
4292ca60 200 // get the kinematic efficiency
549b5f40 201 TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
4292ca60 202 kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
549b5f40 203 TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
4292ca60 204 kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
205 // suppress the errors
51e6bc5a 206 for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
207 kinematicEfficiencyIn->SetBinError(1+i, 0.);
4292ca60 208 kinematicEfficiencyOut->SetBinError(1+i, 0.);
51e6bc5a 209 }
53547ff2
RAB		 210 TH1D* jetFindingEfficiency(0x0);
211 if(fJetFindingEff) {
212 jetFindingEfficiency = ProtectHeap(fJetFindingEff);
213 jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
214 jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
215 }
ad04a83c 216 // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
549b5f40
RAB		 217 TH1D* unfoldedJetSpectrumIn(0x0);
218 TH1D* unfoldedJetSpectrumOut(0x0);
ad04a83c 219 fActiveDir->cd(); // select active dir
220 TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
221 dirIn->cd(); // select inplane subdir
486fb24e 222 // do the inplane unfolding
223 unfoldedJetSpectrumIn = UnfoldWrapper(
224 measuredJetSpectrumIn,
225 resizedResponseIn,
226 kinematicEfficiencyIn,
227 measuredJetSpectrumTrueBinsIn,
228 TString("in"),
229 jetFindingEfficiency);
549b5f40 230 resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
f3ba6c8e 231 resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
232 resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
549b5f40
RAB		 233 resizedResponseIn = ProtectHeap(resizedResponseIn);
234 resizedResponseIn->Write();
4292ca60 235 kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
236 kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
ad04a83c 237 kinematicEfficiencyIn->Write();
238 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
4292ca60 239 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
ad04a83c 240 fDetectorResponse->Write();
4292ca60 241 // optional histograms
242 if(fSaveFull) {
243 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
244 fSpectrumIn->Write();
245 fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
246 fDptInDist->Write();
247 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
248 fDptIn->Write();
249 fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
250 fFullResponseIn->Write();
251 }
ad04a83c 252 fActiveDir->cd();
486fb24e 253 if(fDphiUnfolding) {
5e11c41c 254 TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
255 dirOut->cd();
486fb24e 256 // do the out of plane unfolding
257 unfoldedJetSpectrumOut = UnfoldWrapper(
5e11c41c 258 measuredJetSpectrumOut,
259 resizedResponseOut,
260 kinematicEfficiencyOut,
261 measuredJetSpectrumTrueBinsOut,
262 TString("out"),
263 jetFindingEfficiency);
5e11c41c 264 resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
f3ba6c8e 265 resizedResponseOut->SetXTitle("p_{T, jet}^{true} [GeV/c]");
266 resizedResponseOut->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
5e11c41c 267 resizedResponseOut = ProtectHeap(resizedResponseOut);
268 resizedResponseOut->Write();
269 kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
270 kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
271 kinematicEfficiencyOut->Write();
272 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
273 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
274 fDetectorResponse->Write();
275 if(jetFindingEfficiency) jetFindingEfficiency->Write();
276 // optional histograms
277 if(fSaveFull) {
278 fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
279 fSpectrumOut->Write();
280 fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
281 fDptOutDist->Write();
282 fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
283 fDptOut->Write();
284 fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
285 fFullResponseOut->Write();
ef12d5a5 286 }
5e11c41c 287
288 // write general output histograms to file
289 fActiveDir->cd();
290 if(unfoldedJetSpectrumIn && unfoldedJetSpectrumOut && unfoldedJetSpectrumIn && unfoldedJetSpectrumOut) {
291 TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out")));
292 if(ratio) {
293 ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
f3ba6c8e 294 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
5e11c41c 295 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
296 ratio = ProtectHeap(ratio);
297 ratio->Write();
298 // write histo values to RMS files if both routines converged
299 // input values are weighted by their uncertainty
300 for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
301 if(unfoldedJetSpectrumIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumIn->GetBinError(i+1), 2.));
302 if(unfoldedJetSpectrumOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), unfoldedJetSpectrumOut->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumOut->GetBinError(i+1), 2.));
303 if(unfoldedJetSpectrumOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1) / unfoldedJetSpectrumOut->GetBinContent(i+1));
304 }
305 }
306 TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2")));
307 if(v2) {
308 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
f3ba6c8e 309 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
5e11c41c 310 v2->GetYaxis()->SetTitle("v_{2}");
311 v2 = ProtectHeap(v2);
312 v2->Write();
313 }
314 } else if (unfoldedJetSpectrumOut && unfoldedJetSpectrumIn) {
315 TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
316 if(ratio) {
317 ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
f3ba6c8e 318 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
5e11c41c 319 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
320 ratio = ProtectHeap(ratio);
321 ratio->Write();
322 }
323 TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2")));
324 if(v2) {
325 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
f3ba6c8e 326 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
5e11c41c 327 v2->GetYaxis()->SetTitle("v_{2}");
328 v2 = ProtectHeap(v2);
329 v2->Write();
330 }
4292ca60 331 }
486fb24e 332 } // end of if(fDphiUnfolding)
ad04a83c 333 fDeltaPtDeltaPhi->Write();
9f892925 334 unfoldedJetSpectrumIn->Sumw2();
335 ProtectHeap(unfoldedJetSpectrumIn, kFALSE);
336 unfoldedJetSpectrumIn->Write();
337 unfoldedJetSpectrumOut->Sumw2();
338 ProtectHeap(unfoldedJetSpectrumOut, kFALSE);
339 unfoldedJetSpectrumOut->Write();
ad04a83c 340 fJetPtDeltaPhi->Write();
549b5f40
RAB		 341 // save the current state of the unfolding object
342 SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
f3ba6c8e 343 TH1D* unfoldedJetSpectrumInForSub((TH1D*)unfoldedJetSpectrumIn->Clone("forSubIn"));
344 TH1D* unfoldedJetSpectrumOutForSub((TH1D*)unfoldedJetSpectrumOut->Clone("forSubOut"));
345 unfoldedJetSpectrumInForSub->Add(unfoldedJetSpectrumOutForSub, -1.);
346 unfoldedJetSpectrumInForSub= ProtectHeap(unfoldedJetSpectrumInForSub);
347 unfoldedJetSpectrumInForSub->Write();
348
dc1455ee 349}
350//_____________________________________________________________________________
486fb24e 351TH1D* AliJetFlowTools::UnfoldWrapper(
352 const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
353 const TH2D* resizedResponse, // response matrix
354 const TH1D* kinematicEfficiency, // kinematic efficiency
355 const TH1D* measuredJetSpectrumTrueBins, // unfolding template: same binning is pt gen of response
356 const TString suffix, // suffix (in or out of plane)
357 const TH1D* jetFindingEfficiency) // jet finding efficiency
358{
359 // wrapper function to call specific unfolding routine
360 TH1D* (AliJetFlowTools::*myFunction)(const TH1D*, const TH2D*, const TH1D*, const TH1D*, const TString, const TH1D*);
361 // initialize functon pointer
87233f72 362 if(fUnfoldingAlgorithm == kChi2) myFunction = &AliJetFlowTools::UnfoldSpectrumChi2;
363 else if(fUnfoldingAlgorithm == kBayesian) myFunction = &AliJetFlowTools::UnfoldSpectrumBayesian;
364 else if(fUnfoldingAlgorithm == kBayesianAli) myFunction = &AliJetFlowTools::UnfoldSpectrumBayesianAli;
365 else if(fUnfoldingAlgorithm == kSVD) myFunction = &AliJetFlowTools::UnfoldSpectrumSVD;
366 else if(fUnfoldingAlgorithm == kNone) {
367 TH1D* clone((TH1D*)measuredJetSpectrum->Clone("clone"));
368 clone->SetNameTitle(Form("MeasuredJetSpectrum%s", suffix.Data()), Form("measuredJetSpectrum %s", suffix.Data()));
369 return clone;
486fb24e 370 }
87233f72 371 else return 0x0;
486fb24e 372 // do the actual unfolding with the selected function
373 return (this->*myFunction)( measuredJetSpectrum, resizedResponse, kinematicEfficiency, measuredJetSpectrumTrueBins, suffix, jetFindingEfficiency);
374}
375//_____________________________________________________________________________
549b5f40
RAB		 376TH1D* AliJetFlowTools::UnfoldSpectrumChi2(
377 const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
486fb24e 378 const TH2D* resizedResponse, // response matrix
549b5f40
RAB		 379 const TH1D* kinematicEfficiency, // kinematic efficiency
380 const TH1D* measuredJetSpectrumTrueBins, // unfolding template: same binning is pt gen of response
381 const TString suffix, // suffix (in or out of plane)
382 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
dc1455ee 383{
51e6bc5a 384 // unfold the spectrum using chi2 minimization
385
4292ca60 386 // step 0) setup the static members of AliUnfolding
387 ResetAliUnfolding(); // reset from previous iteration
388 // also deletes and re-creates the global TVirtualFitter
389 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
390 if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
391 else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
392 if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
393 else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
394 AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
dc1455ee 395
549b5f40
RAB		 396 // step 1) clone all input histograms. the histograms are cloned to make sure that the original histograms
397 // stay intact. a local copy of a histogram (which only exists in the scope of this function) is
398 // denoted by the suffix 'Local'
4292ca60 399
549b5f40
RAB		 400 // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
401 TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
4292ca60 402 // unfolded local will be filled with the result of the unfolding
403 TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
404
dc1455ee 405 // full response matrix and kinematic efficiency
549b5f40 406 TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
51e6bc5a 407 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
d7ec324f 408
4292ca60 409 // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
549b5f40
RAB		 410 TH1D *priorLocal = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("priorLocal_%s", suffix.Data()));
411 // optionally, the prior can be smoothened by extrapolating the spectrum using a power law fit
412 if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
ef12d5a5 413
4292ca60 414 // step 2) start the unfolding
51e6bc5a 415 Int_t status(-1), i(0);
416 while(status < 0 && i < 100) {
4292ca60 417 // i > 0 means that the first iteration didn't converge. in that case, the result of the first
418 // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
419 if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
420 status = AliUnfolding::Unfold(
421 resizedResponseLocal, // response matrix
422 kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
549b5f40 423 measuredJetSpectrumLocal, // measured spectrum
4292ca60 424 priorLocal, // initial conditions (set NULL to use measured spectrum)
425 unfoldedLocal); // results
426 // status holds the minuit fit status (where 0 means convergence)
51e6bc5a 427 i++;
428 }
4292ca60 429 // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
430 if(status == 0 && gMinuit->fISW[1] == 3) {
431 // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
ad04a83c 432 TVirtualFitter *fitter(TVirtualFitter::GetFitter());
4292ca60 433 if(gMinuit) gMinuit->Command("SET COV");
434 TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
435 TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
436 pearson->Print();
ad04a83c 437 TH2D *hPearson(new TH2D(*pearson));
4292ca60 438 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
439 hPearson = ProtectHeap(hPearson);
ad04a83c 440 hPearson->Write();
4292ca60 441 } else status = -1;
d7ec324f 442
4292ca60 443 // step 3) refold the unfolded spectrum and save the ratio measured / refolded
444 TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
445 foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
446 unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
549b5f40 447 TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
4292ca60 448 if(ratio) {
53547ff2
RAB		 449 ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
450 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
4292ca60 451 ratio = ProtectHeap(ratio);
452 ratio->Write();
453 }
d7ec324f 454
455 // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
456 // objects are cloned using 'ProtectHeap()'
549b5f40
RAB		 457 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
458 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
459 measuredJetSpectrumLocal->Write();
d7ec324f 460
4292ca60 461 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
462 resizedResponseLocal->Write();
d7ec324f 463
4292ca60 464 unfoldedLocal = ProtectHeap(unfoldedLocal);
53547ff2 465 if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
51e6bc5a 466 unfoldedLocal->Write();
d7ec324f 467
4292ca60 468 foldedLocal = ProtectHeap(foldedLocal);
469 foldedLocal->Write();
d7ec324f 470
ef12d5a5 471 priorLocal = ProtectHeap(priorLocal);
472 priorLocal->Write();
d7ec324f 473
474 // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
53547ff2 475 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 476 fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
477 fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
478 fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
479 fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
480 fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
481 fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
53547ff2
RAB		 482 fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
483 fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
d7ec324f 484 fitStatus->Write();
549b5f40 485 return unfoldedLocal;
dc1455ee 486}
487//_____________________________________________________________________________
549b5f40
RAB		 488TH1D* AliJetFlowTools::UnfoldSpectrumSVD(
489 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
490 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
491 const TH1D* kinematicEfficiency, // kinematic efficiency
492 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
493 const TString suffix, // suffix (in, out)
494 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
51e6bc5a 495{
549b5f40
RAB		 496
497 TH1D* priorLocal( GetPrior(
498 measuredJetSpectrum, // jet pt in pt rec bins
499 resizedResponse, // full response matrix, normalized in slides of pt true
500 kinematicEfficiency, // kinematic efficiency
501 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
502 suffix, // suffix (in, out)
503 jetFindingEfficiency)); // jet finding efficiency (optional)
504 if(!priorLocal) {
505 printf(" > couldn't find prior ! < \n");
506 return 0x0;
507 } else printf(" 1) retrieved prior \n");
508
509 // go back to the 'root' directory of this instance of the SVD unfolding routine
20abfcc4 510 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
511
4292ca60 512 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
549b5f40
RAB		 513 // measured jets in pt rec binning
514 TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
515 // local copie of the response matrix
516 TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
517 // local copy of response matrix, all true slides normalized to 1
518 // this response matrix will eventually be used in the re-folding routine
519 TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
520 resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
ef12d5a5 521 // kinematic efficiency
51e6bc5a 522 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
ef12d5a5 523 // place holder histos
d7ec324f 524 TH1D *unfoldedLocalSVD(0x0);
525 TH1D *foldedLocalSVD(0x0);
4292ca60 526 cout << " 2) setup necessary input " << endl;
51e6bc5a 527 // 3) configure routine
528 RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
d7ec324f 529 cout << " step 3) configured routine " << endl;
530
4292ca60 531 // 4) get transpose matrices
549b5f40
RAB		 532 // a) get the transpose of the full response matrix
533 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
4292ca60 534 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
549b5f40
RAB		 535 // normalize it with the prior. this will ensure that high statistics bins will constrain the
536 // end result most strenuously than bins with limited number of counts
537 responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
538 cout << " 4) retrieved first transpose matrix " << endl;
4292ca60 539
540 // 5) get response for SVD unfolding
541 RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
4292ca60 542 cout << " 5) retrieved roo unfold response object " << endl;
549b5f40 543
4292ca60 544 // 6) actualy unfolding loop
549b5f40 545 RooUnfoldSvd unfoldSVD(&responseSVD, measuredJetSpectrumLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
20abfcc4 546 unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
549b5f40
RAB		 547 // correct the spectrum for the kinematic efficiency
548 unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
549
550 // get the pearson coefficients from the covariance matrix
20abfcc4 551 TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
552 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
4292ca60 553 if(pearson) {
549b5f40 554 TH2D* hPearson(new TH2D(*pearson));
4292ca60 555 pearson->Print();
d7ec324f 556 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
557 hPearson = ProtectHeap(hPearson);
4292ca60 558 hPearson->Write();
549b5f40 559 } else return 0x0; // return if unfolding didn't converge
4292ca60 560
561 // plot singular values and d_i vector
20abfcc4 562 TSVDUnfold* svdUnfold(unfoldSVD.Impl());
51e6bc5a 563 TH1* hSVal(svdUnfold->GetSV());
564 TH1D* hdi(svdUnfold->GetD());
4292ca60 565 hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
51e6bc5a 566 hSVal->SetXTitle("singular values");
567 hSVal->Write();
4292ca60 568 hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
51e6bc5a 569 hdi->SetXTitle("|d_{i}^{kreg}|");
570 hdi->Write();
4292ca60 571 cout << " plotted singular values and d_i vector " << endl;
572
573 // 7) refold the unfolded spectrum
549b5f40
RAB		 574 foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, resizedResponseLocalNorm, kinematicEfficiencyLocal);
575 TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalSVD, "ratio measured / re-folded", kTRUE));
53547ff2 576 ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
4292ca60 577 ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
578 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
579 ratio->Write();
580 cout << " 7) refolded the unfolded spectrum " << endl;
581
549b5f40
RAB		 582 // write the measured, unfolded and re-folded spectra to the output directory
583 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
584 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
585 measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
586 measuredJetSpectrumLocal->Write(); // input spectrum
d7ec324f 587 unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
588 unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
53547ff2 589 if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
4292ca60 590 unfoldedLocalSVD->Write(); // unfolded spectrum
d7ec324f 591 foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
592 foldedLocalSVD = ProtectHeap(foldedLocalSVD);
4292ca60 593 foldedLocalSVD->Write(); // re-folded spectrum
d7ec324f 594
549b5f40
RAB		 595 // save more general bookkeeeping histograms to the output directory
596 responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
f3ba6c8e 597 responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
598 responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
4292ca60 599 responseMatrixLocalTransposePrior->Write();
549b5f40
RAB		 600 priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
601 priorLocal->SetXTitle("p_{t} [GeV/c]");
602 priorLocal = ProtectHeap(priorLocal);
603 priorLocal->Write();
604 resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
605 resizedResponseLocalNorm->Write();
53547ff2
RAB		 606
607 // save some info
608 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));
609 fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
610 fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
611 fitStatus->Write();
612
549b5f40 613 return unfoldedLocalSVD;
51e6bc5a 614}
615//_____________________________________________________________________________
549b5f40
RAB		 616TH1D* AliJetFlowTools::UnfoldSpectrumBayesianAli(
617 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
618 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
619 const TH1D* kinematicEfficiency, // kinematic efficiency
620 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
621 const TString suffix, // suffix (in, out)
622 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
d7ec324f 623{
549b5f40
RAB		 624 // unfold the spectrum using the bayesian unfolding impelmented in AliUnfolding
625 // FIXME careful, not tested yet ! (06122013) FIXME
626
627 // step 0) setup the static members of AliUnfolding
628 ResetAliUnfolding(); // reset from previous iteration
629 // also deletes and re-creates the global TVirtualFitter
630 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
631 if(!strcmp("in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
632 else if(!strcmp("out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
633 else if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
634 else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
4e4f12b6
RAB		 635 AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
636
549b5f40
RAB		 637 // 1) get a prior for unfolding and clone all the input histograms
638 TH1D* priorLocal( GetPrior(
639 measuredJetSpectrum, // jet pt in pt rec bins
640 resizedResponse, // full response matrix, normalized in slides of pt true
641 kinematicEfficiency, // kinematic efficiency
642 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
643 suffix, // suffix (in, out)
644 jetFindingEfficiency)); // jet finding efficiency (optional)
645 if(!priorLocal) {
646 printf(" > couldn't find prior ! < \n");
647 return 0x0;
648 } else printf(" 1) retrieved prior \n");
649 // switch back to root dir of this unfolding procedure
650 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
651
652 // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
653 TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
654 // unfolded local will be filled with the result of the unfolding
655 TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
656
657 // full response matrix and kinematic efficiency
658 TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
659 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
660
661 // step 2) start the unfolding
662 Int_t status(-1), i(0);
663 while(status < 0 && i < 100) {
664 // i > 0 means that the first iteration didn't converge. in that case, the result of the first
665 // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
666 if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
667 status = AliUnfolding::Unfold(
668 resizedResponseLocal, // response matrix
669 kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
670 measuredJetSpectrumLocal, // measured spectrum
671 priorLocal, // initial conditions (set NULL to use measured spectrum)
672 unfoldedLocal); // results
673 // status holds the minuit fit status (where 0 means convergence)
674 i++;
675 }
676 // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
677 if(status == 0 && gMinuit->fISW[1] == 3) {
678 // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
679 TVirtualFitter *fitter(TVirtualFitter::GetFitter());
680 if(gMinuit) gMinuit->Command("SET COV");
681 TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
682 TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
683 pearson->Print();
684 TH2D *hPearson(new TH2D(*pearson));
685 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
686 hPearson = ProtectHeap(hPearson);
687 hPearson->Write();
688 } else status = -1;
689
690 // step 3) refold the unfolded spectrum and save the ratio measured / refolded
691 TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
692 foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
693 unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
694 TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
695 if(ratio) {
696 ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
697 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
698 ratio = ProtectHeap(ratio);
699 ratio->Write();
d7ec324f 700 }
549b5f40
RAB		 701
702 // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
703 // objects are cloned using 'ProtectHeap()'
704 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
705 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
706 measuredJetSpectrumLocal->Write();
707
708 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
709 resizedResponseLocal->Write();
710
711 unfoldedLocal = ProtectHeap(unfoldedLocal);
712 if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
713 unfoldedLocal->Write();
714
715 foldedLocal = ProtectHeap(foldedLocal);
716 foldedLocal->Write();
717
718 priorLocal = ProtectHeap(priorLocal);
719 priorLocal->Write();
720
721 // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
722 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));
723 fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
724 fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
725 fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
726 fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
727 fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
728 fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
729 fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
730 fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
731 fitStatus->Write();
732 return unfoldedLocal;
733}
734//_____________________________________________________________________________
735TH1D* AliJetFlowTools::UnfoldSpectrumBayesian(
736 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
737 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
738 const TH1D* kinematicEfficiency, // kinematic efficiency
739 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
740 const TString suffix, // suffix (in, out)
741 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
742{
743 // use bayesian unfolding from the RooUnfold package to unfold jet spectra
744
745 // 1) get a prior for unfolding.
746 TH1D* priorLocal( GetPrior(
747 measuredJetSpectrum, // jet pt in pt rec bins
748 resizedResponse, // full response matrix, normalized in slides of pt true
749 kinematicEfficiency, // kinematic efficiency
750 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
751 suffix, // suffix (in, out)
752 jetFindingEfficiency)); // jet finding efficiency (optional)
753 if(!priorLocal) {
754 printf(" > couldn't find prior ! < \n");
755 return 0x0;
756 } else printf(" 1) retrieved prior \n");
d7ec324f 757 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
d7ec324f 758
759 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
549b5f40
RAB		 760 // measured jets in pt rec binning
761 TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
762 // local copie of the response matrix
763 TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
764 // local copy of response matrix, all true slides normalized to 1
765 // this response matrix will eventually be used in the re-folding routine
766 TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
767 resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
d7ec324f 768 // kinematic efficiency
769 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
770 // place holder histos
549b5f40
RAB		 771 TH1D *unfoldedLocalBayes(0x0);
772 TH1D *foldedLocalBayes(0x0);
d7ec324f 773 cout << " 2) setup necessary input " << endl;
549b5f40
RAB		 774 // 4) get transpose matrices
775 // a) get the transpose of the full response matrix
776 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
d7ec324f 777 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
549b5f40
RAB		 778 // normalize it with the prior. this will ensure that high statistics bins will constrain the
779 // end result most strenuously than bins with limited number of counts
780 responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
781 // 3) get response for Bayesian unfolding
782 RooUnfoldResponse responseBayes(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedBayes_%s", suffix.Data()), Form("respCombinedBayes_%s", suffix.Data()));
d7ec324f 783
549b5f40
RAB		 784 // 4) actualy unfolding loop
785 RooUnfoldBayes unfoldBayes(&responseBayes, measuredJetSpectrumLocal, (!strcmp("in", suffix.Data())) ? fBayesianIterIn : fBayesianIterOut);
786 RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
787 unfoldedLocalBayes = (TH1D*)unfoldBayes.Hreco(errorTreatment);
788 // correct the spectrum for the kinematic efficiency
789 unfoldedLocalBayes->Divide(kinematicEfficiencyLocal);
790 // get the pearson coefficients from the covariance matrix
791 TMatrixD covarianceMatrix = unfoldBayes.Ereco(errorTreatment);
d7ec324f 792 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
d7ec324f 793 if(pearson) {
549b5f40 794 TH2D* hPearson(new TH2D(*pearson));
d7ec324f 795 pearson->Print();
796 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
797 hPearson = ProtectHeap(hPearson);
798 hPearson->Write();
549b5f40 799 } else return 0x0; // return if unfolding didn't converge
d7ec324f 800
549b5f40
RAB		 801 // 5) refold the unfolded spectrum
802 foldedLocalBayes = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalBayes, resizedResponseLocalNorm, kinematicEfficiencyLocal);
803 TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalBayes, "ratio measured / re-folded", kTRUE));
804 ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
d7ec324f 805 ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
806 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
807 ratio->Write();
808 cout << " 7) refolded the unfolded spectrum " << endl;
809
549b5f40
RAB		 810 // write the measured, unfolded and re-folded spectra to the output directory
811 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
812 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
813 measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
814 measuredJetSpectrumLocal->Write(); // input spectrum
815 unfoldedLocalBayes->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
816 unfoldedLocalBayes = ProtectHeap(unfoldedLocalBayes);
817 if(jetFindingEfficiency) unfoldedLocalBayes->Divide(jetFindingEfficiency);
818 unfoldedLocalBayes->Write(); // unfolded spectrum
819 foldedLocalBayes->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
820 foldedLocalBayes = ProtectHeap(foldedLocalBayes);
821 foldedLocalBayes->Write(); // re-folded spectrum
822
823 // save more general bookkeeeping histograms to the output directory
824 responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
f3ba6c8e 825 responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
826 responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
d7ec324f 827 responseMatrixLocalTransposePrior->Write();
549b5f40
RAB		 828 priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
829 priorLocal->SetXTitle("p_{t} [GeV/c]");
830 priorLocal = ProtectHeap(priorLocal);
831 priorLocal->Write();
832 resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
833 resizedResponseLocalNorm->Write();
53547ff2
RAB		 834
835 // save some info
836 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		 837 fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fBayesianIterIn : fBayesianIterOut);
838 fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fBayesianIterIn" : "fBayesianIterOut");
53547ff2
RAB		 839 fitStatus->Write();
840
549b5f40 841 return unfoldedLocalBayes;
d7ec324f 842}
843//_____________________________________________________________________________
dc1455ee 844Bool_t AliJetFlowTools::PrepareForUnfolding()
845{
846 // prepare for unfolding
4292ca60 847 if(fRawInputProvided) return kTRUE;
dc1455ee 848 if(!fInputList) {
849 printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
850 return kFALSE;
851 }
549b5f40 852 if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
4292ca60 853 // check if the pt bin for true and rec have been set
854 if(!fBinsTrue || !fBinsRec) {
855 printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
856 return kFALSE;
dc1455ee 857 }
486fb24e 858 if(!fRMSSpectrumIn && fDphiUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
4292ca60 859 // procedures, these profiles will be nonsensical, user is responsible
860 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
861 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
862 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
dc1455ee 863 }
ef12d5a5 864 if(!fTrainPower) {
549b5f40 865 // clear minuit state to avoid constraining the fit with the results of the previous iteration
ef12d5a5 866 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
867 }
dc1455ee 868 // extract the spectra
869 TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin));
ad04a83c 870 fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
871 if(!fJetPtDeltaPhi) {
dc1455ee 872 printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
873 return kFALSE;
874 }
4292ca60 875 fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
dc1455ee 876 // in plane spectrum
486fb24e 877 if(!fDphiUnfolding) {
5e11c41c 878 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e");
879 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40, "e");
ef12d5a5 880 } else {
5e11c41c 881 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10, "e");
882 fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40, "e"));
ef12d5a5 883 fSpectrumIn = ProtectHeap(fSpectrumIn);
884 // out of plane spectrum
5e11c41c 885 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30, "e");
ef12d5a5 886 fSpectrumOut = ProtectHeap(fSpectrumOut);
887 }
20abfcc4 888 // normalize spectra to event count if requested
889 if(fNormalizeSpectra) {
890 TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityBin)));
ef12d5a5 891 if(!rho) return 0x0;
892 Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
893 if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
20abfcc4 894 if(fEventCount > 0) {
4292ca60 895 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
896 fSpectrumOut->Sumw2();
ef12d5a5 897 Double_t pt(0);
d7ec324f 898 Double_t error(0); // lots of issues with the errors here ...
ef12d5a5 899 for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
4292ca60 900 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount; // normalized count
d7ec324f 901 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
4292ca60 902 fSpectrumIn->SetBinContent(1+i, pt);
ef12d5a5 903 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
d7ec324f 904 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
ef12d5a5 905 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
4292ca60 906 }
ef12d5a5 907 for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
908 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount; // normalized count
d7ec324f 909 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
ef12d5a5 910 fSpectrumOut->SetBinContent(1+i, pt);
911 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
d7ec324f 912 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
ef12d5a5 913 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
4292ca60 914 }
20abfcc4 915 }
ef12d5a5 916 if(normalizeToFullSpectrum) fEventCount = -1;
20abfcc4 917 }
dc1455ee 918 // extract the delta pt matrices
486fb24e 919 TString deltaptName("");
920 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityBin);
ad04a83c 921 fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
922 if(!fDeltaPtDeltaPhi) {
dc1455ee 923 printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
549b5f40
RAB		 924 printf(" > may be ok, depending no what you want to do < \n");
925 fRefreshInput = kTRUE;
926 return kTRUE;
dc1455ee 927 }
4292ca60 928 fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
929 // in plane delta pt distribution
486fb24e 930 if(!fDphiUnfolding) {
5e11c41c 931 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e");
932 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40, "e");
ef12d5a5 933 } else {
5e11c41c 934 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10, "e");
935 fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40, "e"));
ef12d5a5 936 // out of plane delta pt distribution
5e11c41c 937 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30, "e");
ef12d5a5 938 fDptInDist = ProtectHeap(fDptInDist);
939 fDptOutDist = ProtectHeap(fDptOutDist);
940 // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
941 }
4292ca60 942
dc1455ee 943 // create a rec - true smeared response matrix
944 TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
945 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
ad04a83c 946 Bool_t skip = kFALSE;
dc1455ee 947 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
ad04a83c 948 (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
949 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
dc1455ee 950 }
951 }
952 TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
953 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
ad04a83c 954 Bool_t skip = kFALSE;
dc1455ee 955 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
ad04a83c 956 (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
957 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
dc1455ee 958 }
959 }
960 fDptIn = new TH2D(*rfIn);
961 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
f3ba6c8e 962 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
963 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
4292ca60 964 fDptIn = ProtectHeap(fDptIn);
dc1455ee 965 fDptOut = new TH2D(*rfOut);
966 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
f3ba6c8e 967 fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
968 fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
4292ca60 969 fDptOut = ProtectHeap(fDptOut);
970
971 fRefreshInput = kTRUE; // force cloning of the input
dc1455ee 972 return kTRUE;
973}
974//_____________________________________________________________________________
486fb24e 975Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) {
976 // prepare for unfoldingUA - more robust method to extract input spectra from file
977 // will replace PrepareForUnfolding eventually (09012014)
978 if(!fInputList) {
979 printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
980 return kFALSE;
981 }
982 if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
983 // check if the pt bin for true and rec have been set
984 if(!fBinsTrue || !fBinsRec) {
985 printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
986 return kFALSE;
987 }
988 if(!fTrainPower) {
989 // clear minuit state to avoid constraining the fit with the results of the previous iteration
990 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
991 }
992 // extract the spectra
993 TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin));
994 fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
995 if(!fJetPtDeltaPhi) {
996 printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
997 return kFALSE;
998 }
999 fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
1000 // in plane spectrum
1001 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e");
1002 // extract the delta pt matrices
1003 TString deltaptName("");
1004 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityBin);
1005 fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1006 if(!fDeltaPtDeltaPhi) {
1007 printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
1008 printf(" > may be ok, depending no what you want to do < \n");
1009 fRefreshInput = kTRUE;
1010 return kTRUE;
1011 }
1012 fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
1013 // in plane delta pt distribution
1014 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), low, up, "e");
1015 // create a rec - true smeared response matrix
1016 TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1017 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
1018 Bool_t skip = kFALSE;
1019 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
1020 (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1021 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1022 }
1023 }
1024 fDptIn = new TH2D(*rfIn);
1025 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
f3ba6c8e 1026 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1027 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
486fb24e 1028 fDptIn = ProtectHeap(fDptIn);
1029
1030 return kTRUE;
1031}
1032//_____________________________________________________________________________
549b5f40
RAB		 1033TH1D* AliJetFlowTools::GetPrior(
1034 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
1035 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
1036 const TH1D* kinematicEfficiency, // kinematic efficiency
1037 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
1038 const TString suffix, // suffix (in, out)
1039 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
1040{
1041 // 1) get a prior for unfolding.
1042 // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
1043 TH1D* unfolded(0x0);
1044 TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
1045 dirOut->cd();
1046 switch (fPrior) { // select the prior for unfolding
1047 case kPriorChi2 : {
1048 if(fBinsTruePrior && fBinsRecPrior) { // if set, use different binning for the prior
1049 TArrayD* tempArrayTrue(fBinsTrue); // temporarily cache the original binning
1050 fBinsTrue = fBinsTruePrior; // switch binning schemes (will be used in UnfoldSpectrumChi2())
1051 TArrayD* tempArrayRec(fBinsRec);
1052 fBinsRec = fBinsRecPrior;
1053 TH1D* measuredJetSpectrumChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
1054 TH1D* measuredJetSpectrumTrueBinsChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
1055 TH2D* resizedResponseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
1056 TH1D* kinematicEfficiencyChi2(resizedResponseChi2->ProjectionX());
1057 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
1058 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
1059 unfolded= UnfoldSpectrumChi2(
1060 measuredJetSpectrumChi2,
1061 resizedResponseChi2,
1062 kinematicEfficiencyChi2,
1063 measuredJetSpectrumTrueBinsChi2, // prior for chi2 unfolding (measured)
1064 TString(Form("prior_%s", suffix.Data())));
1065 if(!unfolded) {
1066 printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1067 printf(" probably Chi2 unfolding did not converge < \n");
1068 return 0x0;
1069 }
1070 fBinsTrue = tempArrayTrue; // reset bins borders
1071 fBinsRec = tempArrayRec;
1072 // if the chi2 prior has a different binning, rebin to the true binning for the unfolding
1073 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data()))); // rebin unfolded
1074 } else {
1075 unfolded = UnfoldSpectrumChi2(
1076 measuredJetSpectrum,
1077 resizedResponse,
1078 kinematicEfficiency,
1079 measuredJetSpectrumTrueBins, // prior for chi2 unfolding (measured)
1080 TString(Form("prior_%s", suffix.Data())));
1081 if(!unfolded) {
1082 printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1083 printf(" probably Chi2 unfolding did not converge < \n");
1084 return 0x0;
1085 }
1086 }
1087 break;
1088 }
1089 case kPriorMeasured : {
1090 unfolded = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("kPriorMeasured_%s", suffix.Data())); // copy template to unfolded to use as prior
1091 }
1092 default : break;
1093 }
1094 // it can be important that the prior is smooth, this can be achieved by
1095 // extrapolating the spectrum with a fitted power law when bins are sparsely filed
1096 if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
1097 TH1D *priorLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
1098 if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
1099 return priorLocal;
1100}
1101//_____________________________________________________________________________
dc1455ee 1102TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
1103 // resize the x-axis of a th1d
1104 if(!histo) {
1105 printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1106 return NULL;
1107 }
1108 // see how many bins we need to copy
1109 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);
1110 // low is the bin number of the first new bin
1111 Int_t l = histo->GetXaxis()->FindBin(low);
1112 // set the values
1113 Double_t _x(0), _xx(0);
1114 for(Int_t i(0); i < up-low; i++) {
1115 _x = histo->GetBinContent(l+i);
1116 _xx=histo->GetBinError(l+i);
1117 resized->SetBinContent(i+1, _x);
1118 resized->SetBinError(i+1, _xx);
1119 }
1120 return resized;
1121}
1122//_____________________________________________________________________________
1123TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
1124 // resize the y-axis of a th2d
1125 if(!histo) {
1126 printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1127 return NULL;
1128 }
1129 // see how many bins we need to copy
1130 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());
1131 // assume only the y-axis has changed
1132 // low is the bin number of the first new bin
1133 Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1134 // set the values
1135 Double_t _x(0), _xx(0);
1136 for(Int_t i(0); i < x->GetSize(); i++) {
1137 for(Int_t j(0); j < y->GetSize(); j++) {
1138 _x = histo->GetBinContent(i, low+j);
1139 _xx=histo->GetBinError(i, low+1+j);
1140 resized->SetBinContent(i, j, _x);
1141 resized->SetBinError(i, j, _xx);
1142 }
1143 }
1144 return resized;
1145}
1146//_____________________________________________________________________________
512ced40 1147TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo, Bool_t noError) {
dc1455ee 1148 // general method to normalize all vertical slices of a th2 to unity
1149 // i.e. get a probability matrix
1150 if(!histo) {
1151 printf(" > NormalizeTH2D:: fatal error, NULL pointer passed < \n");
1152 return NULL;
1153 }
1154 Int_t binsX = histo->GetXaxis()->GetNbins();
1155 Int_t binsY = histo->GetYaxis()->GetNbins();
1156
1157 // normalize all slices in x
1158 for(Int_t i(0); i < binsX; i++) { // for each vertical slice
1159 Double_t weight = 0;
1160 for(Int_t j(0); j < binsY; j++) { // loop over all the horizontal components
1161 weight+=histo->GetBinContent(i+1, j+1);
1162 } // now we know the total weight
1163 for(Int_t j(0); j < binsY; j++) {
1164 if (weight <= 0 ) continue;
1165 histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
512ced40
RAB		 1166 if(noError) histo->SetBinError( 1+i, j+1, 0.);
1167 else histo->SetBinError( 1+i, j+1, histo->GetBinError( 1+i, j+1)/weight);
dc1455ee 1168 }
1169 }
1170 return histo;
1171}
1172//_____________________________________________________________________________
53547ff2 1173TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
dc1455ee 1174 // return a TH1D with the supplied histogram rebinned to the supplied bins
53547ff2 1175 // the returned histogram is new, the original is deleted from the heap if kill is true
dc1455ee 1176 if(!histo || !bins) {
53547ff2 1177 printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
dc1455ee 1178 return NULL;
1179 }
1180 // create the output histo
1181 TString name = histo->GetName();
1182 name+="_template";
1183 name+=suffix;
1184 TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1185 for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
4292ca60 1186 // loop over the bins of the old histo and fill the new one with its data
1187 rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
dc1455ee 1188 }
53547ff2 1189 if(kill) delete histo;
dc1455ee 1190 return rebinned;
1191}
1192//_____________________________________________________________________________
4292ca60 1193TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
53547ff2 1194 // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
4292ca60 1195 if(!fResponseMaker || !binsTrue || !binsRec) {
1196 printf(" > RebinTH2D:: function called with NULL arguments < \n");
1197 return 0x0;
dc1455ee 1198 }
4292ca60 1199 TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1200 return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
dc1455ee 1201}
1202//_____________________________________________________________________________
d7ec324f 1203TH2D* AliJetFlowTools::MatrixMultiplication(TH2D* a, TH2D* b, TString name)
1204{
1205 // multiply two matrices
1206 if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1207 TH2D* c = (TH2D*)a->Clone("c");
1208 for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1209 for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1210 Double_t val = 0;
1211 for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1212 Int_t y2 = x1;
1213 val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1214 }
1215 c->SetBinContent(x2, y1, val);
512ced40 1216 c->SetBinError(x2, y1, 0.);
dc1455ee 1217 }
1218 }
d7ec324f 1219 if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1220 return c;
dc1455ee 1221}
1222//_____________________________________________________________________________
d7ec324f 1223TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale)
1224{
dc1455ee 1225 // normalize a th1d to a certain scale
4292ca60 1226 histo->Sumw2();
1227 Double_t integral = histo->Integral()*scale;
1228 if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1229 else if (scale != 1.) histo->Scale(1./scale, "width");
1230 else printf(" > Histogram integral < 0, cannot normalize \n");
1231 return histo;
dc1455ee 1232}
1233//_____________________________________________________________________________
51e6bc5a 1234TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
dc1455ee 1235{
1236 // Calculate pearson coefficients from covariance matrix
51e6bc5a 1237 TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1238 Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
ad04a83c 1239 Double_t pearson(0.);
dc1455ee 1240 if(nrows==0 && ncols==0) return 0x0;
ef12d5a5 1241 for(Int_t row = 0; row < nrows; row++) {
1242 for(Int_t col = 0; col<ncols; col++) {
51e6bc5a 1243 if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1244 (*pearsonCoefficients)(row,col) = pearson;
dc1455ee 1245 }
1246 }
51e6bc5a 1247 return pearsonCoefficients;
dc1455ee 1248}
1249//_____________________________________________________________________________
549b5f40 1250TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
d7ec324f 1251 // smoothen the spectrum using a user defined function
1252 // returns a clone of the original spectrum if fitting failed
1253 // if kill is kTRUE the input spectrum will be deleted from the heap
1254 // if 'count' is selected, bins are filled with integers (necessary if the
1255 // histogram is interpreted in a routine which accepts only counts)
549b5f40
RAB		 1256 if(!spectrum || !function) return 0x0;
1257 if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
d7ec324f 1258 TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1259 temp->Sumw2(); // if already called on the original, this will give off a warning but do nothing
87233f72 1260 TFitResultPtr r = temp->Fit(function, "", "", min, max);
d7ec324f 1261 if((int)r == 0) { // MINUIT status
1262 for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1263 if(temp->GetBinCenter(i) > start) { // from this pt value use extrapolation
486fb24e 1264 (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 1265 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1266 }
1267 }
1268 }
1269 if(kill) delete spectrum;
1270 return temp;
1271}
1272//_____________________________________________________________________________
1273void AliJetFlowTools::Style(TCanvas* c, TString style)
1274{
1275 // set a default style for a canvas
1276 if(!strcmp(style.Data(), "PEARSON")) {
1277 printf(" > style PEARSON canvas < \n");
1278 gStyle->SetOptStat(0);
1279 c->SetGridx();
1280 c->SetGridy();
1281 c->SetTicks();
1282 return;
1283 } else if(!strcmp(style.Data(), "SPECTRUM")) {
1284 printf(" > style SPECTRUM canvas < \n");
1285 gStyle->SetOptStat(0);
1286 c->SetLogy();
1287 c->SetGridx();
1288 c->SetGridy();
1289 c->SetTicks();
1290 return;
1291 } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1292}
1293//_____________________________________________________________________________
1294void AliJetFlowTools::Style(TVirtualPad* c, TString style)
1295{
1296 // set a default style for a canvas
1297 if(!strcmp(style.Data(), "PEARSON")) {
1298 printf(" > style PEARSON pad < \n");
1299 gStyle->SetOptStat(0);
1300 c->SetGridx();
1301 c->SetGridy();
1302 c->SetTicks();
1303 return;
1304 } else if(!strcmp(style.Data(), "SPECTRUM")) {
1305 printf(" > style SPECTRUM pad < \n");
1306 gStyle->SetOptStat(0);
1307 c->SetLogy();
1308 c->SetGridx();
1309 c->SetGridy();
1310 c->SetTicks();
1311 return;
1312 } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1313}
1314//_____________________________________________________________________________
53547ff2
RAB		 1315void AliJetFlowTools::Style(TLegend* l)
1316{
1317 // set a default style for a legend
1318 l->SetTextSize(.06);
1319 l->SetFillColor(0);
1320}
1321//_____________________________________________________________________________
1322void AliJetFlowTools::Style(TH1* h, EColor col, histoType type)
1323{
1324 // style a histo
1325 h->SetLineColor(col);
1326 h->SetMarkerColor(col);
1327 h->SetLineWidth(2.);
1328 h->SetMarkerSize(2.);
1329 h->SetTitleSize(0.06);
1330 h->GetYaxis()->SetLabelSize(0.06);
1331 h->GetXaxis()->SetLabelSize(0.06);
1332 h->GetYaxis()->SetTitleSize(0.06);
1333 h->GetXaxis()->SetTitleSize(0.06);
1334 h->GetYaxis()->SetTitleOffset(1.);
1335 h->GetXaxis()->SetTitleOffset(.9);
1336 switch (type) {
1337 case kInPlaneSpectrum : {
1338 h->SetTitle("IN PLANE");
f3ba6c8e 1339 h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
53547ff2
RAB		 1340 h->GetYaxis()->SetTitle("counts");
1341 } break;
1342 case kOutPlaneSpectrum : {
1343 h->SetTitle("OUT OF PLANE");
f3ba6c8e 1344 h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
53547ff2
RAB		 1345 h->GetYaxis()->SetTitle("counts");
1346 } break;
1347 case kUnfoldedSpectrum : {
1348 h->SetTitle("UNFOLDED");
f3ba6c8e 1349 h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
53547ff2
RAB		 1350 h->GetYaxis()->SetTitle("counts");
1351 } break;
1352 case kFoldedSpectrum : {
1353 h->SetTitle("FOLDED");
f3ba6c8e 1354 h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
53547ff2
RAB		 1355 h->GetYaxis()->SetTitle("counts");
1356 } break;
1357 case kMeasuredSpectrum : {
1358 h->SetTitle("MEASURED");
f3ba6c8e 1359 h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
53547ff2
RAB		 1360 h->GetYaxis()->SetTitle("counts");
1361 } break;
1362 default : break;
1363 }
1364}
1365//_____________________________________________________________________________
1366void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type)
1367{
1368 // style a histo
1369 h->SetLineColor(col);
1370 h->SetMarkerColor(col);
1371 h->SetLineWidth(2.);
1372 h->SetMarkerSize(2.);
1373 h->GetYaxis()->SetLabelSize(0.06);
1374 h->GetXaxis()->SetLabelSize(0.06);
1375 h->GetYaxis()->SetTitleSize(0.06);
1376 h->GetXaxis()->SetTitleSize(0.06);
1377 h->GetYaxis()->SetTitleOffset(1.);
1378 h->GetXaxis()->SetTitleOffset(.9);
1379 switch (type) {
1380 case kInPlaneSpectrum : {
1381 h->SetTitle("IN PLANE");
f3ba6c8e 1382 h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
53547ff2
RAB		 1383 h->GetYaxis()->SetTitle("counts");
1384 } break;
1385 case kOutPlaneSpectrum : {
1386 h->SetTitle("OUT OF PLANE");
f3ba6c8e 1387 h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
53547ff2
RAB		 1388 h->GetYaxis()->SetTitle("counts");
1389 } break;
1390 case kUnfoldedSpectrum : {
1391 h->SetTitle("UNFOLDED");
f3ba6c8e 1392 h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
53547ff2
RAB		 1393 h->GetYaxis()->SetTitle("counts");
1394 } break;
1395 case kFoldedSpectrum : {
1396 h->SetTitle("FOLDED");
f3ba6c8e 1397 h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
53547ff2
RAB		 1398 h->GetYaxis()->SetTitle("counts");
1399 } break;
1400 case kMeasuredSpectrum : {
1401 h->SetTitle("MEASURED");
f3ba6c8e 1402 h->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
53547ff2
RAB		 1403 h->GetYaxis()->SetTitle("counts");
1404 } break;
1405 default : break;
1406 }
1407}
1408//_____________________________________________________________________________
f3ba6c8e 1409void AliJetFlowTools::DoSystematics(Int_t nominalIn, Int_t nominalOut, TArrayI* variationsIn, TArrayI* variationsOut, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
1410{
1411 // do systematics. structure of this function is similar to that of PostProcess, however
1412 // in this case a nomial index can be supplied as well as an array of systematic checks
1413 //
1414 // FIXME make separate canvas for nominal value
1415 //
1416 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
1417 TFile readMe(in.Data(), "READ"); // open file read-only
1418 if(readMe.IsZombie()) {
1419 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1420 return;
1421 }
1422 printf("\n\n\n\t\t DOSYSTEMATICS \n > Recovered the following file structure : \n <");
1423 readMe.ls();
1424 TList* listOfKeys((TList*)readMe.GetListOfKeys());
1425 if(!listOfKeys) {
1426 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
1427 return;
1428 }
1429 // check input params
1430 if(variationsIn->GetSize() != variationsOut->GetSize()) {
1431 printf(" > DoSystematics: fatal error, input arrays have different sizes ! < \n ");
1432 return;
1433 }
1434 TDirectoryFile* defRootDirIn(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(nominalIn)->GetName())));
1435 TDirectoryFile* defRootDirOut(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(nominalOut)->GetName())));
1436 if(!(defRootDirIn && defRootDirOut)) {
1437 printf(" > DoSystematics: fatal error, couldn't retrieve nominal values ! < \n ");
1438 return;
1439 }
1440 TString defIn(defRootDirIn->GetName());
1441 TString defOut(defRootDirOut->GetName());
1442 // prepare necessary canvasses
1443 TCanvas* canvasIn(new TCanvas("SYSTPearsonIn", "SYSTPearsonIn"));
1444 TCanvas* canvasOut(0x0);
1445 if(fDphiUnfolding) canvasOut = new TCanvas("SYSTPearsonOut", "SYSTPearsonOut");
1446 TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("SYSTRefoldedIn", "SYSTRefoldedIn"));
1447 TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
1448 if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("SYSTRefoldedOut", "SYSTRefoldedOut");
1449 TCanvas* canvasSpectraIn(new TCanvas("SYSTSpectraIn", "SYSTSpectraIn"));
1450 TCanvas* canvasSpectraOut(0x0);
1451 if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SYSTSpectraOut", "SYSTSpectraOut");
1452 TCanvas* canvasRatio(0x0);
1453 if(fDphiUnfolding) canvasRatio = new TCanvas("SYSTRatio", "SYSTRatio");
1454 TCanvas* canvasV2(0x0);
1455 if(fDphiUnfolding) canvasV2 = new TCanvas("SYSTV2", "SYSTV2");
1456 TCanvas* canvasMISC(new TCanvas("SYSTMISC", "SYSTMISC"));
1457 TCanvas* canvasMasterIn(new TCanvas("SYSTdefaultIn", "SYSTdefaultIn"));
1458 TCanvas* canvasMasterOut(0x0);
1459 if(fDphiUnfolding) canvasMasterOut = new TCanvas("SYSTdefaultOut", "SYSTdefaultOut");
1460 (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
1461
1462 TCanvas* canvasProfiles(new TCanvas("SYSTcanvasProfiles", "SYSTcanvasProfiles"));
1463 canvasProfiles->Divide(2);
1464 TProfile* ratioProfile(new TProfile("SYSTratioProfile", "SYSTratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1465 TProfile* v2Profile(new TProfile("SYSTv2Profile", "SYSTv2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1466 // get an estimate of the number of outputs and find the default set
1467 Int_t rows(TMath::Floor(variationsIn->GetSize()/(float)columns)+((variationsIn->GetSize()%columns)>0));
1468 canvasIn->Divide(columns, rows);
1469 if(canvasOut) canvasOut->Divide(columns, rows);
1470 canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
1471 if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
1472 canvasSpectraIn->Divide(columns, rows);
1473 if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
1474 if(canvasRatio) canvasRatio->Divide(columns, rows);
1475 if(canvasV2) canvasV2->Divide(columns, rows);
1476 canvasMasterIn->Divide(columns, rows);
1477 if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
1478
1479 // prepare a separate set of canvases to hold the nominal points
1480 TCanvas* canvasNominalIn(new TCanvas("NominalPearsonIn", "NominalPearsonIn"));
1481 TCanvas* canvasNominalOut(0x0);
1482 if(fDphiUnfolding) canvasNominalOut = new TCanvas("NominalPearsonOut", "NominalPearsonOut");
1483 TCanvas* canvasNominalRatioMeasuredRefoldedIn(new TCanvas("NominalRefoldedIn", "NominalRefoldedIn"));
1484 TCanvas* canvasNominalRatioMeasuredRefoldedOut(0x0);
1485 if(fDphiUnfolding) canvasNominalRatioMeasuredRefoldedOut = new TCanvas("NominalRefoldedOut", "NominalRefoldedOut");
1486 TCanvas* canvasNominalSpectraIn(new TCanvas("NominalSpectraIn", "NominalSpectraIn"));
1487 TCanvas* canvasNominalSpectraOut(0x0);
1488 if(fDphiUnfolding) canvasNominalSpectraOut = new TCanvas("NominalSpectraOut", "NominalSpectraOut");
1489 TCanvas* canvasNominalRatio(0x0);
1490 if(fDphiUnfolding) canvasNominalRatio = new TCanvas("NominalRatio", "NominalRatio");
1491 TCanvas* canvasNominalV2(0x0);
1492 if(fDphiUnfolding) canvasNominalV2 = new TCanvas("NominalV2", "NominalV2");
1493 TCanvas* canvasNominalMISC(new TCanvas("NominalMISC", "NominalMISC"));
1494 TCanvas* canvasNominalMasterIn(new TCanvas("NominaldefaultIn", "NominaldefaultIn"));
1495 TCanvas* canvasNominalMasterOut(0x0);
1496 if(fDphiUnfolding) canvasNominalMasterOut = new TCanvas("NominaldefaultOut", "NominaldefaultOut");
1497 (fDphiUnfolding) ? canvasNominalMISC->Divide(4, 2) : canvasNominalMISC->Divide(4, 1);
1498
1499 columns = 1; rows = 1;
1500 canvasNominalIn->Divide(columns, rows);
1501 if(canvasNominalOut) canvasNominalOut->Divide(columns, rows);
1502 canvasNominalRatioMeasuredRefoldedIn->Divide(columns, rows);
1503 if(canvasNominalRatioMeasuredRefoldedOut) canvasNominalRatioMeasuredRefoldedOut->Divide(columns, rows);
1504 canvasNominalSpectraIn->Divide(columns, rows);
1505 if(canvasNominalSpectraOut) canvasNominalSpectraOut->Divide(columns, rows);
1506 if(canvasNominalRatio) canvasNominalRatio->Divide(columns, rows);
1507 if(canvasNominalV2) canvasNominalV2->Divide(columns, rows);
1508
1509 canvasNominalMasterIn->Divide(columns, rows);
1510 if(canvasNominalMasterOut) canvasNominalMasterOut->Divide(columns, rows);
1511
1512 // extract the default output
1513 TH1D* defaultUnfoldedJetSpectrumIn(0x0);
1514 TH1D* defaultUnfoldedJetSpectrumOut(0x0);
1515 TDirectoryFile* defDirIn = (TDirectoryFile*)defRootDirIn->Get(Form("InPlane___%s", defIn.Data()));
1516 TDirectoryFile* defDirOut = (TDirectoryFile*)defRootDirOut->Get(Form("OutOfPlane___%s", defOut.Data()));
1517 if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", defIn.Data()));
1518 if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", defOut.Data()));
1519 printf(" > succesfully extracted default results < \n");
1520
1521 // loop through the directories, only plot the graphs if the deconvolution converged
1522 TDirectoryFile* tempDirIn(0x0);
1523 TDirectoryFile* tempDirOut(0x0);
1524 TDirectoryFile* tempIn(0x0);
1525 TDirectoryFile* tempOut(0x0);
1526 TH1D* unfoldedSpectrumInForRatio(0x0);
1527 TH1D* unfoldedSpectrumOutForRatio(0x0);
1528 for(Int_t i(0), j(-1); i < variationsIn->GetSize(); i++) {
1529 tempDirIn = (dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(variationsIn->At(i))->GetName())));
1530 tempDirOut = (dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(variationsOut->At(i))->GetName())));
1531 if(!(tempDirIn && tempDirOut)) {
1532 printf(" > DoSystematics: couldn't get a set of variations < \n");
1533 continue;
1534 }
1535 TString dirNameIn(tempDirIn->GetName());
1536 TString dirNameOut(tempDirOut->GetName());
1537 // try to read the in- and out of plane subdirs
1538 tempIn = (TDirectoryFile*)tempDirIn->Get(Form("InPlane___%s", dirNameIn.Data()));
1539 tempOut = (TDirectoryFile*)tempDirOut->Get(Form("OutOfPlane___%s", dirNameOut.Data()));
1540 j++;
1541 if(tempIn) {
1542 // to see if the unfolding converged try to extract the pearson coefficients
1543 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirNameIn.Data())));
1544 if(pIn) {
1545 printf(" - %s in plane converged \n", dirNameIn.Data());
1546 canvasIn->cd(j);
1547 if(i==0) canvasNominalIn->cd(j);
1548 Style(gPad, "PEARSON");
1549 pIn->DrawCopy("colz");
1550 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirNameIn.Data())));
1551 if(rIn) {
1552 Style(rIn);
1553 printf(" > found RatioRefoldedMeasured < \n");
1554 canvasRatioMeasuredRefoldedIn->cd(j);
1555 if(i==0) canvasNominalRatioMeasuredRefoldedIn->cd(j);
1556 rIn->SetFillColor(kRed);
1557 rIn->Draw("ap");
1558 }
1559 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
1560 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
1561 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirNameIn.Data())));
1562 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirNameIn.Data())));
1563 if(dvector && avalue && rm && eff) {
1564 Style(dvector);
1565 Style(avalue);
1566 Style(rm);
1567 Style(eff);
1568 canvasMISC->cd(1);
1569 if(i==0) canvasNominalMISC->cd(1);
1570 Style(gPad, "SPECTRUM");
1571 dvector->DrawCopy();
1572 canvasMISC->cd(2);
1573 if(i==0) canvasNominalMISC->cd(2);
1574 Style(gPad, "SPECTRUM");
1575 avalue->DrawCopy();
1576 canvasMISC->cd(3);
1577 if(i==0) canvasNominalMISC->cd(3);
1578 Style(gPad, "PEARSON");
1579 rm->DrawCopy("colz");
1580 canvasMISC->cd(4);
1581 if(i==0) canvasNominalMISC->cd(4);
1582 eff->DrawCopy();
1583 } else if(rm && eff) {
1584 Style(rm);
1585 Style(eff);
1586 canvasMISC->cd(3);
1587 if(i==0) canvasNominalMISC->cd(3);
1588 Style(gPad, "PEARSON");
1589 rm->DrawCopy("colz");
1590 canvasMISC->cd(4);
1591 if(i==0) canvasNominalMISC->cd(4);
1592 eff->DrawCopy();
1593 }
1594 }
1595 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirNameIn.Data())));
1596 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirNameIn.Data())));
1597 unfoldedSpectrumInForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
1598 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirNameIn.Data())));
1599 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1600 if(defaultUnfoldedJetSpectrumIn) {
1601 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
1602 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirNameIn.Data())));
1603 temp->Divide(unfoldedSpectrum);
1604 temp->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
1605 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1606 temp->GetYaxis()->SetTitle("ratio");
1607 canvasMasterIn->cd(j);
1608 if(i==0) canvasNominalMasterIn->cd(j);
1609 temp->GetYaxis()->SetRangeUser(0., 2);
1610 temp->DrawCopy();
1611 }
1612 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirNameIn.Data())));
1613 canvasSpectraIn->cd(j);
1614 if(i==0) canvasNominalSpectraIn->cd(j);
1615 Style(gPad);
1616 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
1617 unfoldedSpectrum->DrawCopy();
1618 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
1619 inputSpectrum->DrawCopy("same");
1620 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
1621 refoldedSpectrum->DrawCopy("same");
1622 TLegend* l(AddLegend(gPad));
1623 Style(l);
1624 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1625 Float_t chi(fitStatus->GetBinContent(1));
1626 Float_t pen(fitStatus->GetBinContent(2));
1627 Int_t dof((int)fitStatus->GetBinContent(3));
1628 Float_t beta(fitStatus->GetBinContent(4));
1629 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1630 } else if (fitStatus) { // only available in SVD fit
1631 Int_t reg((int)fitStatus->GetBinContent(1));
1632 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
1633 }
1634 }
1635 }
1636 if(tempOut) {
1637 TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirNameOut.Data())));
1638 if(pOut) {
1639 printf(" - %s out of plane converged \n", dirNameOut.Data());
1640 canvasOut->cd(j);
1641 if(i==0) canvasNominalOut->cd(j);
1642 Style(gPad, "PEARSON");
1643 pOut->DrawCopy("colz");
1644 TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirNameOut.Data())));
1645 if(rOut) {
1646 Style(rOut);
1647 printf(" > found RatioRefoldedMeasured < \n");
1648 canvasRatioMeasuredRefoldedOut->cd(j);
1649 if(i==0) canvasNominalRatioMeasuredRefoldedOut->cd(j);
1650 rOut->SetFillColor(kRed);
1651 rOut->Draw("ap");
1652 }
1653 TH1D* dvector((TH1D*)tempOut->Get("dVector"));
1654 TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
1655 TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirNameOut.Data())));
1656 TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirNameOut.Data())));
1657 if(dvector && avalue && rm && eff) {
1658 Style(dvector);
1659 Style(avalue);
1660 Style(rm);
1661 Style(eff);
1662 canvasMISC->cd(5);
1663 if(i==0) canvasNominalMISC->cd(5);
1664 Style(gPad, "SPECTRUM");
1665 dvector->DrawCopy();
1666 canvasMISC->cd(6);
1667 if(i==0) canvasNominalMISC->cd(6);
1668 Style(gPad, "SPECTRUM");
1669 avalue->DrawCopy();
1670 canvasMISC->cd(7);
1671 if(i==0) canvasNominalMISC->cd(7);
1672 Style(gPad, "PEARSON");
1673 rm->DrawCopy("colz");
1674 canvasMISC->cd(8);
1675 if(i==0) canvasNominalMISC->cd(8);
1676 eff->DrawCopy();
1677 } else if(rm && eff) {
1678 Style(rm);
1679 Style(eff);
1680 canvasMISC->cd(7);
1681 if(i==0) canvasNominalMISC->cd(7);
1682 Style(gPad, "PEARSON");
1683 rm->DrawCopy("colz");
1684 canvasMISC->cd(8);
1685 if(i==0) canvasNominalMISC->cd(8);
1686 eff->DrawCopy();
1687 }
1688 }
1689 TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirNameOut.Data())));
1690 TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirNameOut.Data())));
1691 unfoldedSpectrumOutForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
1692 TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirNameOut.Data())));
1693 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1694 if(defaultUnfoldedJetSpectrumOut) {
1695 Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
1696 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirNameOut.Data())));
1697 temp->Divide(unfoldedSpectrum);
1698 temp->SetTitle(Form("ratio nominal [%s] / [%s]", defOut.Data(), dirNameOut.Data()));
1699 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
1700 temp->GetYaxis()->SetTitle("ratio");
1701 canvasMasterOut->cd(j);
1702 if(i==0) canvasNominalMasterOut->cd(j);
1703 temp->GetYaxis()->SetRangeUser(0., 2.);
1704 temp->DrawCopy();
1705 }
1706 TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirNameOut.Data())));
1707 canvasSpectraOut->cd(j);
1708 if(i==0) canvasNominalSpectraOut->cd(j);
1709 Style(gPad);
1710 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
1711 unfoldedSpectrum->DrawCopy();
1712 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
1713 inputSpectrum->DrawCopy("same");
1714 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
1715 refoldedSpectrum->DrawCopy("same");
1716 TLegend* l(AddLegend(gPad));
1717 Style(l);
1718 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1719 Float_t chi(fitStatus->GetBinContent(1));
1720 Float_t pen(fitStatus->GetBinContent(2));
1721 Int_t dof((int)fitStatus->GetBinContent(3));
1722 Float_t beta(fitStatus->GetBinContent(4));
1723 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1724 } else if (fitStatus) { // only available in SVD fit
1725 Int_t reg((int)fitStatus->GetBinContent(1));
1726 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
1727 }
1728 }
1729 }
1730 if(canvasRatio && canvasV2) {
1731 canvasRatio->cd(j);
1732 if(i==0) canvasNominalRatio->cd(j);
1733 TGraphErrors* ratioYield(GetRatio(unfoldedSpectrumInForRatio, unfoldedSpectrumOutForRatio, TString(Form("ratio [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
1734 Double_t _tempx(0), _tempy(0);
1735 if(ratioYield) {
1736 Style(ratioYield);
1737 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
1738 ratioYield->GetPoint(b+1, _tempx, _tempy);
1739 ratioProfile->Fill(_tempx, _tempy);
1740 }
1741 ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
1742 ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1743 ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
1744 ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1745 ratioYield->SetFillColor(kRed);
1746 ratioYield->Draw("ap");
1747 }
1748 canvasV2->cd(j);
1749 if(i==0) canvasNominalV2->cd(j);
1750 TGraphErrors* ratioV2(GetV2(unfoldedSpectrumInForRatio,unfoldedSpectrumOutForRatio, .7, TString(Form("v_{2} [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
1751 if(ratioV2) {
1752 Style(ratioV2);
1753 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
1754 ratioV2->GetPoint(b+1, _tempx, _tempy);
1755 v2Profile->Fill(_tempx, _tempy);
1756
1757 }
1758 v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
1759 v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1760 ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
1761 ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1762 ratioV2->SetFillColor(kRed);
1763 ratioV2->Draw("ap");
1764 }
1765 }
1766 delete unfoldedSpectrumInForRatio;
1767 delete unfoldedSpectrumOutForRatio;
1768 }
1769 TFile output(out.Data(), "RECREATE");
1770 // save the canvasses
1771 canvasProfiles->cd(1);
1772 Style(ratioProfile);
1773 ratioProfile->DrawCopy();
1774 canvasProfiles->cd(2);
1775 Style(v2Profile);
1776 v2Profile->DrawCopy();
1777 SavePadToPDF(canvasProfiles);
1778 canvasProfiles->Write();
1779 SavePadToPDF(canvasIn);
1780 canvasIn->Write();
1781 if(canvasOut) {
1782 SavePadToPDF(canvasOut);
1783 canvasOut->Write();
1784 }
1785 SavePadToPDF(canvasRatioMeasuredRefoldedIn);
1786 canvasRatioMeasuredRefoldedIn->Write();
1787 if(canvasRatioMeasuredRefoldedOut) {
1788 SavePadToPDF(canvasRatioMeasuredRefoldedOut);
1789 canvasRatioMeasuredRefoldedOut->Write();
1790 }
1791 SavePadToPDF(canvasSpectraIn);
1792 canvasSpectraIn->Write();
1793 if(canvasSpectraOut) {
1794 SavePadToPDF(canvasSpectraOut);
1795 canvasSpectraOut->Write();
1796 SavePadToPDF(canvasRatio);
1797 canvasRatio->Write();
1798 SavePadToPDF(canvasV2);
1799 canvasV2->Write();
1800 }
1801 SavePadToPDF(canvasMasterIn);
1802 canvasMasterIn->Write();
1803 if(canvasMasterOut) {
1804 SavePadToPDF(canvasMasterOut);
1805 canvasMasterOut->Write();
1806 }
1807 SavePadToPDF(canvasMISC);
1808 canvasMISC->Write();
1809 // save the nomial canvasses
1810 SavePadToPDF(canvasNominalIn);
1811 canvasNominalIn->Write();
1812 if(canvasNominalOut) {
1813 SavePadToPDF(canvasNominalOut);
1814 canvasNominalOut->Write();
1815 }
1816 SavePadToPDF(canvasNominalRatioMeasuredRefoldedIn);
1817 canvasNominalRatioMeasuredRefoldedIn->Write();
1818 if(canvasNominalRatioMeasuredRefoldedOut) {
1819 SavePadToPDF(canvasNominalRatioMeasuredRefoldedOut);
1820 canvasNominalRatioMeasuredRefoldedOut->Write();
1821 }
1822 SavePadToPDF(canvasNominalSpectraIn);
1823 canvasNominalSpectraIn->Write();
1824 if(canvasNominalSpectraOut) {
1825 SavePadToPDF(canvasNominalSpectraOut);
1826 canvasNominalSpectraOut->Write();
1827 SavePadToPDF(canvasNominalRatio);
1828 canvasNominalRatio->Write();
1829 SavePadToPDF(canvasNominalV2);
1830 canvasNominalV2->Write();
1831 }
1832 SavePadToPDF(canvasNominalMasterIn);
1833 canvasNominalMasterIn->Write();
1834 if(canvasNominalMasterOut) {
1835 SavePadToPDF(canvasNominalMasterOut);
1836 canvasNominalMasterOut->Write();
1837 }
1838 SavePadToPDF(canvasNominalMISC);
1839 canvasNominalMISC->Write();
1840
1841 output.Write();
1842 output.Close();
1843}
1844//_____________________________________________________________________________
87233f72 1845void AliJetFlowTools::PostProcess(TString def, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
d7ec324f 1846{
1847 // go through the output file and perform post processing routines
1848 // can either be performed in one go with the unfolding, or at a later stage
53547ff2 1849 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
d7ec324f 1850 TFile readMe(in.Data(), "READ"); // open file read-only
1851 if(readMe.IsZombie()) {
1852 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1853 return;
1854 }
1855 printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
1856 readMe.ls();
1857 TList* listOfKeys((TList*)readMe.GetListOfKeys());
1858 if(!listOfKeys) {
1859 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
1860 return;
1861 }
1862 // prepare necessary canvasses
53547ff2 1863 TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
5e11c41c 1864 TCanvas* canvasOut(0x0);
486fb24e 1865 if(fDphiUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
53547ff2 1866 TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
5e11c41c 1867 TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
486fb24e 1868 if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
53547ff2 1869 TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
5e11c41c 1870 TCanvas* canvasSpectraOut(0x0);
486fb24e 1871 if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
5e11c41c 1872 TCanvas* canvasRatio(0x0);
486fb24e 1873 if(fDphiUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
5e11c41c 1874 TCanvas* canvasV2(0x0);
486fb24e 1875 if(fDphiUnfolding) canvasV2 = new TCanvas("V2", "V2");
53547ff2
RAB		 1876 TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
1877 TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
5e11c41c 1878 TCanvas* canvasMasterOut(0x0);
486fb24e 1879 if(fDphiUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
1880 (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
d7ec324f 1881 TDirectoryFile* defDir(0x0);
87233f72 1882 TCanvas* canvasProfiles(new TCanvas("canvasProfiles", "canvasProfiles"));
1883 canvasProfiles->Divide(2);
1884 TProfile* ratioProfile(new TProfile("ratioProfile", "ratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1885 TProfile* v2Profile(new TProfile("v2Profile", "v2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
d7ec324f 1886 // get an estimate of the number of outputs and find the default set
1887 Int_t cacheMe(0);
1888 for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
1889 if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
1890 if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1891 cacheMe++;
1892 }
1893 }
f3ba6c8e 1894 Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%columns)>0));
53547ff2 1895 canvasIn->Divide(columns, rows);
5e11c41c 1896 if(canvasOut) canvasOut->Divide(columns, rows);
53547ff2 1897 canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
5e11c41c 1898 if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
53547ff2 1899 canvasSpectraIn->Divide(columns, rows);
5e11c41c 1900 if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
1901 if(canvasRatio) canvasRatio->Divide(columns, rows);
1902 if(canvasV2) canvasV2->Divide(columns, rows);
d7ec324f 1903
53547ff2 1904 canvasMasterIn->Divide(columns, rows);
5e11c41c 1905 if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
d7ec324f 1906 // extract the default output
ab474fb0 1907 TH1D* defaultUnfoldedJetSpectrumIn(0x0);
1908 TH1D* defaultUnfoldedJetSpectrumOut(0x0);
d7ec324f 1909 if(defDir) {
1910 TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
1911 TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
ab474fb0 1912 if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
1913 if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
d7ec324f 1914 printf(" > succesfully extracted default results < \n");
1915 }
1916
1917 // loop through the directories, only plot the graphs if the deconvolution converged
1918 TDirectoryFile* tempDir(0x0);
1919 TDirectoryFile* tempIn(0x0);
1920 TDirectoryFile* tempOut(0x0);
1921 for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
ab474fb0 1922 // read the directory by its name
d7ec324f 1923 tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1924 if(!tempDir) continue;
1925 TString dirName(tempDir->GetName());
ab474fb0 1926 // try to read the in- and out of plane subdirs
d7ec324f 1927 tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
1928 tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
1929 j++;
486fb24e 1930 if(!tempIn) { // attempt to get MakeAU output
1931 TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
1932 TCanvas* tempPearson(new TCanvas(Form("pearson_%s", dirName.Data()), Form("pearson_%s", dirName.Data())));
1933 tempPearson->Divide(4,2);
1934 TCanvas* tempRatio(new TCanvas(Form("ratio_%s", dirName.Data()), Form("ratio_%s", dirName.Data())));
1935 tempRatio->Divide(4,2);
1936 TCanvas* tempResult(new TCanvas(Form("result_%s", dirName.Data()), Form("result_%s", dirName.Data())));
1937 tempResult->Divide(4,2);
1938 for(Int_t q(0); q < 8; q++) {
1939 tempIn = (TDirectoryFile*)tempDir->Get(Form("%s___%s", stringArray[q].Data(), dirName.Data()));
1940 if(tempIn) {
1941 // to see if the unfolding converged try to extract the pearson coefficients
1942 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
1943 if(pIn) {
1944 printf(" - %s in plane converged \n", dirName.Data());
1945 tempPearson->cd(1+q);
1946 Style(gPad, "PEARSON");
1947 pIn->DrawCopy("colz");
1948 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
1949 if(rIn) {
1950 Style(rIn);
1951 printf(" > found RatioRefoldedMeasured < \n");
1952 tempRatio->cd(q+1);
1953 rIn->SetFillColor(kRed);
1954 rIn->Draw("ap");
1955 }
1956 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
1957 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
1958 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
1959 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
1960 if(dvector && avalue && rm && eff) {
1961 Style(dvector);
1962 Style(avalue);
1963 Style(rm);
1964 Style(eff);
1965 canvasMISC->cd(1);
1966 Style(gPad, "SPECTRUM");
1967 dvector->DrawCopy();
1968 canvasMISC->cd(2);
1969 Style(gPad, "SPECTRUM");
1970 avalue->DrawCopy();
1971 canvasMISC->cd(3);
1972 Style(gPad, "PEARSON");
1973 rm->DrawCopy("colz");
1974 canvasMISC->cd(4);
1975 eff->DrawCopy();
1976 } else if(rm && eff) {
1977 Style(rm);
1978 Style(eff);
1979 canvasMISC->cd(3);
1980 Style(gPad, "PEARSON");
1981 rm->DrawCopy("colz");
1982 canvasMISC->cd(4);
1983 eff->DrawCopy();
1984 }
1985 }
1986 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
1987 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
1988 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
1989 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
1990 if(defaultUnfoldedJetSpectrumIn) {
1991 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
1992 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
1993 temp->Divide(unfoldedSpectrum);
1994 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
f3ba6c8e 1995 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
486fb24e 1996 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
1997 canvasMasterIn->cd(j);
1998 temp->GetYaxis()->SetRangeUser(0., 2);
1999 temp->DrawCopy();
2000 }
2001 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
2002 tempResult->cd(q+1);
2003 Style(gPad);
2004 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2005 unfoldedSpectrum->DrawCopy();
2006 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2007 inputSpectrum->DrawCopy("same");
2008 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2009 refoldedSpectrum->DrawCopy("same");
2010 TLegend* l(AddLegend(gPad));
2011 Style(l);
2012 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2013 Float_t chi(fitStatus->GetBinContent(1));
2014 Float_t pen(fitStatus->GetBinContent(2));
2015 Int_t dof((int)fitStatus->GetBinContent(3));
2016 Float_t beta(fitStatus->GetBinContent(4));
2017 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2018 } else if (fitStatus) { // only available in SVD fit
2019 Int_t reg((int)fitStatus->GetBinContent(1));
2020 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2021 }
2022 }
2023 }
2024 }
2025 }
d7ec324f 2026 if(tempIn) {
2027 // to see if the unfolding converged try to extract the pearson coefficients
2028 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
2029 if(pIn) {
2030 printf(" - %s in plane converged \n", dirName.Data());
2031 canvasIn->cd(j);
2032 Style(gPad, "PEARSON");
2033 pIn->DrawCopy("colz");
2034 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2035 if(rIn) {
53547ff2 2036 Style(rIn);
d7ec324f 2037 printf(" > found RatioRefoldedMeasured < \n");
2038 canvasRatioMeasuredRefoldedIn->cd(j);
ab474fb0 2039 rIn->SetFillColor(kRed);
2040 rIn->Draw("ap");
d7ec324f 2041 }
2042 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2043 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2044 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
2045 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
2046 if(dvector && avalue && rm && eff) {
53547ff2
RAB		 2047 Style(dvector);
2048 Style(avalue);
2049 Style(rm);
2050 Style(eff);
d7ec324f 2051 canvasMISC->cd(1);
2052 Style(gPad, "SPECTRUM");
2053 dvector->DrawCopy();
2054 canvasMISC->cd(2);
2055 Style(gPad, "SPECTRUM");
2056 avalue->DrawCopy();
2057 canvasMISC->cd(3);
2058 Style(gPad, "PEARSON");
2059 rm->DrawCopy("colz");
2060 canvasMISC->cd(4);
2061 eff->DrawCopy();
53547ff2
RAB		 2062 } else if(rm && eff) {
2063 Style(rm);
2064 Style(eff);
2065 canvasMISC->cd(3);
2066 Style(gPad, "PEARSON");
2067 rm->DrawCopy("colz");
2068 canvasMISC->cd(4);
2069 eff->DrawCopy();
d7ec324f 2070 }
2071 }
2072 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
2073 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
2074 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
2075 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
ab474fb0 2076 if(defaultUnfoldedJetSpectrumIn) {
2077 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2078 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
d7ec324f 2079 temp->Divide(unfoldedSpectrum);
2080 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
f3ba6c8e 2081 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
d7ec324f 2082 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
2083 canvasMasterIn->cd(j);
ab474fb0 2084 temp->GetYaxis()->SetRangeUser(0., 2);
d7ec324f 2085 temp->DrawCopy();
2086 }
2087 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
2088 canvasSpectraIn->cd(j);
2089 Style(gPad);
53547ff2 2090 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
d7ec324f 2091 unfoldedSpectrum->DrawCopy();
53547ff2 2092 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
d7ec324f 2093 inputSpectrum->DrawCopy("same");
53547ff2 2094 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
d7ec324f 2095 refoldedSpectrum->DrawCopy("same");
2096 TLegend* l(AddLegend(gPad));
53547ff2
RAB		 2097 Style(l);
2098 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2099 Float_t chi(fitStatus->GetBinContent(1));
2100 Float_t pen(fitStatus->GetBinContent(2));
d7ec324f 2101 Int_t dof((int)fitStatus->GetBinContent(3));
53547ff2
RAB		 2102 Float_t beta(fitStatus->GetBinContent(4));
2103 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2104 } else if (fitStatus) { // only available in SVD fit
2105 Int_t reg((int)fitStatus->GetBinContent(1));
2106 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
d7ec324f 2107 }
2108 }
2109 }
2110 if(tempOut) {
2111 TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
2112 if(pOut) {
2113 printf(" - %s out of plane converged \n", dirName.Data());
2114 canvasOut->cd(j);
2115 Style(gPad, "PEARSON");
2116 pOut->DrawCopy("colz");
2117 TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
2118 if(rOut) {
53547ff2 2119 Style(rOut);
d7ec324f 2120 printf(" > found RatioRefoldedMeasured < \n");
2121 canvasRatioMeasuredRefoldedOut->cd(j);
ab474fb0 2122 rOut->SetFillColor(kRed);
2123 rOut->Draw("ap");
d7ec324f 2124 }
2125 TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2126 TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2127 TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
2128 TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
2129 if(dvector && avalue && rm && eff) {
53547ff2
RAB		 2130 Style(dvector);
2131 Style(avalue);
2132 Style(rm);
2133 Style(eff);
d7ec324f 2134 canvasMISC->cd(5);
2135 Style(gPad, "SPECTRUM");
2136 dvector->DrawCopy();
2137 canvasMISC->cd(6);
2138 Style(gPad, "SPECTRUM");
2139 avalue->DrawCopy();
2140 canvasMISC->cd(7);
2141 Style(gPad, "PEARSON");
2142 rm->DrawCopy("colz");
2143 canvasMISC->cd(8);
2144 eff->DrawCopy();
53547ff2
RAB		 2145 } else if(rm && eff) {
2146 Style(rm);
2147 Style(eff);
549b5f40 2148 canvasMISC->cd(7);
53547ff2
RAB		 2149 Style(gPad, "PEARSON");
2150 rm->DrawCopy("colz");
549b5f40 2151 canvasMISC->cd(8);
53547ff2 2152 eff->DrawCopy();
d7ec324f 2153 }
2154 }
2155 TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
2156 TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
2157 TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
2158 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
ab474fb0 2159 if(defaultUnfoldedJetSpectrumOut) {
2160 Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
2161 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirName.Data())));
d7ec324f 2162 temp->Divide(unfoldedSpectrum);
2163 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
f3ba6c8e 2164 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
d7ec324f 2165 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
2166 canvasMasterOut->cd(j);
ab474fb0 2167 temp->GetYaxis()->SetRangeUser(0., 2.);
d7ec324f 2168 temp->DrawCopy();
2169 }
2170 TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
2171 canvasSpectraOut->cd(j);
2172 Style(gPad);
53547ff2 2173 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
d7ec324f 2174 unfoldedSpectrum->DrawCopy();
53547ff2 2175 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
d7ec324f 2176 inputSpectrum->DrawCopy("same");
53547ff2 2177 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
d7ec324f 2178 refoldedSpectrum->DrawCopy("same");
2179 TLegend* l(AddLegend(gPad));
53547ff2
RAB		 2180 Style(l);
2181 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2182 Float_t chi(fitStatus->GetBinContent(1));
2183 Float_t pen(fitStatus->GetBinContent(2));
2184 Int_t dof((int)fitStatus->GetBinContent(3));
2185 Float_t beta(fitStatus->GetBinContent(4));
2186 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2187 } else if (fitStatus) { // only available in SVD fit
2188 Int_t reg((int)fitStatus->GetBinContent(1));
2189 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
d7ec324f 2190 }
2191 }
2192 }
5e11c41c 2193 if(canvasRatio && canvasV2) {
2194 canvasRatio->cd(j);
2195 TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
87233f72 2196 Double_t _tempx(0), _tempy(0);
5e11c41c 2197 if(ratioYield) {
2198 Style(ratioYield);
87233f72 2199 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2200 ratioYield->GetPoint(b+1, _tempx, _tempy);
2201 ratioProfile->Fill(_tempx, _tempy);
2202 }
9f892925 2203 ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
87233f72 2204 ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
9f892925 2205 ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
87233f72 2206 ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
ab474fb0 2207 ratioYield->SetFillColor(kRed);
2208 ratioYield->Draw("ap");
5e11c41c 2209 }
2210 canvasV2->cd(j);
2211 TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
2212 if(ratioV2) {
2213 Style(ratioV2);
87233f72 2214 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2215 ratioV2->GetPoint(b+1, _tempx, _tempy);
2216 v2Profile->Fill(_tempx, _tempy);
2217
2218 }
9f892925 2219 v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
87233f72 2220 v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
5e11c41c 2221 ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
87233f72 2222 ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
ab474fb0 2223 ratioV2->SetFillColor(kRed);
2224 ratioV2->Draw("ap");
5e11c41c 2225 }
d7ec324f 2226 }
2227 }
2228 TFile output(out.Data(), "RECREATE");
87233f72 2229 canvasProfiles->cd(1);
2230 Style(ratioProfile);
2231 ratioProfile->DrawCopy();
2232 canvasProfiles->cd(2);
2233 Style(v2Profile);
2234 v2Profile->DrawCopy();
2235 SavePadToPDF(canvasProfiles);
2236 canvasProfiles->Write();
512ced40 2237 SavePadToPDF(canvasIn);
d7ec324f 2238 canvasIn->Write();
5e11c41c 2239 if(canvasOut) {
2240 SavePadToPDF(canvasOut);
2241 canvasOut->Write();
2242 }
512ced40 2243 SavePadToPDF(canvasRatioMeasuredRefoldedIn);
d7ec324f 2244 canvasRatioMeasuredRefoldedIn->Write();
5e11c41c 2245 if(canvasRatioMeasuredRefoldedOut) {
2246 SavePadToPDF(canvasRatioMeasuredRefoldedOut);
2247 canvasRatioMeasuredRefoldedOut->Write();
2248 }
512ced40 2249 SavePadToPDF(canvasSpectraIn);
d7ec324f 2250 canvasSpectraIn->Write();
5e11c41c 2251 if(canvasSpectraOut) {
2252 SavePadToPDF(canvasSpectraOut);
2253 canvasSpectraOut->Write();
2254 SavePadToPDF(canvasRatio);
2255 canvasRatio->Write();
2256 SavePadToPDF(canvasV2);
2257 canvasV2->Write();
2258 }
512ced40 2259 SavePadToPDF(canvasMasterIn);
d7ec324f 2260 canvasMasterIn->Write();
5e11c41c 2261 if(canvasMasterOut) {
2262 SavePadToPDF(canvasMasterOut);
2263 canvasMasterOut->Write();
2264 }
512ced40 2265 SavePadToPDF(canvasMISC);
d7ec324f 2266 canvasMISC->Write();
2267 output.Write();
2268 output.Close();
2269}
2270//_____________________________________________________________________________
4292ca60 2271Bool_t AliJetFlowTools::SetRawInput (
2272 TH2D* detectorResponse, // detector response matrix
2273 TH1D* jetPtIn, // in plane jet spectrum
2274 TH1D* jetPtOut, // out of plane jet spectrum
2275 TH1D* dptIn, // in plane delta pt distribution
2276 TH1D* dptOut, // out of plane delta pt distribution
2277 Int_t eventCount) {
2278 // set input histograms manually
2279 fDetectorResponse = detectorResponse;
2280 fSpectrumIn = jetPtIn;
2281 fSpectrumOut = jetPtOut;
2282 fDptInDist = dptIn;
2283 fDptOutDist = dptOut;
2284 fRawInputProvided = kTRUE;
2285 // check if all data is provided
2286 if(!fDetectorResponse) {
2287 printf(" fDetectorResponse not found \n ");
2288 return kFALSE;
2289 }
2290 // check if the pt bin for true and rec have been set
2291 if(!fBinsTrue || !fBinsRec) {
2292 printf(" No true or rec bins set, please set binning ! \n");
2293 return kFALSE;
2294 }
2295 if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
2296 // procedures, these profiles will be nonsensical, user is responsible
2297 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2298 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2299 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2300 }
2301 // normalize spectra to event count if requested
2302 if(fNormalizeSpectra) {
2303 fEventCount = eventCount;
2304 if(fEventCount > 0) {
2305 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
2306 fSpectrumOut->Sumw2();
2307 fSpectrumIn->Scale(1./((double)fEventCount));
2308 fSpectrumOut->Scale(1./((double)fEventCount));
2309 }
2310 }
2311 if(!fNormalizeSpectra && fEventCount > 0) {
2312 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
2313 fSpectrumOut->Sumw2();
2314 fSpectrumIn->Scale(1./((double)fEventCount));
2315 fSpectrumOut->Scale(1./((double)fEventCount));
2316 }
2317 fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
2318 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
f3ba6c8e 2319 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
2320 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
4292ca60 2321 fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
2322 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
f3ba6c8e 2323 fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
2324 fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
4292ca60 2325
2326 return kTRUE;
2327}
2328//_____________________________________________________________________________
2329TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
dc1455ee 2330{
ef12d5a5 2331 // return ratio of h1 / h2
2332 // histograms can have different binning. errors are propagated as uncorrelated
20abfcc4 2333 if(!(h1 && h2) ) {
2334 printf(" GetRatio called with NULL argument(s) \n ");
2335 return 0x0;
2336 }
ad04a83c 2337 Int_t j(0);
4292ca60 2338 TGraphErrors *gr = new TGraphErrors();
f3ba6c8e 2339 gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
ef12d5a5 2340 Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
7bf39895 2341 TH1* dud((TH1*)h1->Clone("dud"));
486fb24e 2342 dud->Sumw2();
2343 h1->Sumw2();
2344 h2->Sumw2();
2345 if(!dud->Divide(h1, h2)) {
2346 printf(" > ROOT failed to divide two histogams, dividing manually \n < ");
7bf39895 2347 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
2348 binCent = h1->GetXaxis()->GetBinCenter(i);
2349 if(xmax > 0. && binCent > xmax) continue;
2350 j = h2->FindBin(binCent);
2351 binWidth = h1->GetXaxis()->GetBinWidth(i);
2352 if(h2->GetBinContent(j) > 0.) {
2353 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
486fb24e 2354 /* original propagation of uncertainty changed 08012014
dc1455ee 2355 Double_t A = 1./h2->GetBinContent(j)*h1->GetBinError(i);
4292ca60 2356 Double_t B = 0.;
2357 if(h2->GetBinError(j)>0.) {
2358 B = -1.*h1->GetBinContent(i)/(h2->GetBinContent(j)*h2->GetBinContent(j))*h2->GetBinError(j);
2359 error2 = A*A + B*B;
ab474fb0 2360 } else error2 = A*A; */
7bf39895 2361 Double_t A = h1->GetBinError(i)/h1->GetBinContent(i);
2362 Double_t B = h2->GetBinError(i)/h2->GetBinContent(i);
2363 error2 = ratio*ratio*A*A+ratio*ratio*B*B;
2364 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
486fb24e 2365 gr->SetPoint(gr->GetN(), binCent, ratio);
2366 gr->SetPointError(gr->GetN()-1, 0.5*binWidth, error2);
7bf39895 2367 }
2368 }
2369 } else {
486fb24e 2370 printf( " > using ROOT to divide two histograms < \n");
7bf39895 2371 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
2372 binCent = dud->GetXaxis()->GetBinCenter(i);
2373 if(xmax > 0. && binCent > xmax) continue;
2374 binWidth = dud->GetXaxis()->GetBinWidth(i);
2375 gr->SetPoint(gr->GetN(),binCent,dud->GetBinContent(i));
2376 gr->SetPointError(gr->GetN()-1,0.5*binWidth,dud->GetBinError(i));
dc1455ee 2377 }
2378 }
7bf39895 2379
ad04a83c 2380 if(appendFit) {
4292ca60 2381 TF1* fit(new TF1("lin", "pol0", 10, 100));
ad04a83c 2382 gr->Fit(fit);
2383 }
4292ca60 2384 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
7bf39895 2385 if(dud) delete dud;
dc1455ee 2386 return gr;
2387}
2388//_____________________________________________________________________________
ef12d5a5 2389TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
2390{
2391 // get v2 from difference of in plane, out of plane yield
2392 // h1 must hold the in-plane yield, h2 holds the out of plane yield
549b5f40
RAB		 2393 // different binning is allowed but will mean that the error
2394 // propagation is unreliable
ef12d5a5 2395 // r is the event plane resolution for the chosen centrality
2396 if(!(h1 && h2) ) {
2397 printf(" GetV2 called with NULL argument(s) \n ");
2398 return 0x0;
2399 }
2400 Int_t j(0);
2401 TGraphErrors *gr = new TGraphErrors();
f3ba6c8e 2402 gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
ef12d5a5 2403 Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
2404 Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
2405 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
2406 binCent = h1->GetXaxis()->GetBinCenter(i);
2407 j = h2->FindBin(binCent);
2408 binWidth = h1->GetXaxis()->GetBinWidth(i);
2409 if(h2->GetBinContent(j) > 0.) {
2410 in = h1->GetBinContent(i);
2411 ein = h1->GetBinError(i);
2412 out = h2->GetBinContent(j);
2413 eout = h2->GetBinError(j);
2414 ratio = pre*((in-out)/(in+out));
549b5f40 2415 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 2416 if(error2 > 0) error2 = TMath::Sqrt(error2);
2417 gr->SetPoint(gr->GetN(),binCent,ratio);
2418 gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
2419 }
2420 }
2421 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
2422 return gr;
2423}
2424//_____________________________________________________________________________
549b5f40
RAB		 2425void AliJetFlowTools::WriteObject(TObject* object, TString suffix, Bool_t kill) {
2426 // write object, if a unique identifier is given the object is cloned
2427 // and the clone is saved. setting kill to true will delete the original obect from the heap
4292ca60 2428 if(!object) {
2429 printf(" > WriteObject:: called with NULL arguments \n ");
2430 return;
549b5f40
RAB		 2431 } else if(!strcmp("", suffix.Data())) object->Write();
2432 else {
2433 TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
2434 newObject->Write();
2435 }
2436 if(kill) delete object;
4292ca60 2437}
2438//_____________________________________________________________________________
2439TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
2440 // construt a delta pt response matrix from supplied dpt distribution
2441 // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
2442 // do a weighted rebinning to a (coarser) dpt distribution
2443 // be careful with the binning of the dpt response: it should be equal to that
2444 // of the response matrix, otherwise dpt and response matrices cannot be multiplied
2445 //
2446 // the response matrix will be square and have the same binning
2447 // (min, max and granularity) of the input histogram
2448 Int_t bins(dpt->GetXaxis()->GetNbins()); // number of bins, will also be no of rows, columns
2449 Double_t _bins[bins+1]; // prepare array with bin borders
2450 for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
2451 _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1); // get upper edge
2452 TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
2453 for(Int_t j(0); j < bins+1; j++) { // loop on pt true slices j
2454 Bool_t skip = kFALSE;
2455 for(Int_t k(0); k < bins+1; k++) { // loop on pt gen slices k
2456 (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
2457 if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
2458 }
2459 }
2460 return res;
2461}
ef12d5a5 2462//_____________________________________________________________________________
2463TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
2464 if(!binsTrue || !binsRec) {
2465 printf(" > GetUnityResponse:: function called with NULL arguments < \n");
2466 return 0x0;
2467 }
2468 TString name(Form("unityResponse_%s", suffix.Data()));
2469 TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
2470 for(Int_t i(0); i < binsTrue->GetSize(); i++) {
2471 for(Int_t j(0); j < binsRec->GetSize(); j++) {
2472 if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
2473 }
2474 }
2475 return unity;
2476}
4292ca60 2477//_____________________________________________________________________________
549b5f40 2478void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
4292ca60 2479 // save configuration parameters to histogram
549b5f40 2480 TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
ef12d5a5 2481 summary->SetBinContent(1, fBetaIn);
2482 summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
2483 summary->SetBinContent(2, fBetaOut);
2484 summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
2485 summary->SetBinContent(3, fCentralityBin);
2486 summary->GetXaxis()->SetBinLabel(3, "fCentralityBin");
2487 summary->SetBinContent(4, (int)convergedIn);
2488 summary->GetXaxis()->SetBinLabel(4, "convergedIn");
2489 summary->SetBinContent(5, (int)convergedOut);
2490 summary->GetXaxis()->SetBinLabel(5, "convergedOut");
2491 summary->SetBinContent(6, (int)fAvoidRoundingError);
2492 summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
2493 summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
2494 summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
2495 summary->SetBinContent(8, (int)fPrior);
2496 summary->GetXaxis()->SetBinLabel(8, "fPrior");
2497 summary->SetBinContent(9, fSVDRegIn);
2498 summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
2499 summary->SetBinContent(10, fSVDRegOut);
2500 summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
2501 summary->SetBinContent(11, (int)fSVDToy);
2502 summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
2503 summary->SetBinContent(12, fJetRadius);
2504 summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
2505 summary->SetBinContent(13, (int)fNormalizeSpectra);
2506 summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
549b5f40
RAB		 2507 summary->SetBinContent(14, (int)fSmoothenPrior);
2508 summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
ef12d5a5 2509 summary->SetBinContent(15, (int)fTestMode);
2510 summary->GetXaxis()->SetBinLabel(15, "fTestMode");
2511 summary->SetBinContent(16, (int)fUseDetectorResponse);
2512 summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
549b5f40
RAB		 2513 summary->SetBinContent(17, fBayesianIterIn);
2514 summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
2515 summary->SetBinContent(18, fBayesianIterOut);
2516 summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
2517 summary->SetBinContent(19, fBayesianSmoothIn);
2518 summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
2519 summary->SetBinContent(20, fBayesianSmoothOut);
2520 summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
4292ca60 2521}
2522//_____________________________________________________________________________
2523void AliJetFlowTools::ResetAliUnfolding() {
2524 // ugly function: reset all unfolding parameters
2525 TVirtualFitter* fitter(TVirtualFitter::GetFitter());
2526 if(fitter) {
2527 printf(" > Found fitter, will delete it < \n");
2528 delete fitter;
2529 }
d7ec324f 2530 if(gMinuit) {
2531 printf(" > Found gMinuit, will re-create it < \n");
2532 delete gMinuit;
2533 gMinuit = new TMinuit;
2534 }
4292ca60 2535 AliUnfolding::fgCorrelationMatrix = 0;
2536 AliUnfolding::fgCorrelationMatrixSquared = 0;
2537 AliUnfolding::fgCorrelationCovarianceMatrix = 0;
2538 AliUnfolding::fgCurrentESDVector = 0;
2539 AliUnfolding::fgEntropyAPriori = 0;
2540 AliUnfolding::fgEfficiency = 0;
2541 AliUnfolding::fgUnfoldedAxis = 0;
2542 AliUnfolding::fgMeasuredAxis = 0;
2543 AliUnfolding::fgFitFunction = 0;
2544 AliUnfolding::fgMaxInput = -1;
2545 AliUnfolding::fgMaxParams = -1;
2546 AliUnfolding::fgOverflowBinLimit = -1;
2547 AliUnfolding::fgRegularizationWeight = 10000;
2548 AliUnfolding::fgSkipBinsBegin = 0;
2549 AliUnfolding::fgMinuitStepSize = 0.1;
2550 AliUnfolding::fgMinuitPrecision = 1e-6;
2551 AliUnfolding::fgMinuitMaxIterations = 1000000;
2552 AliUnfolding::fgMinuitStrategy = 1.;
2553 AliUnfolding::fgMinimumInitialValue = kFALSE;
2554 AliUnfolding::fgMinimumInitialValueFix = -1;
2555 AliUnfolding::fgNormalizeInput = kFALSE;
2556 AliUnfolding::fgNotFoundEvents = 0;
2557 AliUnfolding::fgSkipBin0InChi2 = kFALSE;
2558 AliUnfolding::fgBayesianSmoothing = 1;
2559 AliUnfolding::fgBayesianIterations = 10;
2560 AliUnfolding::fgDebug = kFALSE;
2561 AliUnfolding::fgCallCount = 0;
2562 AliUnfolding::fgPowern = 5;
2563 AliUnfolding::fChi2FromFit = 0.;
2564 AliUnfolding::fPenaltyVal = 0.;
2565 AliUnfolding::fAvgResidual = 0.;
2566 AliUnfolding::fgPrintChi2Details = 0;
2567 AliUnfolding::fgCanvas = 0;
2568 AliUnfolding::fghUnfolded = 0;
2569 AliUnfolding::fghCorrelation = 0;
2570 AliUnfolding::fghEfficiency = 0;
2571 AliUnfolding::fghMeasured = 0;
2572 AliUnfolding::SetMinuitStepSize(1.);
2573 AliUnfolding::SetMinuitPrecision(1e-6);
2574 AliUnfolding::SetMinuitMaxIterations(100000);
2575 AliUnfolding::SetMinuitStrategy(2.);
2576 AliUnfolding::SetDebug(1);
2577}
2578//_____________________________________________________________________________
f3ba6c8e 2579TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) const {
4292ca60 2580 // protect heap by adding unique qualifier to name
2581 if(!protect) return 0x0;
2582 TH1D* p = (TH1D*)protect->Clone();
ef12d5a5 2583 TString tempString(fActiveString);
2584 tempString+=suffix;
2585 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
2586 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 2587 if(kill) delete protect;
2588 return p;
2589}
2590//_____________________________________________________________________________
f3ba6c8e 2591TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) const {
4292ca60 2592 // protect heap by adding unique qualifier to name
2593 if(!protect) return 0x0;
2594 TH2D* p = (TH2D*)protect->Clone();
ef12d5a5 2595 TString tempString(fActiveString);
2596 tempString+=suffix;
2597 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
2598 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 2599 if(kill) delete protect;
2600 return p;
2601}
2602//_____________________________________________________________________________
f3ba6c8e 2603TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) const {
4292ca60 2604 // protect heap by adding unique qualifier to name
2605 if(!protect) return 0x0;
2606 TGraphErrors* p = (TGraphErrors*)protect->Clone();
ef12d5a5 2607 TString tempString(fActiveString);
2608 tempString+=suffix;
2609 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
2610 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 2611 if(kill) delete protect;
2612 return p;
2613}
2614//_____________________________________________________________________________
486fb24e 2615void AliJetFlowTools::MakeAU() {
2616 // === azimuthal unfolding ===
2617 //
2618 // unfolds the spectrum in delta phi bins, extracts the yield per bin, and does a fit
2619 // in transverse momentum and azimuthal correlation space to extract v2 params
2620 // settings are equal to the ones used for 'Make()'
2621 //
2622 // basic steps that are followed:
2623 // 1) rebin the raw output of the jet task to the desired binnings
2624 // 2) calls the unfolding routine
2625 // 3) writes output to file
2626 // can be repeated multiple times with different configurations
2627
2628 Int_t low[] = {1, 6, 11, 16, 21, 26, 31, 36};
2629 Int_t up[] = {5, 10, 15, 20, 25, 30, 35, 40};
2630 TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
2631 TH1D* dPtdPhi[6];
2632 for(Int_t i(0); i < 6; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), 8, 0, TMath::Pi());
2633
2634 for(Int_t i(0); i < 8; i++) {
2635 // 1) manipulation of input histograms
2636 // check if the input variables are present
2637 if(!PrepareForUnfolding(low[i], up[i])) return;
2638 // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
2639 // parts of the spectrum can end up in over or underflow bins
2640 TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, Form("resized_%s", stringArray[i].Data()), kFALSE);
2641
2642 // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
2643 // the template will be used as a prior for the chi2 unfolding
2644 TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, stringArray[i], kFALSE);
2645
2646 // get the full response matrix from the dpt and the detector response
2647 fDetectorResponse = NormalizeTH2D(fDetectorResponse);
2648 // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
2649 // so that unfolding should return the initial spectrum
2650 if(fUseDptResponse && fUseDetectorResponse) fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
2651 else if (fUseDptResponse && !fUseDetectorResponse) fFullResponseIn = fDptIn;
2652 else if (!fUseDptResponse && fUseDetectorResponse) fFullResponseIn = fDetectorResponse;
2653 else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) return;
2654 // normalize each slide of the response to one
2655 NormalizeTH2D(fFullResponseIn);
2656 // resize to desired binning scheme
2657 TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, stringArray[i]);
2658 // get the kinematic efficiency
2659 TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
2660 kinematicEfficiencyIn->SetNameTitle(Form("kin_eff_%s", stringArray[i].Data()), Form("kin_eff_%s", stringArray[i].Data()));
2661 // suppress the errors
2662 for(Int_t j(0); j < kinematicEfficiencyIn->GetXaxis()->GetNbins(); j++) kinematicEfficiencyIn->SetBinError(1+j, 0.);
2663 TH1D* jetFindingEfficiency(0x0);
2664 if(fJetFindingEff) {
2665 jetFindingEfficiency = ProtectHeap(fJetFindingEff);
2666 jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
2667 jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
2668 }
2669 // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
2670 TH1D* unfoldedJetSpectrumIn(0x0);
2671 fActiveDir->cd(); // select active dir
2672 TDirectoryFile* dirIn = new TDirectoryFile(Form("%s___%s", stringArray[i].Data(), fActiveString.Data()), Form("%s___%s", stringArray[i].Data(), fActiveString.Data()));
2673 dirIn->cd(); // select inplane subdir
2674 // select the unfolding method
2675 unfoldedJetSpectrumIn = UnfoldWrapper(
2676 measuredJetSpectrumIn,
2677 resizedResponseIn,
2678 kinematicEfficiencyIn,
2679 measuredJetSpectrumTrueBinsIn,
2680 TString("dPtdPhiUnfolding"),
2681 jetFindingEfficiency);
2682 if(i==5) {
2683 resizedResponseIn->SetNameTitle(Form("ResponseMatrix_%s", stringArray[i].Data()), Form("response matrix %s", stringArray[i].Data()));
f3ba6c8e 2684 resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
2685 resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
486fb24e 2686 resizedResponseIn = ProtectHeap(resizedResponseIn);
2687 resizedResponseIn->Write();
2688 kinematicEfficiencyIn->SetNameTitle(Form("KinematicEfficiency_%s", stringArray[i].Data()), Form("Kinematic efficiency, %s", stringArray[i].Data()));
2689 kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
2690 kinematicEfficiencyIn->Write();
2691 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
2692 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
2693 fDetectorResponse->Write();
2694 // optional histograms
2695 if(fSaveFull) {
2696 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", Form("[INPUT] Jet spectrum, %s", stringArray[i].Data()));
2697 fSpectrumIn->Write();
2698 fDptInDist->SetNameTitle("[ORIG]DeltaPt", Form("#delta p_{T} distribution, %s", stringArray[i].Data()));
2699 fDptInDist->Write();
2700 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix", Form("#delta p_{T} matrix, %s", stringArray[i].Data()));
2701 fDptIn->Write();
2702 fFullResponseIn->SetNameTitle("ResponseMatrix", Form("Response matrix, %s", stringArray[i].Data()));
2703 fFullResponseIn->Write();
2704 }
2705 }
2706 fActiveDir->cd();
2707 fDeltaPtDeltaPhi->Write();
2708 fJetPtDeltaPhi->Write();
2709
2710 TH1D* dud(ProtectHeap(unfoldedJetSpectrumIn, kTRUE, stringArray[i]));;
2711 Double_t integralError(0);
2712 for(Int_t j(0); j < 6; j++) {
2713 // get the integrated
2714 Double_t integral(dud->IntegralAndError(2*j+1, 2*j+3, integralError));
2715 dPtdPhi[j]->SetBinContent(i+1, integral);
2716 dPtdPhi[j]->SetBinError(i+1, integralError);
2717 }
2718 dud->Write();
2719 // save the current state of the unfolding object
2720 SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, kFALSE);
2721 }
2722 TF1* fourier = new TF1("fourier", "[0]*(1.+0.5*[1]*(TMath::Cos(2.*x)))", 0, TMath::Pi());
2723 TH1D* v2(new TH1D("v2FromFit", "v2FromFit", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2724 for(Int_t i(0); i < 6; i++) {
2725 dPtdPhi[i]->Fit(fourier, "VI");
2726 v2->SetBinContent(1+i, fourier->GetParameter(1));
2727 v2->SetBinError(1+i, fourier->GetParError(1));
2728 dPtdPhi[i]->Write();
2729 }
2730 v2->Write();
2731}
2732//_____________________________________________________________________________