]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Tasks/AliJetFlowTools.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[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
ef12d5a5 28// libraries must be present on the system (see http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html).
29// A test mode is available in which the spectrum is unfolded with a generated unity response
30// matrix.
4292ca60 31//
32// The weak spot of this class is the function PrepareForUnfolding, which will read
33// output from two output files and expects histograms with certain names and binning.
34// Unfolding methods itself are general and should be able to handle any input, therefore one
35// can forgo the PrepareForUnfolding method, and supply necessary input information via the
36// SetRawInput() method
37//
38// to see an example of how to use this class, see $ALICE_ROOT/PWGCF/FLOW/macros/jetFlowTools.C
dc1455ee 39
40// root includes
41#include "TF1.h"
42#include "TH1D.h"
43#include "TH2D.h"
d7ec324f 44#include "THStack.h"
53547ff2 45#include "TGraph.h"
4292ca60 46#include "TGraphErrors.h"
d7ec324f 47#include "TCanvas.h"
48#include "TLegend.h"
dc1455ee 49#include "TArrayD.h"
50#include "TList.h"
51#include "TMinuit.h"
52#include "TVirtualFitter.h"
53#include "TLegend.h"
54#include "TCanvas.h"
55#include "TStyle.h"
dc1455ee 56#include "TLine.h"
ad04a83c 57#include "TMath.h"
4292ca60 58#include "TVirtualFitter.h"
ef12d5a5 59#include "TFitResultPtr.h"
4292ca60 60// aliroot includes
dc1455ee 61#include "AliUnfolding.h"
62#include "AliAnaChargedJetResponseMaker.h"
4292ca60 63// class includes
dc1455ee 64#include "AliJetFlowTools.h"
51e6bc5a 65// roo unfold includes (make sure you have these available on your system)
66#include "RooUnfold.h"
67#include "RooUnfoldResponse.h"
68#include "RooUnfoldSvd.h"
549b5f40 69#include "RooUnfoldBayes.h"
51e6bc5a 70#include "TSVDUnfold.h"
dc1455ee 71
72using namespace std;
dc1455ee 73//_____________________________________________________________________________
74AliJetFlowTools::AliJetFlowTools() :
4292ca60 75 fResponseMaker (new AliAnaChargedJetResponseMaker()),
ef12d5a5 76 fPower (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
4292ca60 77 fSaveFull (kFALSE),
dc1455ee 78 fActiveString (""),
ad04a83c 79 fActiveDir (0x0),
dc1455ee 80 fInputList (0x0),
4292ca60 81 fRefreshInput (kTRUE),
dc1455ee 82 fOutputFileName ("UnfoldedSpectra.root"),
ad04a83c 83 fOutputFile (0x0),
dc1455ee 84 fCentralityBin (0),
dc1455ee 85 fDetectorResponse (0x0),
53547ff2 86 fJetFindingEff (0x0),
4292ca60 87 fBetaIn (.1),
88 fBetaOut (.1),
549b5f40
RAB		 89 fBayesianIterIn (4),
90 fBayesianIterOut (4),
91 fBayesianSmoothIn (0.),
92 fBayesianSmoothOut (0.),
ef12d5a5 93 fAvoidRoundingError (kFALSE),
51e6bc5a 94 fUnfoldingAlgorithm (kChi2),
95 fPrior (kPriorMeasured),
dc1455ee 96 fBinsTrue (0x0),
97 fBinsRec (0x0),
ef12d5a5 98 fBinsTruePrior (0x0),
99 fBinsRecPrior (0x0),
4292ca60 100 fSVDRegIn (5),
101 fSVDRegOut (5),
51e6bc5a 102 fSVDToy (kTRUE),
103 fJetRadius (0.3),
20abfcc4 104 fEventCount (-1),
549b5f40
RAB		 105 fNormalizeSpectra (kFALSE),
106 fSmoothenPrior (kFALSE),
4292ca60 107 fFitMin (60.),
549b5f40 108 fFitMax (300.),
4292ca60 109 fFitStart (75.),
549b5f40 110 fSmoothenCounts (kTRUE),
ef12d5a5 111 fTestMode (kFALSE),
112 fNoDphi (kFALSE),
4292ca60 113 fRawInputProvided (kFALSE),
ef12d5a5 114 fEventPlaneRes (.63),
115 fUseDetectorResponse(kTRUE),
549b5f40 116 fUseDptResponse (kTRUE),
ef12d5a5 117 fTrainPower (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
RAB		 131 fFullResponseOut (0x0) { // class constructor
132 // create response maker weight function (tuned to PYTHIA spectrum)
4292ca60 133 fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0 ,200));
134}
dc1455ee 135//_____________________________________________________________________________
136void AliJetFlowTools::Make() {
ad04a83c 137 // core function of the class
4292ca60 138 // 1) rebin the raw output of the jet task to the desired binnings
ad04a83c 139 // 2) calls the unfolding routine
140 // 3) writes output to file
4292ca60 141 // can be repeated multiple times with different configurations
142
ad04a83c 143 // 1) manipulation of input histograms
dc1455ee 144 // check if the input variables are present
4292ca60 145 if(fRefreshInput) {
146 if(!PrepareForUnfolding()) {
147 printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
148 return;
149 }
dc1455ee 150 }
4292ca60 151 // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
152 // parts of the spectrum can end up in over or underflow bins
549b5f40
RAB		 153 TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
154 TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec, TString("resized_out_"), kFALSE);
4292ca60 155
549b5f40 156 // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
4292ca60 157 // the template will be used as a prior for the chi2 unfolding
549b5f40
RAB		 158 TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
159 TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
dc1455ee 160
4292ca60 161 // get the full response matrix from the dpt and the detector response
dc1455ee 162 fDetectorResponse = NormalizeTH2D(fDetectorResponse);
ef12d5a5 163 // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
164 // so that unfolding should return the initial spectrum
165 if(!fTestMode) {
549b5f40
RAB		 166 if(fUseDptResponse && fUseDetectorResponse) {
167 fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
168 fFullResponseOut = MatrixMultiplication(fDptOut, fDetectorResponse);
169 } else if (fUseDptResponse && !fUseDetectorResponse) {
170 fFullResponseIn = fDptIn;
171 fFullResponseOut = fDptOut;
172 } else if (!fUseDptResponse && fUseDetectorResponse) {
173 fFullResponseIn = fDetectorResponse;
174 fFullResponseOut = fDetectorResponse;
175 } else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) {
176 printf(" > No response, exiting ! < \n" );
177 return;
178 }
ef12d5a5 179 } else {
180 fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
181 fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
182 }
4292ca60 183 // normalize each slide of the response to one
dc1455ee 184 NormalizeTH2D(fFullResponseIn);
185 NormalizeTH2D(fFullResponseOut);
4292ca60 186 // resize to desired binning scheme
549b5f40
RAB		 187 TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
188 TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
4292ca60 189 // get the kinematic efficiency
549b5f40 190 TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
4292ca60 191 kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
549b5f40 192 TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
4292ca60 193 kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
194 // suppress the errors
51e6bc5a 195 for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
196 kinematicEfficiencyIn->SetBinError(1+i, 0.);
4292ca60 197 kinematicEfficiencyOut->SetBinError(1+i, 0.);
51e6bc5a 198 }
53547ff2
RAB		 199 TH1D* jetFindingEfficiency(0x0);
200 if(fJetFindingEff) {
201 jetFindingEfficiency = ProtectHeap(fJetFindingEff);
202 jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
203 jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
204 }
ad04a83c 205 // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
549b5f40
RAB		 206 TH1D* unfoldedJetSpectrumIn(0x0);
207 TH1D* unfoldedJetSpectrumOut(0x0);
ad04a83c 208 fActiveDir->cd(); // select active dir
209 TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
210 dirIn->cd(); // select inplane subdir
51e6bc5a 211 // select the unfolding method
212 switch (fUnfoldingAlgorithm) {
213 case kChi2 : {
549b5f40
RAB		 214 unfoldedJetSpectrumIn = UnfoldSpectrumChi2( // do the inplane unfolding
215 measuredJetSpectrumIn,
216 resizedResponseIn,
51e6bc5a 217 kinematicEfficiencyIn,
549b5f40 218 measuredJetSpectrumTrueBinsIn,
53547ff2
RAB		 219 TString("in"),
220 jetFindingEfficiency);
20abfcc4 221 printf(" > Spectrum (in plane) unfolded using kChi2 unfolding < \n");
51e6bc5a 222 } break;
549b5f40
RAB		 223 case kBayesian : {
224 unfoldedJetSpectrumIn = UnfoldSpectrumBayesian( // do the inplane unfolding
225 measuredJetSpectrumIn,
226 resizedResponseIn,
51e6bc5a 227 kinematicEfficiencyIn,
549b5f40 228 measuredJetSpectrumTrueBinsIn,
53547ff2
RAB		 229 TString("in"),
230 jetFindingEfficiency);
549b5f40
RAB		 231 printf(" > Spectrum (in plane) unfolded using kBayesian unfolding < \n");
232 } break;
233 case kBayesianAli : {
234 unfoldedJetSpectrumIn = UnfoldSpectrumBayesianAli( // do the inplane unfolding
235 measuredJetSpectrumIn,
236 resizedResponseIn,
237 kinematicEfficiencyIn,
238 measuredJetSpectrumTrueBinsIn,
239 TString("in"),
240 jetFindingEfficiency);
241 printf(" > Spectrum (in plane) unfolded using kBayesianAli unfolding < \n");
51e6bc5a 242 } break;
549b5f40
RAB		 243 case kSVD : {
244 unfoldedJetSpectrumIn = UnfoldSpectrumSVD( // do the inplane unfolding
245 measuredJetSpectrumIn,
246 resizedResponseIn,
d7ec324f 247 kinematicEfficiencyIn,
549b5f40 248 measuredJetSpectrumTrueBinsIn,
53547ff2
RAB		 249 TString("in"),
250 jetFindingEfficiency);
d7ec324f 251 printf(" > Spectrum (in plane) unfolded using kSVD unfolding < \n");
252 } break;
549b5f40
RAB		 253 case kNone : { // do nothing
254 resizedResponseIn->SetNameTitle("measuredSpectrumIn", "measured spectrum, in plane");
255 unfoldedJetSpectrumIn = ProtectHeap(measuredJetSpectrumIn, kTRUE, TString("in"));
ef12d5a5 256 } break;
257
51e6bc5a 258 default : {
259 printf(" > Selected unfolding method is not implemented yet ! \n");
260 return;
261 }
262 }
549b5f40
RAB		 263 resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
264 resizedResponseIn->SetXTitle("p_{T}^{true} [GeV/c]");
265 resizedResponseIn->SetYTitle("p_{T}^{rec} [GeV/c]");
266 resizedResponseIn = ProtectHeap(resizedResponseIn);
267 resizedResponseIn->Write();
4292ca60 268 kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
269 kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
ad04a83c 270 kinematicEfficiencyIn->Write();
271 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
4292ca60 272 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
ad04a83c 273 fDetectorResponse->Write();
4292ca60 274 // optional histograms
275 if(fSaveFull) {
276 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
277 fSpectrumIn->Write();
278 fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
279 fDptInDist->Write();
280 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
281 fDptIn->Write();
282 fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
283 fFullResponseIn->Write();
284 }
ad04a83c 285 fActiveDir->cd();
286 TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
287 dirOut->cd();
51e6bc5a 288 switch (fUnfoldingAlgorithm) {
289 case kChi2 : {
549b5f40
RAB		 290 unfoldedJetSpectrumOut = UnfoldSpectrumChi2(
291 measuredJetSpectrumOut,
292 resizedResponseOut,
51e6bc5a 293 kinematicEfficiencyOut,
549b5f40 294 measuredJetSpectrumTrueBinsOut,
53547ff2
RAB		 295 TString("out"),
296 jetFindingEfficiency);
20abfcc4 297 printf(" > Spectrum (out of plane) unfolded using kChi2 < \n");
51e6bc5a 298 } break;
549b5f40
RAB		 299 case kBayesian : {
300 unfoldedJetSpectrumOut = UnfoldSpectrumBayesian(
301 measuredJetSpectrumOut,
302 resizedResponseOut,
51e6bc5a 303 kinematicEfficiencyOut,
549b5f40 304 measuredJetSpectrumTrueBinsOut,
53547ff2
RAB		 305 TString("out"),
306 jetFindingEfficiency);
549b5f40 307 printf(" > Spectrum (out of plane) unfolded using kBayesian < \n");
51e6bc5a 308 } break;
549b5f40
RAB		 309 case kBayesianAli : {
310 unfoldedJetSpectrumOut = UnfoldSpectrumBayesianAli(
311 measuredJetSpectrumOut,
312 resizedResponseOut,
d7ec324f 313 kinematicEfficiencyOut,
549b5f40
RAB		 314 measuredJetSpectrumTrueBinsOut,
315 TString("out"),
316 jetFindingEfficiency);
317 printf(" > Spectrum (out of plane) unfolded using kBayesianAli < \n");
318 } break;
319 case kSVD : {
320 unfoldedJetSpectrumOut = UnfoldSpectrumSVD(
321 measuredJetSpectrumOut,
322 resizedResponseOut,
323 kinematicEfficiencyOut,
324 measuredJetSpectrumTrueBinsOut,
53547ff2
RAB		 325 TString("out"),
326 jetFindingEfficiency);
d7ec324f 327 printf(" > Spectrum (out of plane) unfolded using kSVD < \n");
328 } break;
549b5f40
RAB		 329 case kNone : { // do nothing
330 resizedResponseOut->SetNameTitle("measuredSpectrumOut", "measured spectrum, out plane");
331 unfoldedJetSpectrumOut = ProtectHeap(measuredJetSpectrumOut, kTRUE, TString("out"));
ef12d5a5 332 } break;
51e6bc5a 333 default : {
334 printf(" > Selected unfolding method is not implemented yet ! \n");
335 return;
336 }
337 }
549b5f40
RAB		 338 resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
339 resizedResponseOut->SetXTitle("p_{T}^{true} [GeV/c]");
340 resizedResponseOut->SetYTitle("p_{T}^{rec} [GeV/c]");
341 resizedResponseOut = ProtectHeap(resizedResponseOut);
342 resizedResponseOut->Write();
4292ca60 343 kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
344 kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
ad04a83c 345 kinematicEfficiencyOut->Write();
346 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
4292ca60 347 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
ad04a83c 348 fDetectorResponse->Write();
53547ff2 349 if(jetFindingEfficiency) jetFindingEfficiency->Write();
4292ca60 350 // optional histograms
351 if(fSaveFull) {
ef12d5a5 352 fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
4292ca60 353 fSpectrumOut->Write();
354 fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
355 fDptOutDist->Write();
356 fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
357 fDptOut->Write();
358 fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
359 fFullResponseOut->Write();
360 }
53547ff2 361
51e6bc5a 362 // write general output histograms to file
ad04a83c 363 fActiveDir->cd();
549b5f40
RAB		 364 if(unfoldedJetSpectrumIn && unfoldedJetSpectrumOut && unfoldedJetSpectrumIn && unfoldedJetSpectrumOut) {
365 TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out")));
4292ca60 366 if(ratio) {
367 ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
368 ratio->GetXaxis()->SetTitle("p_{T} [GeV/c]");
369 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
ef12d5a5 370 ratio = ProtectHeap(ratio);
4292ca60 371 ratio->Write();
372 // write histo values to RMS files if both routines converged
373 // input values are weighted by their uncertainty
374 for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
549b5f40
RAB		 375 if(unfoldedJetSpectrumIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumIn->GetBinError(i+1), 2.));
376 if(unfoldedJetSpectrumOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), unfoldedJetSpectrumOut->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumOut->GetBinError(i+1), 2.));
377 if(unfoldedJetSpectrumOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1) / unfoldedJetSpectrumOut->GetBinContent(i+1));
ef12d5a5 378 }
379 }
549b5f40 380 TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2")));
ef12d5a5 381 if(v2) {
382 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
383 v2->GetXaxis()->SetTitle("p_{T} [GeV/c]");
384 v2->GetYaxis()->SetTitle("v_{2}");
385 v2 = ProtectHeap(v2);
386 v2->Write();
387 }
549b5f40
RAB		 388 } else if (unfoldedJetSpectrumOut && unfoldedJetSpectrumIn) {
389 TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
ef12d5a5 390 if(ratio) {
391 ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
392 ratio->GetXaxis()->SetTitle("p_{T} [GeV/c]");
393 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
394 ratio = ProtectHeap(ratio);
395 ratio->Write();
396 }
549b5f40 397 TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2")));
ef12d5a5 398 if(v2) {
399 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
400 v2->GetXaxis()->SetTitle("p_{T} [GeV/c]");
401 v2->GetYaxis()->SetTitle("v_{2}");
402 v2 = ProtectHeap(v2);
403 v2->Write();
4292ca60 404 }
20abfcc4 405 }
ad04a83c 406 fDeltaPtDeltaPhi->Write();
407 fJetPtDeltaPhi->Write();
549b5f40
RAB		 408 // save the current state of the unfolding object
409 SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
dc1455ee 410}
411//_____________________________________________________________________________
549b5f40
RAB		 412TH1D* AliJetFlowTools::UnfoldSpectrumChi2(
413 const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
414 const TH2D* resizedResponse, // response matrix
415 const TH1D* kinematicEfficiency, // kinematic efficiency
416 const TH1D* measuredJetSpectrumTrueBins, // unfolding template: same binning is pt gen of response
417 const TString suffix, // suffix (in or out of plane)
418 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
dc1455ee 419{
51e6bc5a 420 // unfold the spectrum using chi2 minimization
421
4292ca60 422 // step 0) setup the static members of AliUnfolding
423 ResetAliUnfolding(); // reset from previous iteration
424 // also deletes and re-creates the global TVirtualFitter
425 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
426 if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
427 else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
428 if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
429 else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
430 AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
dc1455ee 431
549b5f40
RAB		 432 // step 1) clone all input histograms. the histograms are cloned to make sure that the original histograms
433 // stay intact. a local copy of a histogram (which only exists in the scope of this function) is
434 // denoted by the suffix 'Local'
4292ca60 435
549b5f40
RAB		 436 // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
437 TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
4292ca60 438 // unfolded local will be filled with the result of the unfolding
439 TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
440
dc1455ee 441 // full response matrix and kinematic efficiency
549b5f40 442 TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
51e6bc5a 443 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
d7ec324f 444
4292ca60 445 // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
549b5f40
RAB		 446 TH1D *priorLocal = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("priorLocal_%s", suffix.Data()));
447 // optionally, the prior can be smoothened by extrapolating the spectrum using a power law fit
448 if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
ef12d5a5 449
4292ca60 450 // step 2) start the unfolding
51e6bc5a 451 Int_t status(-1), i(0);
452 while(status < 0 && i < 100) {
4292ca60 453 // i > 0 means that the first iteration didn't converge. in that case, the result of the first
454 // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
455 if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
456 status = AliUnfolding::Unfold(
457 resizedResponseLocal, // response matrix
458 kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
549b5f40 459 measuredJetSpectrumLocal, // measured spectrum
4292ca60 460 priorLocal, // initial conditions (set NULL to use measured spectrum)
461 unfoldedLocal); // results
462 // status holds the minuit fit status (where 0 means convergence)
51e6bc5a 463 i++;
464 }
4292ca60 465 // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
466 if(status == 0 && gMinuit->fISW[1] == 3) {
467 // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
ad04a83c 468 TVirtualFitter *fitter(TVirtualFitter::GetFitter());
4292ca60 469 if(gMinuit) gMinuit->Command("SET COV");
470 TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
471 TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
472 pearson->Print();
ad04a83c 473 TH2D *hPearson(new TH2D(*pearson));
4292ca60 474 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
475 hPearson = ProtectHeap(hPearson);
ad04a83c 476 hPearson->Write();
4292ca60 477 } else status = -1;
d7ec324f 478
4292ca60 479 // step 3) refold the unfolded spectrum and save the ratio measured / refolded
480 TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
481 foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
482 unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
549b5f40 483 TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
4292ca60 484 if(ratio) {
53547ff2
RAB		 485 ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
486 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
4292ca60 487 ratio = ProtectHeap(ratio);
488 ratio->Write();
489 }
d7ec324f 490
491 // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
492 // objects are cloned using 'ProtectHeap()'
549b5f40
RAB		 493 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
494 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
495 measuredJetSpectrumLocal->Write();
d7ec324f 496
4292ca60 497 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
498 resizedResponseLocal->Write();
d7ec324f 499
4292ca60 500 unfoldedLocal = ProtectHeap(unfoldedLocal);
53547ff2 501 if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
51e6bc5a 502 unfoldedLocal->Write();
d7ec324f 503
4292ca60 504 foldedLocal = ProtectHeap(foldedLocal);
505 foldedLocal->Write();
d7ec324f 506
ef12d5a5 507 priorLocal = ProtectHeap(priorLocal);
508 priorLocal->Write();
d7ec324f 509
510 // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
53547ff2 511 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 512 fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
513 fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
514 fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
515 fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
516 fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
517 fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
53547ff2
RAB		 518 fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
519 fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
d7ec324f 520 fitStatus->Write();
549b5f40 521 return unfoldedLocal;
dc1455ee 522}
523//_____________________________________________________________________________
549b5f40
RAB		 524TH1D* AliJetFlowTools::UnfoldSpectrumSVD(
525 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
526 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
527 const TH1D* kinematicEfficiency, // kinematic efficiency
528 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
529 const TString suffix, // suffix (in, out)
530 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
51e6bc5a 531{
549b5f40
RAB		 532
533 TH1D* priorLocal( GetPrior(
534 measuredJetSpectrum, // jet pt in pt rec bins
535 resizedResponse, // full response matrix, normalized in slides of pt true
536 kinematicEfficiency, // kinematic efficiency
537 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
538 suffix, // suffix (in, out)
539 jetFindingEfficiency)); // jet finding efficiency (optional)
540 if(!priorLocal) {
541 printf(" > couldn't find prior ! < \n");
542 return 0x0;
543 } else printf(" 1) retrieved prior \n");
544
545 // go back to the 'root' directory of this instance of the SVD unfolding routine
20abfcc4 546 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
547
4292ca60 548 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
549b5f40
RAB		 549 // measured jets in pt rec binning
550 TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
551 // local copie of the response matrix
552 TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
553 // local copy of response matrix, all true slides normalized to 1
554 // this response matrix will eventually be used in the re-folding routine
555 TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
556 resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
ef12d5a5 557 // kinematic efficiency
51e6bc5a 558 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
ef12d5a5 559 // place holder histos
d7ec324f 560 TH1D *unfoldedLocalSVD(0x0);
561 TH1D *foldedLocalSVD(0x0);
4292ca60 562 cout << " 2) setup necessary input " << endl;
51e6bc5a 563 // 3) configure routine
564 RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
d7ec324f 565 cout << " step 3) configured routine " << endl;
566
4292ca60 567 // 4) get transpose matrices
549b5f40
RAB		 568 // a) get the transpose of the full response matrix
569 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
4292ca60 570 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
549b5f40
RAB		 571 // normalize it with the prior. this will ensure that high statistics bins will constrain the
572 // end result most strenuously than bins with limited number of counts
573 responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
574 cout << " 4) retrieved first transpose matrix " << endl;
4292ca60 575
576 // 5) get response for SVD unfolding
577 RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
4292ca60 578 cout << " 5) retrieved roo unfold response object " << endl;
549b5f40 579
4292ca60 580 // 6) actualy unfolding loop
549b5f40 581 RooUnfoldSvd unfoldSVD(&responseSVD, measuredJetSpectrumLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
20abfcc4 582 unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
549b5f40
RAB		 583 // correct the spectrum for the kinematic efficiency
584 unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
585
586 // get the pearson coefficients from the covariance matrix
20abfcc4 587 TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
588 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
4292ca60 589 if(pearson) {
549b5f40 590 TH2D* hPearson(new TH2D(*pearson));
4292ca60 591 pearson->Print();
d7ec324f 592 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
593 hPearson = ProtectHeap(hPearson);
4292ca60 594 hPearson->Write();
549b5f40 595 } else return 0x0; // return if unfolding didn't converge
4292ca60 596
597 // plot singular values and d_i vector
20abfcc4 598 TSVDUnfold* svdUnfold(unfoldSVD.Impl());
51e6bc5a 599 TH1* hSVal(svdUnfold->GetSV());
600 TH1D* hdi(svdUnfold->GetD());
4292ca60 601 hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
51e6bc5a 602 hSVal->SetXTitle("singular values");
603 hSVal->Write();
4292ca60 604 hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
51e6bc5a 605 hdi->SetXTitle("|d_{i}^{kreg}|");
606 hdi->Write();
4292ca60 607 cout << " plotted singular values and d_i vector " << endl;
608
609 // 7) refold the unfolded spectrum
549b5f40
RAB		 610 foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, resizedResponseLocalNorm, kinematicEfficiencyLocal);
611 TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalSVD, "ratio measured / re-folded", kTRUE));
53547ff2 612 ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
4292ca60 613 ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
614 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
615 ratio->Write();
616 cout << " 7) refolded the unfolded spectrum " << endl;
617
549b5f40
RAB		 618 // write the measured, unfolded and re-folded spectra to the output directory
619 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
620 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
621 measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
622 measuredJetSpectrumLocal->Write(); // input spectrum
d7ec324f 623 unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
624 unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
53547ff2 625 if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
4292ca60 626 unfoldedLocalSVD->Write(); // unfolded spectrum
d7ec324f 627 foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
628 foldedLocalSVD = ProtectHeap(foldedLocalSVD);
4292ca60 629 foldedLocalSVD->Write(); // re-folded spectrum
d7ec324f 630
549b5f40
RAB		 631 // save more general bookkeeeping histograms to the output directory
632 responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
4292ca60 633 responseMatrixLocalTransposePrior->SetXTitle("p_{T}^{true} [GeV/c]");
634 responseMatrixLocalTransposePrior->SetYTitle("p_{T}^{rec} [GeV/c]");
635 responseMatrixLocalTransposePrior->Write();
549b5f40
RAB		 636 priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
637 priorLocal->SetXTitle("p_{t} [GeV/c]");
638 priorLocal = ProtectHeap(priorLocal);
639 priorLocal->Write();
640 resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
641 resizedResponseLocalNorm->Write();
53547ff2
RAB		 642
643 // save some info
644 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));
645 fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
646 fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
647 fitStatus->Write();
648
549b5f40 649 return unfoldedLocalSVD;
51e6bc5a 650}
651//_____________________________________________________________________________
549b5f40
RAB		 652TH1D* AliJetFlowTools::UnfoldSpectrumBayesianAli(
653 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
654 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
655 const TH1D* kinematicEfficiency, // kinematic efficiency
656 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
657 const TString suffix, // suffix (in, out)
658 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
d7ec324f 659{
549b5f40
RAB		 660 // unfold the spectrum using the bayesian unfolding impelmented in AliUnfolding
661 // FIXME careful, not tested yet ! (06122013) FIXME
662
663 // step 0) setup the static members of AliUnfolding
664 ResetAliUnfolding(); // reset from previous iteration
665 // also deletes and re-creates the global TVirtualFitter
666 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
667 if(!strcmp("in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
668 else if(!strcmp("out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
669 else if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
670 else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
4e4f12b6
RAB		 671 AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
672
549b5f40
RAB		 673 // 1) get a prior for unfolding and clone all the input histograms
674 TH1D* priorLocal( GetPrior(
675 measuredJetSpectrum, // jet pt in pt rec bins
676 resizedResponse, // full response matrix, normalized in slides of pt true
677 kinematicEfficiency, // kinematic efficiency
678 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
679 suffix, // suffix (in, out)
680 jetFindingEfficiency)); // jet finding efficiency (optional)
681 if(!priorLocal) {
682 printf(" > couldn't find prior ! < \n");
683 return 0x0;
684 } else printf(" 1) retrieved prior \n");
685 // switch back to root dir of this unfolding procedure
686 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
687
688 // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
689 TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
690 // unfolded local will be filled with the result of the unfolding
691 TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
692
693 // full response matrix and kinematic efficiency
694 TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
695 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
696
697 // step 2) start the unfolding
698 Int_t status(-1), i(0);
699 while(status < 0 && i < 100) {
700 // i > 0 means that the first iteration didn't converge. in that case, the result of the first
701 // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
702 if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
703 status = AliUnfolding::Unfold(
704 resizedResponseLocal, // response matrix
705 kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
706 measuredJetSpectrumLocal, // measured spectrum
707 priorLocal, // initial conditions (set NULL to use measured spectrum)
708 unfoldedLocal); // results
709 // status holds the minuit fit status (where 0 means convergence)
710 i++;
711 }
712 // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
713 if(status == 0 && gMinuit->fISW[1] == 3) {
714 // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
715 TVirtualFitter *fitter(TVirtualFitter::GetFitter());
716 if(gMinuit) gMinuit->Command("SET COV");
717 TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
718 TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
719 pearson->Print();
720 TH2D *hPearson(new TH2D(*pearson));
721 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
722 hPearson = ProtectHeap(hPearson);
723 hPearson->Write();
724 } else status = -1;
725
726 // step 3) refold the unfolded spectrum and save the ratio measured / refolded
727 TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
728 foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
729 unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
730 TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
731 if(ratio) {
732 ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
733 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
734 ratio = ProtectHeap(ratio);
735 ratio->Write();
d7ec324f 736 }
549b5f40
RAB		 737
738 // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
739 // objects are cloned using 'ProtectHeap()'
740 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
741 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
742 measuredJetSpectrumLocal->Write();
743
744 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
745 resizedResponseLocal->Write();
746
747 unfoldedLocal = ProtectHeap(unfoldedLocal);
748 if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
749 unfoldedLocal->Write();
750
751 foldedLocal = ProtectHeap(foldedLocal);
752 foldedLocal->Write();
753
754 priorLocal = ProtectHeap(priorLocal);
755 priorLocal->Write();
756
757 // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
758 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));
759 fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
760 fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
761 fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
762 fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
763 fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
764 fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
765 fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
766 fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
767 fitStatus->Write();
768 return unfoldedLocal;
769}
770//_____________________________________________________________________________
771TH1D* AliJetFlowTools::UnfoldSpectrumBayesian(
772 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
773 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
774 const TH1D* kinematicEfficiency, // kinematic efficiency
775 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
776 const TString suffix, // suffix (in, out)
777 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
778{
779 // use bayesian unfolding from the RooUnfold package to unfold jet spectra
780
781 // 1) get a prior for unfolding.
782 TH1D* priorLocal( GetPrior(
783 measuredJetSpectrum, // jet pt in pt rec bins
784 resizedResponse, // full response matrix, normalized in slides of pt true
785 kinematicEfficiency, // kinematic efficiency
786 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
787 suffix, // suffix (in, out)
788 jetFindingEfficiency)); // jet finding efficiency (optional)
789 if(!priorLocal) {
790 printf(" > couldn't find prior ! < \n");
791 return 0x0;
792 } else printf(" 1) retrieved prior \n");
d7ec324f 793 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
d7ec324f 794
795 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
549b5f40
RAB		 796 // measured jets in pt rec binning
797 TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
798 // local copie of the response matrix
799 TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
800 // local copy of response matrix, all true slides normalized to 1
801 // this response matrix will eventually be used in the re-folding routine
802 TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
803 resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
d7ec324f 804 // kinematic efficiency
805 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
806 // place holder histos
549b5f40
RAB		 807 TH1D *unfoldedLocalBayes(0x0);
808 TH1D *foldedLocalBayes(0x0);
d7ec324f 809 cout << " 2) setup necessary input " << endl;
549b5f40
RAB		 810 // 4) get transpose matrices
811 // a) get the transpose of the full response matrix
812 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
d7ec324f 813 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
549b5f40
RAB		 814 // normalize it with the prior. this will ensure that high statistics bins will constrain the
815 // end result most strenuously than bins with limited number of counts
816 responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
817 // 3) get response for Bayesian unfolding
818 RooUnfoldResponse responseBayes(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedBayes_%s", suffix.Data()), Form("respCombinedBayes_%s", suffix.Data()));
d7ec324f 819
549b5f40
RAB		 820 // 4) actualy unfolding loop
821 RooUnfoldBayes unfoldBayes(&responseBayes, measuredJetSpectrumLocal, (!strcmp("in", suffix.Data())) ? fBayesianIterIn : fBayesianIterOut);
822 RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
823 unfoldedLocalBayes = (TH1D*)unfoldBayes.Hreco(errorTreatment);
824 // correct the spectrum for the kinematic efficiency
825 unfoldedLocalBayes->Divide(kinematicEfficiencyLocal);
826 // get the pearson coefficients from the covariance matrix
827 TMatrixD covarianceMatrix = unfoldBayes.Ereco(errorTreatment);
d7ec324f 828 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
d7ec324f 829 if(pearson) {
549b5f40 830 TH2D* hPearson(new TH2D(*pearson));
d7ec324f 831 pearson->Print();
832 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
833 hPearson = ProtectHeap(hPearson);
834 hPearson->Write();
549b5f40 835 } else return 0x0; // return if unfolding didn't converge
d7ec324f 836
549b5f40
RAB		 837 // 5) refold the unfolded spectrum
838 foldedLocalBayes = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalBayes, resizedResponseLocalNorm, kinematicEfficiencyLocal);
839 TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalBayes, "ratio measured / re-folded", kTRUE));
840 ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
d7ec324f 841 ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
842 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
843 ratio->Write();
844 cout << " 7) refolded the unfolded spectrum " << endl;
845
549b5f40
RAB		 846 // write the measured, unfolded and re-folded spectra to the output directory
847 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
848 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
849 measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
850 measuredJetSpectrumLocal->Write(); // input spectrum
851 unfoldedLocalBayes->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
852 unfoldedLocalBayes = ProtectHeap(unfoldedLocalBayes);
853 if(jetFindingEfficiency) unfoldedLocalBayes->Divide(jetFindingEfficiency);
854 unfoldedLocalBayes->Write(); // unfolded spectrum
855 foldedLocalBayes->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
856 foldedLocalBayes = ProtectHeap(foldedLocalBayes);
857 foldedLocalBayes->Write(); // re-folded spectrum
858
859 // save more general bookkeeeping histograms to the output directory
860 responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
d7ec324f 861 responseMatrixLocalTransposePrior->SetXTitle("p_{T}^{true} [GeV/c]");
862 responseMatrixLocalTransposePrior->SetYTitle("p_{T}^{rec} [GeV/c]");
863 responseMatrixLocalTransposePrior->Write();
549b5f40
RAB		 864 priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
865 priorLocal->SetXTitle("p_{t} [GeV/c]");
866 priorLocal = ProtectHeap(priorLocal);
867 priorLocal->Write();
868 resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
869 resizedResponseLocalNorm->Write();
53547ff2
RAB		 870
871 // save some info
872 TH1F* fitStatus(new TH1F(Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), Form("fitStatus_%s_%s", fActiveString.Data(), suffix.Data()), 1, -0.5, 0.5));
549b5f40
RAB		 873 fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fBayesianIterIn : fBayesianIterOut);
874 fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fBayesianIterIn" : "fBayesianIterOut");
53547ff2
RAB		 875 fitStatus->Write();
876
549b5f40 877 return unfoldedLocalBayes;
d7ec324f 878}
879//_____________________________________________________________________________
dc1455ee 880Bool_t AliJetFlowTools::PrepareForUnfolding()
881{
882 // prepare for unfolding
4292ca60 883 if(fRawInputProvided) return kTRUE;
dc1455ee 884 if(!fInputList) {
885 printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
886 return kFALSE;
887 }
549b5f40 888 if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
4292ca60 889 // check if the pt bin for true and rec have been set
890 if(!fBinsTrue || !fBinsRec) {
891 printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
892 return kFALSE;
dc1455ee 893 }
4292ca60 894 if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
895 // procedures, these profiles will be nonsensical, user is responsible
896 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
897 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
898 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
dc1455ee 899 }
ef12d5a5 900 if(!fTrainPower) {
549b5f40 901 // clear minuit state to avoid constraining the fit with the results of the previous iteration
ef12d5a5 902 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
903 }
dc1455ee 904 // extract the spectra
905 TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin));
ad04a83c 906 fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
907 if(!fJetPtDeltaPhi) {
dc1455ee 908 printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
909 return kFALSE;
910 }
4292ca60 911 fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
dc1455ee 912 // in plane spectrum
ef12d5a5 913 if(fNoDphi) {
914 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40);
915 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40);
916 } else {
917 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10);
918 fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40));
919 fSpectrumIn = ProtectHeap(fSpectrumIn);
920 // out of plane spectrum
921 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30);
922 fSpectrumOut = ProtectHeap(fSpectrumOut);
923 }
20abfcc4 924 // normalize spectra to event count if requested
925 if(fNormalizeSpectra) {
926 TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityBin)));
ef12d5a5 927 if(!rho) return 0x0;
928 Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
929 if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
20abfcc4 930 if(fEventCount > 0) {
4292ca60 931 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
932 fSpectrumOut->Sumw2();
ef12d5a5 933 Double_t pt(0);
d7ec324f 934 Double_t error(0); // lots of issues with the errors here ...
ef12d5a5 935 for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
4292ca60 936 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount; // normalized count
d7ec324f 937 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
4292ca60 938 fSpectrumIn->SetBinContent(1+i, pt);
ef12d5a5 939 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
d7ec324f 940 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
ef12d5a5 941 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
4292ca60 942 }
ef12d5a5 943 for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
944 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount; // normalized count
d7ec324f 945 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
ef12d5a5 946 fSpectrumOut->SetBinContent(1+i, pt);
947 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
d7ec324f 948 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
ef12d5a5 949 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
4292ca60 950 }
20abfcc4 951 }
ef12d5a5 952 if(normalizeToFullSpectrum) fEventCount = -1;
20abfcc4 953 }
dc1455ee 954 // extract the delta pt matrices
512ced40 955 TString deltaptName(Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityBin));
ad04a83c 956 fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
957 if(!fDeltaPtDeltaPhi) {
dc1455ee 958 printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
549b5f40
RAB		 959 printf(" > may be ok, depending no what you want to do < \n");
960 fRefreshInput = kTRUE;
961 return kTRUE;
dc1455ee 962 }
4292ca60 963 fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
964 // in plane delta pt distribution
ef12d5a5 965 if(fNoDphi) {
966 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40);
967 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40);
968 } else {
969 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10);
970 fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40));
971 // out of plane delta pt distribution
972 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30);
973 fDptInDist = ProtectHeap(fDptInDist);
974 fDptOutDist = ProtectHeap(fDptOutDist);
975 // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
976 }
4292ca60 977
dc1455ee 978 // create a rec - true smeared response matrix
979 TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
980 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
ad04a83c 981 Bool_t skip = kFALSE;
dc1455ee 982 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
ad04a83c 983 (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
984 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
dc1455ee 985 }
986 }
987 TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
988 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
ad04a83c 989 Bool_t skip = kFALSE;
dc1455ee 990 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
ad04a83c 991 (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
992 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
dc1455ee 993 }
994 }
995 fDptIn = new TH2D(*rfIn);
996 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
997 fDptIn->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
998 fDptIn->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
4292ca60 999 fDptIn = ProtectHeap(fDptIn);
dc1455ee 1000 fDptOut = new TH2D(*rfOut);
1001 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
1002 fDptOut->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1003 fDptOut->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
4292ca60 1004 fDptOut = ProtectHeap(fDptOut);
1005
1006 fRefreshInput = kTRUE; // force cloning of the input
dc1455ee 1007 return kTRUE;
1008}
1009//_____________________________________________________________________________
549b5f40
RAB		 1010TH1D* AliJetFlowTools::GetPrior(
1011 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
1012 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
1013 const TH1D* kinematicEfficiency, // kinematic efficiency
1014 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
1015 const TString suffix, // suffix (in, out)
1016 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
1017{
1018 // 1) get a prior for unfolding.
1019 // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
1020 TH1D* unfolded(0x0);
1021 TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
1022 dirOut->cd();
1023 switch (fPrior) { // select the prior for unfolding
1024 case kPriorChi2 : {
1025 if(fBinsTruePrior && fBinsRecPrior) { // if set, use different binning for the prior
1026 TArrayD* tempArrayTrue(fBinsTrue); // temporarily cache the original binning
1027 fBinsTrue = fBinsTruePrior; // switch binning schemes (will be used in UnfoldSpectrumChi2())
1028 TArrayD* tempArrayRec(fBinsRec);
1029 fBinsRec = fBinsRecPrior;
1030 TH1D* measuredJetSpectrumChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
1031 TH1D* measuredJetSpectrumTrueBinsChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
1032 TH2D* resizedResponseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
1033 TH1D* kinematicEfficiencyChi2(resizedResponseChi2->ProjectionX());
1034 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
1035 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
1036 unfolded= UnfoldSpectrumChi2(
1037 measuredJetSpectrumChi2,
1038 resizedResponseChi2,
1039 kinematicEfficiencyChi2,
1040 measuredJetSpectrumTrueBinsChi2, // prior for chi2 unfolding (measured)
1041 TString(Form("prior_%s", suffix.Data())));
1042 if(!unfolded) {
1043 printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1044 printf(" probably Chi2 unfolding did not converge < \n");
1045 return 0x0;
1046 }
1047 fBinsTrue = tempArrayTrue; // reset bins borders
1048 fBinsRec = tempArrayRec;
1049 // if the chi2 prior has a different binning, rebin to the true binning for the unfolding
1050 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data()))); // rebin unfolded
1051 } else {
1052 unfolded = UnfoldSpectrumChi2(
1053 measuredJetSpectrum,
1054 resizedResponse,
1055 kinematicEfficiency,
1056 measuredJetSpectrumTrueBins, // prior for chi2 unfolding (measured)
1057 TString(Form("prior_%s", suffix.Data())));
1058 if(!unfolded) {
1059 printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1060 printf(" probably Chi2 unfolding did not converge < \n");
1061 return 0x0;
1062 }
1063 }
1064 break;
1065 }
1066 case kPriorMeasured : {
1067 unfolded = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("kPriorMeasured_%s", suffix.Data())); // copy template to unfolded to use as prior
1068 }
1069 default : break;
1070 }
1071 // it can be important that the prior is smooth, this can be achieved by
1072 // extrapolating the spectrum with a fitted power law when bins are sparsely filed
1073 if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
1074 TH1D *priorLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
1075 if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
1076 return priorLocal;
1077}
1078//_____________________________________________________________________________
dc1455ee 1079TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
1080 // resize the x-axis of a th1d
1081 if(!histo) {
1082 printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1083 return NULL;
1084 }
1085 // see how many bins we need to copy
1086 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);
1087 // low is the bin number of the first new bin
1088 Int_t l = histo->GetXaxis()->FindBin(low);
1089 // set the values
1090 Double_t _x(0), _xx(0);
1091 for(Int_t i(0); i < up-low; i++) {
1092 _x = histo->GetBinContent(l+i);
1093 _xx=histo->GetBinError(l+i);
1094 resized->SetBinContent(i+1, _x);
1095 resized->SetBinError(i+1, _xx);
1096 }
1097 return resized;
1098}
1099//_____________________________________________________________________________
1100TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
1101 // resize the y-axis of a th2d
1102 if(!histo) {
1103 printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1104 return NULL;
1105 }
1106 // see how many bins we need to copy
1107 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());
1108 // assume only the y-axis has changed
1109 // low is the bin number of the first new bin
1110 Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1111 // set the values
1112 Double_t _x(0), _xx(0);
1113 for(Int_t i(0); i < x->GetSize(); i++) {
1114 for(Int_t j(0); j < y->GetSize(); j++) {
1115 _x = histo->GetBinContent(i, low+j);
1116 _xx=histo->GetBinError(i, low+1+j);
1117 resized->SetBinContent(i, j, _x);
1118 resized->SetBinError(i, j, _xx);
1119 }
1120 }
1121 return resized;
1122}
1123//_____________________________________________________________________________
512ced40 1124TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo, Bool_t noError) {
dc1455ee 1125 // general method to normalize all vertical slices of a th2 to unity
1126 // i.e. get a probability matrix
1127 if(!histo) {
1128 printf(" > NormalizeTH2D:: fatal error, NULL pointer passed < \n");
1129 return NULL;
1130 }
1131 Int_t binsX = histo->GetXaxis()->GetNbins();
1132 Int_t binsY = histo->GetYaxis()->GetNbins();
1133
1134 // normalize all slices in x
1135 for(Int_t i(0); i < binsX; i++) { // for each vertical slice
1136 Double_t weight = 0;
1137 for(Int_t j(0); j < binsY; j++) { // loop over all the horizontal components
1138 weight+=histo->GetBinContent(i+1, j+1);
1139 } // now we know the total weight
1140 for(Int_t j(0); j < binsY; j++) {
1141 if (weight <= 0 ) continue;
1142 histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
512ced40
RAB		 1143 if(noError) histo->SetBinError( 1+i, j+1, 0.);
1144 else histo->SetBinError( 1+i, j+1, histo->GetBinError( 1+i, j+1)/weight);
dc1455ee 1145 }
1146 }
1147 return histo;
1148}
1149//_____________________________________________________________________________
53547ff2 1150TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
dc1455ee 1151 // return a TH1D with the supplied histogram rebinned to the supplied bins
53547ff2 1152 // the returned histogram is new, the original is deleted from the heap if kill is true
dc1455ee 1153 if(!histo || !bins) {
53547ff2 1154 printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
dc1455ee 1155 return NULL;
1156 }
1157 // create the output histo
1158 TString name = histo->GetName();
1159 name+="_template";
1160 name+=suffix;
1161 TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1162 for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
4292ca60 1163 // loop over the bins of the old histo and fill the new one with its data
1164 rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
dc1455ee 1165 }
53547ff2 1166 if(kill) delete histo;
dc1455ee 1167 return rebinned;
1168}
1169//_____________________________________________________________________________
4292ca60 1170TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
53547ff2 1171 // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
4292ca60 1172 if(!fResponseMaker || !binsTrue || !binsRec) {
1173 printf(" > RebinTH2D:: function called with NULL arguments < \n");
1174 return 0x0;
dc1455ee 1175 }
4292ca60 1176 TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1177 return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
dc1455ee 1178}
1179//_____________________________________________________________________________
d7ec324f 1180TH2D* AliJetFlowTools::MatrixMultiplication(TH2D* a, TH2D* b, TString name)
1181{
1182 // multiply two matrices
1183 if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1184 TH2D* c = (TH2D*)a->Clone("c");
1185 for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1186 for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1187 Double_t val = 0;
1188 for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1189 Int_t y2 = x1;
1190 val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1191 }
1192 c->SetBinContent(x2, y1, val);
512ced40 1193 c->SetBinError(x2, y1, 0.);
dc1455ee 1194 }
1195 }
d7ec324f 1196 if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1197 return c;
dc1455ee 1198}
1199//_____________________________________________________________________________
d7ec324f 1200TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale)
1201{
dc1455ee 1202 // normalize a th1d to a certain scale
4292ca60 1203 histo->Sumw2();
1204 Double_t integral = histo->Integral()*scale;
1205 if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1206 else if (scale != 1.) histo->Scale(1./scale, "width");
1207 else printf(" > Histogram integral < 0, cannot normalize \n");
1208 return histo;
dc1455ee 1209}
1210//_____________________________________________________________________________
51e6bc5a 1211TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
dc1455ee 1212{
1213 // Calculate pearson coefficients from covariance matrix
51e6bc5a 1214 TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1215 Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
ad04a83c 1216 Double_t pearson(0.);
dc1455ee 1217 if(nrows==0 && ncols==0) return 0x0;
ef12d5a5 1218 for(Int_t row = 0; row < nrows; row++) {
1219 for(Int_t col = 0; col<ncols; col++) {
51e6bc5a 1220 if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1221 (*pearsonCoefficients)(row,col) = pearson;
dc1455ee 1222 }
1223 }
51e6bc5a 1224 return pearsonCoefficients;
dc1455ee 1225}
1226//_____________________________________________________________________________
549b5f40 1227TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
d7ec324f 1228 // smoothen the spectrum using a user defined function
1229 // returns a clone of the original spectrum if fitting failed
1230 // if kill is kTRUE the input spectrum will be deleted from the heap
1231 // if 'count' is selected, bins are filled with integers (necessary if the
1232 // histogram is interpreted in a routine which accepts only counts)
549b5f40
RAB		 1233 if(!spectrum || !function) return 0x0;
1234 if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
d7ec324f 1235 TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1236 temp->Sumw2(); // if already called on the original, this will give off a warning but do nothing
549b5f40 1237 TFitResultPtr r = temp->Fit(function, "WLIS", "", min, max);
d7ec324f 1238 if((int)r == 0) { // MINUIT status
1239 for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1240 if(temp->GetBinCenter(i) > start) { // from this pt value use extrapolation
1241 if(counts) temp->SetBinContent(i, (int)function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
1242 else temp->SetBinContent(i, function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
1243 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1244 }
1245 }
1246 }
1247 if(kill) delete spectrum;
1248 return temp;
1249}
1250//_____________________________________________________________________________
1251void AliJetFlowTools::Style(TCanvas* c, TString style)
1252{
1253 // set a default style for a canvas
1254 if(!strcmp(style.Data(), "PEARSON")) {
1255 printf(" > style PEARSON canvas < \n");
1256 gStyle->SetOptStat(0);
1257 c->SetGridx();
1258 c->SetGridy();
1259 c->SetTicks();
1260 return;
1261 } else if(!strcmp(style.Data(), "SPECTRUM")) {
1262 printf(" > style SPECTRUM canvas < \n");
1263 gStyle->SetOptStat(0);
1264 c->SetLogy();
1265 c->SetGridx();
1266 c->SetGridy();
1267 c->SetTicks();
1268 return;
1269 } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1270}
1271//_____________________________________________________________________________
1272void AliJetFlowTools::Style(TVirtualPad* c, TString style)
1273{
1274 // set a default style for a canvas
1275 if(!strcmp(style.Data(), "PEARSON")) {
1276 printf(" > style PEARSON pad < \n");
1277 gStyle->SetOptStat(0);
1278 c->SetGridx();
1279 c->SetGridy();
1280 c->SetTicks();
1281 return;
1282 } else if(!strcmp(style.Data(), "SPECTRUM")) {
1283 printf(" > style SPECTRUM pad < \n");
1284 gStyle->SetOptStat(0);
1285 c->SetLogy();
1286 c->SetGridx();
1287 c->SetGridy();
1288 c->SetTicks();
1289 return;
1290 } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1291}
1292//_____________________________________________________________________________
53547ff2
RAB		 1293void AliJetFlowTools::Style(TLegend* l)
1294{
1295 // set a default style for a legend
1296 l->SetTextSize(.06);
1297 l->SetFillColor(0);
1298}
1299//_____________________________________________________________________________
1300void AliJetFlowTools::Style(TH1* h, EColor col, histoType type)
1301{
1302 // style a histo
1303 h->SetLineColor(col);
1304 h->SetMarkerColor(col);
1305 h->SetLineWidth(2.);
1306 h->SetMarkerSize(2.);
1307 h->SetTitleSize(0.06);
1308 h->GetYaxis()->SetLabelSize(0.06);
1309 h->GetXaxis()->SetLabelSize(0.06);
1310 h->GetYaxis()->SetTitleSize(0.06);
1311 h->GetXaxis()->SetTitleSize(0.06);
1312 h->GetYaxis()->SetTitleOffset(1.);
1313 h->GetXaxis()->SetTitleOffset(.9);
1314 switch (type) {
1315 case kInPlaneSpectrum : {
1316 h->SetTitle("IN PLANE");
1317 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1318 h->GetYaxis()->SetTitle("counts");
1319 } break;
1320 case kOutPlaneSpectrum : {
1321 h->SetTitle("OUT OF PLANE");
1322 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1323 h->GetYaxis()->SetTitle("counts");
1324 } break;
1325 case kUnfoldedSpectrum : {
1326 h->SetTitle("UNFOLDED");
1327 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1328 h->GetYaxis()->SetTitle("counts");
1329 } break;
1330 case kFoldedSpectrum : {
1331 h->SetTitle("FOLDED");
1332 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1333 h->GetYaxis()->SetTitle("counts");
1334 } break;
1335 case kMeasuredSpectrum : {
1336 h->SetTitle("MEASURED");
1337 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1338 h->GetYaxis()->SetTitle("counts");
1339 } break;
1340 default : break;
1341 }
1342}
1343//_____________________________________________________________________________
1344void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type)
1345{
1346 // style a histo
1347 h->SetLineColor(col);
1348 h->SetMarkerColor(col);
1349 h->SetLineWidth(2.);
1350 h->SetMarkerSize(2.);
1351 h->GetYaxis()->SetLabelSize(0.06);
1352 h->GetXaxis()->SetLabelSize(0.06);
1353 h->GetYaxis()->SetTitleSize(0.06);
1354 h->GetXaxis()->SetTitleSize(0.06);
1355 h->GetYaxis()->SetTitleOffset(1.);
1356 h->GetXaxis()->SetTitleOffset(.9);
1357 switch (type) {
1358 case kInPlaneSpectrum : {
1359 h->SetTitle("IN PLANE");
1360 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1361 h->GetYaxis()->SetTitle("counts");
1362 } break;
1363 case kOutPlaneSpectrum : {
1364 h->SetTitle("OUT OF PLANE");
1365 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1366 h->GetYaxis()->SetTitle("counts");
1367 } break;
1368 case kUnfoldedSpectrum : {
1369 h->SetTitle("UNFOLDED");
1370 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1371 h->GetYaxis()->SetTitle("counts");
1372 } break;
1373 case kFoldedSpectrum : {
1374 h->SetTitle("FOLDED");
1375 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1376 h->GetYaxis()->SetTitle("counts");
1377 } break;
1378 case kMeasuredSpectrum : {
1379 h->SetTitle("MEASURED");
1380 h->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1381 h->GetYaxis()->SetTitle("counts");
1382 } break;
1383 default : break;
1384 }
1385}
1386//_____________________________________________________________________________
549b5f40 1387void AliJetFlowTools::PostProcess(TString def, TString in, TString out, Int_t columns) const
d7ec324f 1388{
1389 // go through the output file and perform post processing routines
1390 // can either be performed in one go with the unfolding, or at a later stage
53547ff2 1391 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
d7ec324f 1392 TFile readMe(in.Data(), "READ"); // open file read-only
1393 if(readMe.IsZombie()) {
1394 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1395 return;
1396 }
1397 printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
1398 readMe.ls();
1399 TList* listOfKeys((TList*)readMe.GetListOfKeys());
1400 if(!listOfKeys) {
1401 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
1402 return;
1403 }
1404 // prepare necessary canvasses
53547ff2
RAB		 1405 TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
1406 TCanvas* canvasOut(new TCanvas("PearsonOut", "PearsonOut"));
1407 TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
1408 TCanvas* canvasRatioMeasuredRefoldedOut(new TCanvas("RefoldedOut", "RefoldedOut"));
1409 TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
1410 TCanvas* canvasSpectraOut(new TCanvas("SpectraOut", "SpectraOut"));
1411 TCanvas* canvasRatio(new TCanvas("Ratio", "Ratio"));
1412 TCanvas* canvasV2(new TCanvas("V2", "V2"));
1413 TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
1414 TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
1415 TCanvas* canvasMasterOut(new TCanvas("defaultOut", "defaultOut"));
d7ec324f 1416 canvasMISC->Divide(4, 2);
1417 TDirectoryFile* defDir(0x0);
1418
1419 // get an estimate of the number of outputs and find the default set
1420 Int_t cacheMe(0);
1421 for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
1422 if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
1423 if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1424 cacheMe++;
1425 }
1426 }
53547ff2
RAB		 1427 Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%4)>0));
1428 canvasIn->Divide(columns, rows);
1429 canvasOut->Divide(columns, rows);
1430 canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
1431 canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
1432 canvasSpectraIn->Divide(columns, rows);
1433 canvasSpectraOut->Divide(columns, rows);
1434 canvasRatio->Divide(columns, rows);
1435 canvasV2->Divide(columns, rows);
d7ec324f 1436
53547ff2
RAB		 1437 canvasMasterIn->Divide(columns, rows);
1438 canvasMasterOut->Divide(columns, rows);
d7ec324f 1439 // extract the default output
549b5f40
RAB		 1440 TH1D* deunfoldedJetSpectrumIn(0x0);
1441 TH1D* deunfoldedJetSpectrumOut(0x0);
d7ec324f 1442 THStack stackIn("StackRatioIn","StackRatioIn");
1443 THStack stackOut("StackRatioOut", "StackRatioOut");
1444 if(defDir) {
1445 TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
1446 TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
549b5f40
RAB		 1447 if(defDirIn) deunfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
1448 if(deunfoldedJetSpectrumIn) stackIn.Add(deunfoldedJetSpectrumIn);
1449 if(defDirOut) deunfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
1450 if(deunfoldedJetSpectrumOut) stackOut.Add(deunfoldedJetSpectrumOut);
d7ec324f 1451 printf(" > succesfully extracted default results < \n");
1452 }
1453
1454 // loop through the directories, only plot the graphs if the deconvolution converged
1455 TDirectoryFile* tempDir(0x0);
1456 TDirectoryFile* tempIn(0x0);
1457 TDirectoryFile* tempOut(0x0);
1458 for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
1459 tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
1460 if(!tempDir) continue;
1461 TString dirName(tempDir->GetName());
1462 tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
1463 tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
1464 j++;
1465 if(tempIn) {
1466 // to see if the unfolding converged try to extract the pearson coefficients
1467 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
1468 if(pIn) {
1469 printf(" - %s in plane converged \n", dirName.Data());
1470 canvasIn->cd(j);
1471 Style(gPad, "PEARSON");
1472 pIn->DrawCopy("colz");
1473 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
1474 if(rIn) {
53547ff2 1475 Style(rIn);
d7ec324f 1476 printf(" > found RatioRefoldedMeasured < \n");
1477 canvasRatioMeasuredRefoldedIn->cd(j);
53547ff2 1478 rIn->Draw("ac");
d7ec324f 1479 }
1480 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
1481 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
1482 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
1483 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
1484 if(dvector && avalue && rm && eff) {
53547ff2
RAB		 1485 Style(dvector);
1486 Style(avalue);
1487 Style(rm);
1488 Style(eff);
d7ec324f 1489 canvasMISC->cd(1);
1490 Style(gPad, "SPECTRUM");
1491 dvector->DrawCopy();
1492 canvasMISC->cd(2);
1493 Style(gPad, "SPECTRUM");
1494 avalue->DrawCopy();
1495 canvasMISC->cd(3);
1496 Style(gPad, "PEARSON");
1497 rm->DrawCopy("colz");
1498 canvasMISC->cd(4);
1499 eff->DrawCopy();
53547ff2
RAB		 1500 } else if(rm && eff) {
1501 Style(rm);
1502 Style(eff);
1503 canvasMISC->cd(3);
1504 Style(gPad, "PEARSON");
1505 rm->DrawCopy("colz");
1506 canvasMISC->cd(4);
1507 eff->DrawCopy();
d7ec324f 1508 }
1509 }
1510 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
1511 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
1512 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
1513 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
549b5f40
RAB		 1514 if(deunfoldedJetSpectrumIn) {
1515 TH1D* temp((TH1D*)deunfoldedJetSpectrumIn->Clone(Form("deunfoldedJetSpectrumIn_%s", dirName.Data())));
d7ec324f 1516 temp->Divide(unfoldedSpectrum);
1517 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
1518 temp->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1519 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
1520 canvasMasterIn->cd(j);
1521 temp->GetXaxis()->SetRangeUser(0., 2);
1522 temp->DrawCopy();
1523 }
1524 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
1525 canvasSpectraIn->cd(j);
1526 Style(gPad);
53547ff2 1527 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
d7ec324f 1528 unfoldedSpectrum->DrawCopy();
53547ff2 1529 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
d7ec324f 1530 inputSpectrum->DrawCopy("same");
53547ff2 1531 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
d7ec324f 1532 refoldedSpectrum->DrawCopy("same");
1533 TLegend* l(AddLegend(gPad));
53547ff2
RAB		 1534 Style(l);
1535 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1536 Float_t chi(fitStatus->GetBinContent(1));
1537 Float_t pen(fitStatus->GetBinContent(2));
d7ec324f 1538 Int_t dof((int)fitStatus->GetBinContent(3));
53547ff2
RAB		 1539 Float_t beta(fitStatus->GetBinContent(4));
1540 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1541 } else if (fitStatus) { // only available in SVD fit
1542 Int_t reg((int)fitStatus->GetBinContent(1));
1543 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
d7ec324f 1544 }
1545 }
1546 }
1547 if(tempOut) {
1548 TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
1549 if(pOut) {
1550 printf(" - %s out of plane converged \n", dirName.Data());
1551 canvasOut->cd(j);
1552 Style(gPad, "PEARSON");
1553 pOut->DrawCopy("colz");
1554 TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
1555 if(rOut) {
53547ff2 1556 Style(rOut);
d7ec324f 1557 printf(" > found RatioRefoldedMeasured < \n");
1558 canvasRatioMeasuredRefoldedOut->cd(j);
53547ff2 1559 rOut->Draw("ac");
d7ec324f 1560 }
1561 TH1D* dvector((TH1D*)tempOut->Get("dVector"));
1562 TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
1563 TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
1564 TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
1565 if(dvector && avalue && rm && eff) {
53547ff2
RAB		 1566 Style(dvector);
1567 Style(avalue);
1568 Style(rm);
1569 Style(eff);
d7ec324f 1570 canvasMISC->cd(5);
1571 Style(gPad, "SPECTRUM");
1572 dvector->DrawCopy();
1573 canvasMISC->cd(6);
1574 Style(gPad, "SPECTRUM");
1575 avalue->DrawCopy();
1576 canvasMISC->cd(7);
1577 Style(gPad, "PEARSON");
1578 rm->DrawCopy("colz");
1579 canvasMISC->cd(8);
1580 eff->DrawCopy();
53547ff2
RAB		 1581 } else if(rm && eff) {
1582 Style(rm);
1583 Style(eff);
549b5f40 1584 canvasMISC->cd(7);
53547ff2
RAB		 1585 Style(gPad, "PEARSON");
1586 rm->DrawCopy("colz");
549b5f40 1587 canvasMISC->cd(8);
53547ff2 1588 eff->DrawCopy();
d7ec324f 1589 }
1590 }
1591 TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
1592 TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
1593 TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
1594 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
549b5f40
RAB		 1595 if(deunfoldedJetSpectrumOut) {
1596 TH1D* temp((TH1D*)deunfoldedJetSpectrumOut->Clone(Form("deunfoldedJetSpectrumOut_%s", dirName.Data())));
d7ec324f 1597 temp->Divide(unfoldedSpectrum);
1598 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
1599 temp->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1600 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
1601 canvasMasterOut->cd(j);
1602 temp->GetXaxis()->SetRangeUser(0., 2.);
1603 temp->DrawCopy();
1604 }
1605 TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
1606 canvasSpectraOut->cd(j);
1607 Style(gPad);
53547ff2 1608 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
d7ec324f 1609 unfoldedSpectrum->DrawCopy();
53547ff2 1610 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
d7ec324f 1611 inputSpectrum->DrawCopy("same");
53547ff2 1612 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
d7ec324f 1613 refoldedSpectrum->DrawCopy("same");
1614 TLegend* l(AddLegend(gPad));
53547ff2
RAB		 1615 Style(l);
1616 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
1617 Float_t chi(fitStatus->GetBinContent(1));
1618 Float_t pen(fitStatus->GetBinContent(2));
1619 Int_t dof((int)fitStatus->GetBinContent(3));
1620 Float_t beta(fitStatus->GetBinContent(4));
1621 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
1622 } else if (fitStatus) { // only available in SVD fit
1623 Int_t reg((int)fitStatus->GetBinContent(1));
1624 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
d7ec324f 1625 }
1626 }
1627 }
1628 canvasRatio->cd(j);
1629 TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
1630 if(ratioYield) {
53547ff2 1631 Style(ratioYield);
549b5f40 1632 ratioYield->GetYaxis()->SetRangeUser(-1., 3.);
53547ff2 1633 ratioYield->Draw("ac");
d7ec324f 1634 }
1635 canvasV2->cd(j);
1636 TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
1637 if(ratioV2) {
53547ff2 1638 Style(ratioV2);
549b5f40 1639 ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
53547ff2 1640 ratioV2->Draw("ac");
d7ec324f 1641 }
1642 }
1643 TFile output(out.Data(), "RECREATE");
512ced40 1644 SavePadToPDF(canvasIn);
d7ec324f 1645 canvasIn->Write();
512ced40 1646 SavePadToPDF(canvasOut);
d7ec324f 1647 canvasOut->Write();
512ced40 1648 SavePadToPDF(canvasRatioMeasuredRefoldedIn);
d7ec324f 1649 canvasRatioMeasuredRefoldedIn->Write();
512ced40 1650 SavePadToPDF(canvasRatioMeasuredRefoldedOut);
d7ec324f 1651 canvasRatioMeasuredRefoldedOut->Write();
512ced40 1652 SavePadToPDF(canvasSpectraIn);
d7ec324f 1653 canvasSpectraIn->Write();
512ced40 1654 SavePadToPDF(canvasSpectraOut);
d7ec324f 1655 canvasSpectraOut->Write();
512ced40 1656 SavePadToPDF(canvasRatio);
d7ec324f 1657 canvasRatio->Write();
512ced40 1658 SavePadToPDF(canvasV2);
d7ec324f 1659 canvasV2->Write();
512ced40 1660 SavePadToPDF(canvasMasterIn);
d7ec324f 1661 canvasMasterIn->Write();
512ced40 1662 SavePadToPDF(canvasMasterIn);
d7ec324f 1663 canvasMasterOut->Write();
512ced40 1664 SavePadToPDF(canvasMISC);
d7ec324f 1665 canvasMISC->Write();
1666 output.Write();
1667 output.Close();
1668}
1669//_____________________________________________________________________________
4292ca60 1670Bool_t AliJetFlowTools::SetRawInput (
1671 TH2D* detectorResponse, // detector response matrix
1672 TH1D* jetPtIn, // in plane jet spectrum
1673 TH1D* jetPtOut, // out of plane jet spectrum
1674 TH1D* dptIn, // in plane delta pt distribution
1675 TH1D* dptOut, // out of plane delta pt distribution
1676 Int_t eventCount) {
1677 // set input histograms manually
1678 fDetectorResponse = detectorResponse;
1679 fSpectrumIn = jetPtIn;
1680 fSpectrumOut = jetPtOut;
1681 fDptInDist = dptIn;
1682 fDptOutDist = dptOut;
1683 fRawInputProvided = kTRUE;
1684 // check if all data is provided
1685 if(!fDetectorResponse) {
1686 printf(" fDetectorResponse not found \n ");
1687 return kFALSE;
1688 }
1689 // check if the pt bin for true and rec have been set
1690 if(!fBinsTrue || !fBinsRec) {
1691 printf(" No true or rec bins set, please set binning ! \n");
1692 return kFALSE;
1693 }
1694 if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
1695 // procedures, these profiles will be nonsensical, user is responsible
1696 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1697 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1698 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1699 }
1700 // normalize spectra to event count if requested
1701 if(fNormalizeSpectra) {
1702 fEventCount = eventCount;
1703 if(fEventCount > 0) {
1704 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
1705 fSpectrumOut->Sumw2();
1706 fSpectrumIn->Scale(1./((double)fEventCount));
1707 fSpectrumOut->Scale(1./((double)fEventCount));
1708 }
1709 }
1710 if(!fNormalizeSpectra && fEventCount > 0) {
1711 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
1712 fSpectrumOut->Sumw2();
1713 fSpectrumIn->Scale(1./((double)fEventCount));
1714 fSpectrumOut->Scale(1./((double)fEventCount));
1715 }
1716 fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
1717 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
1718 fDptIn->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1719 fDptIn->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1720 fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
1721 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
1722 fDptOut->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1723 fDptOut->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1724
1725 return kTRUE;
1726}
1727//_____________________________________________________________________________
1728TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
dc1455ee 1729{
ef12d5a5 1730 // return ratio of h1 / h2
1731 // histograms can have different binning. errors are propagated as uncorrelated
20abfcc4 1732 if(!(h1 && h2) ) {
1733 printf(" GetRatio called with NULL argument(s) \n ");
1734 return 0x0;
1735 }
ad04a83c 1736 Int_t j(0);
4292ca60 1737 TGraphErrors *gr = new TGraphErrors();
53547ff2 1738 gr->GetXaxis()->SetTitle("p_{T} [GeV/c]");
ef12d5a5 1739 Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
4292ca60 1740 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
dc1455ee 1741 binCent = h1->GetXaxis()->GetBinCenter(i);
4292ca60 1742 if(xmax > 0. && binCent > xmax) continue;
dc1455ee 1743 j = h2->FindBin(binCent);
1744 binWidth = h1->GetXaxis()->GetBinWidth(i);
ad04a83c 1745 if(h2->GetBinContent(j) > 0.) {
dc1455ee 1746 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
1747 Double_t A = 1./h2->GetBinContent(j)*h1->GetBinError(i);
4292ca60 1748 Double_t B = 0.;
1749 if(h2->GetBinError(j)>0.) {
1750 B = -1.*h1->GetBinContent(i)/(h2->GetBinContent(j)*h2->GetBinContent(j))*h2->GetBinError(j);
1751 error2 = A*A + B*B;
ef12d5a5 1752 } else error2 = A*A;
1753 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
dc1455ee 1754 gr->SetPoint(gr->GetN(),binCent,ratio);
ef12d5a5 1755 gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
dc1455ee 1756 }
1757 }
ad04a83c 1758 if(appendFit) {
4292ca60 1759 TF1* fit(new TF1("lin", "pol0", 10, 100));
ad04a83c 1760 gr->Fit(fit);
1761 }
4292ca60 1762 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
dc1455ee 1763 return gr;
1764}
1765//_____________________________________________________________________________
ef12d5a5 1766TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
1767{
1768 // get v2 from difference of in plane, out of plane yield
1769 // h1 must hold the in-plane yield, h2 holds the out of plane yield
549b5f40
RAB		 1770 // different binning is allowed but will mean that the error
1771 // propagation is unreliable
ef12d5a5 1772 // r is the event plane resolution for the chosen centrality
1773 if(!(h1 && h2) ) {
1774 printf(" GetV2 called with NULL argument(s) \n ");
1775 return 0x0;
1776 }
1777 Int_t j(0);
1778 TGraphErrors *gr = new TGraphErrors();
53547ff2 1779 gr->GetXaxis()->SetTitle("p_{T} [GeV/c]");
ef12d5a5 1780 Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
1781 Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
1782 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
1783 binCent = h1->GetXaxis()->GetBinCenter(i);
1784 j = h2->FindBin(binCent);
1785 binWidth = h1->GetXaxis()->GetBinWidth(i);
1786 if(h2->GetBinContent(j) > 0.) {
1787 in = h1->GetBinContent(i);
1788 ein = h1->GetBinError(i);
1789 out = h2->GetBinContent(j);
1790 eout = h2->GetBinError(j);
1791 ratio = pre*((in-out)/(in+out));
549b5f40 1792 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 1793 if(error2 > 0) error2 = TMath::Sqrt(error2);
1794 gr->SetPoint(gr->GetN(),binCent,ratio);
1795 gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
1796 }
1797 }
1798 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
1799 return gr;
1800}
1801//_____________________________________________________________________________
549b5f40
RAB		 1802void AliJetFlowTools::WriteObject(TObject* object, TString suffix, Bool_t kill) {
1803 // write object, if a unique identifier is given the object is cloned
1804 // and the clone is saved. setting kill to true will delete the original obect from the heap
4292ca60 1805 if(!object) {
1806 printf(" > WriteObject:: called with NULL arguments \n ");
1807 return;
549b5f40
RAB		 1808 } else if(!strcmp("", suffix.Data())) object->Write();
1809 else {
1810 TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
1811 newObject->Write();
1812 }
1813 if(kill) delete object;
4292ca60 1814}
1815//_____________________________________________________________________________
1816TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
1817 // construt a delta pt response matrix from supplied dpt distribution
1818 // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
1819 // do a weighted rebinning to a (coarser) dpt distribution
1820 // be careful with the binning of the dpt response: it should be equal to that
1821 // of the response matrix, otherwise dpt and response matrices cannot be multiplied
1822 //
1823 // the response matrix will be square and have the same binning
1824 // (min, max and granularity) of the input histogram
1825 Int_t bins(dpt->GetXaxis()->GetNbins()); // number of bins, will also be no of rows, columns
1826 Double_t _bins[bins+1]; // prepare array with bin borders
1827 for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
1828 _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1); // get upper edge
1829 TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
1830 for(Int_t j(0); j < bins+1; j++) { // loop on pt true slices j
1831 Bool_t skip = kFALSE;
1832 for(Int_t k(0); k < bins+1; k++) { // loop on pt gen slices k
1833 (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
1834 if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
1835 }
1836 }
1837 return res;
1838}
ef12d5a5 1839//_____________________________________________________________________________
1840TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1841 if(!binsTrue || !binsRec) {
1842 printf(" > GetUnityResponse:: function called with NULL arguments < \n");
1843 return 0x0;
1844 }
1845 TString name(Form("unityResponse_%s", suffix.Data()));
1846 TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
1847 for(Int_t i(0); i < binsTrue->GetSize(); i++) {
1848 for(Int_t j(0); j < binsRec->GetSize(); j++) {
1849 if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
1850 }
1851 }
1852 return unity;
1853}
4292ca60 1854//_____________________________________________________________________________
549b5f40 1855void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
4292ca60 1856 // save configuration parameters to histogram
549b5f40 1857 TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
ef12d5a5 1858 summary->SetBinContent(1, fBetaIn);
1859 summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
1860 summary->SetBinContent(2, fBetaOut);
1861 summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
1862 summary->SetBinContent(3, fCentralityBin);
1863 summary->GetXaxis()->SetBinLabel(3, "fCentralityBin");
1864 summary->SetBinContent(4, (int)convergedIn);
1865 summary->GetXaxis()->SetBinLabel(4, "convergedIn");
1866 summary->SetBinContent(5, (int)convergedOut);
1867 summary->GetXaxis()->SetBinLabel(5, "convergedOut");
1868 summary->SetBinContent(6, (int)fAvoidRoundingError);
1869 summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
1870 summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
1871 summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
1872 summary->SetBinContent(8, (int)fPrior);
1873 summary->GetXaxis()->SetBinLabel(8, "fPrior");
1874 summary->SetBinContent(9, fSVDRegIn);
1875 summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
1876 summary->SetBinContent(10, fSVDRegOut);
1877 summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
1878 summary->SetBinContent(11, (int)fSVDToy);
1879 summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
1880 summary->SetBinContent(12, fJetRadius);
1881 summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
1882 summary->SetBinContent(13, (int)fNormalizeSpectra);
1883 summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
549b5f40
RAB		 1884 summary->SetBinContent(14, (int)fSmoothenPrior);
1885 summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
ef12d5a5 1886 summary->SetBinContent(15, (int)fTestMode);
1887 summary->GetXaxis()->SetBinLabel(15, "fTestMode");
1888 summary->SetBinContent(16, (int)fUseDetectorResponse);
1889 summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
549b5f40
RAB		 1890 summary->SetBinContent(17, fBayesianIterIn);
1891 summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
1892 summary->SetBinContent(18, fBayesianIterOut);
1893 summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
1894 summary->SetBinContent(19, fBayesianSmoothIn);
1895 summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
1896 summary->SetBinContent(20, fBayesianSmoothOut);
1897 summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
4292ca60 1898}
1899//_____________________________________________________________________________
1900void AliJetFlowTools::ResetAliUnfolding() {
1901 // ugly function: reset all unfolding parameters
1902 TVirtualFitter* fitter(TVirtualFitter::GetFitter());
1903 if(fitter) {
1904 printf(" > Found fitter, will delete it < \n");
1905 delete fitter;
1906 }
d7ec324f 1907 if(gMinuit) {
1908 printf(" > Found gMinuit, will re-create it < \n");
1909 delete gMinuit;
1910 gMinuit = new TMinuit;
1911 }
4292ca60 1912 AliUnfolding::fgCorrelationMatrix = 0;
1913 AliUnfolding::fgCorrelationMatrixSquared = 0;
1914 AliUnfolding::fgCorrelationCovarianceMatrix = 0;
1915 AliUnfolding::fgCurrentESDVector = 0;
1916 AliUnfolding::fgEntropyAPriori = 0;
1917 AliUnfolding::fgEfficiency = 0;
1918 AliUnfolding::fgUnfoldedAxis = 0;
1919 AliUnfolding::fgMeasuredAxis = 0;
1920 AliUnfolding::fgFitFunction = 0;
1921 AliUnfolding::fgMaxInput = -1;
1922 AliUnfolding::fgMaxParams = -1;
1923 AliUnfolding::fgOverflowBinLimit = -1;
1924 AliUnfolding::fgRegularizationWeight = 10000;
1925 AliUnfolding::fgSkipBinsBegin = 0;
1926 AliUnfolding::fgMinuitStepSize = 0.1;
1927 AliUnfolding::fgMinuitPrecision = 1e-6;
1928 AliUnfolding::fgMinuitMaxIterations = 1000000;
1929 AliUnfolding::fgMinuitStrategy = 1.;
1930 AliUnfolding::fgMinimumInitialValue = kFALSE;
1931 AliUnfolding::fgMinimumInitialValueFix = -1;
1932 AliUnfolding::fgNormalizeInput = kFALSE;
1933 AliUnfolding::fgNotFoundEvents = 0;
1934 AliUnfolding::fgSkipBin0InChi2 = kFALSE;
1935 AliUnfolding::fgBayesianSmoothing = 1;
1936 AliUnfolding::fgBayesianIterations = 10;
1937 AliUnfolding::fgDebug = kFALSE;
1938 AliUnfolding::fgCallCount = 0;
1939 AliUnfolding::fgPowern = 5;
1940 AliUnfolding::fChi2FromFit = 0.;
1941 AliUnfolding::fPenaltyVal = 0.;
1942 AliUnfolding::fAvgResidual = 0.;
1943 AliUnfolding::fgPrintChi2Details = 0;
1944 AliUnfolding::fgCanvas = 0;
1945 AliUnfolding::fghUnfolded = 0;
1946 AliUnfolding::fghCorrelation = 0;
1947 AliUnfolding::fghEfficiency = 0;
1948 AliUnfolding::fghMeasured = 0;
1949 AliUnfolding::SetMinuitStepSize(1.);
1950 AliUnfolding::SetMinuitPrecision(1e-6);
1951 AliUnfolding::SetMinuitMaxIterations(100000);
1952 AliUnfolding::SetMinuitStrategy(2.);
1953 AliUnfolding::SetDebug(1);
1954}
1955//_____________________________________________________________________________
ef12d5a5 1956TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) {
4292ca60 1957 // protect heap by adding unique qualifier to name
1958 if(!protect) return 0x0;
1959 TH1D* p = (TH1D*)protect->Clone();
ef12d5a5 1960 TString tempString(fActiveString);
1961 tempString+=suffix;
1962 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1963 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 1964 if(kill) delete protect;
1965 return p;
1966}
1967//_____________________________________________________________________________
ef12d5a5 1968TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) {
4292ca60 1969 // protect heap by adding unique qualifier to name
1970 if(!protect) return 0x0;
1971 TH2D* p = (TH2D*)protect->Clone();
ef12d5a5 1972 TString tempString(fActiveString);
1973 tempString+=suffix;
1974 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1975 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 1976 if(kill) delete protect;
1977 return p;
1978}
1979//_____________________________________________________________________________
ef12d5a5 1980TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) {
4292ca60 1981 // protect heap by adding unique qualifier to name
1982 if(!protect) return 0x0;
1983 TGraphErrors* p = (TGraphErrors*)protect->Clone();
ef12d5a5 1984 TString tempString(fActiveString);
1985 tempString+=suffix;
1986 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1987 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 1988 if(kill) delete protect;
1989 return p;
1990}
1991//_____________________________________________________________________________