]> 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
unfolding bugfixes and updates
[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)
18//
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"
4292ca60 44#include "TGraphErrors.h"
dc1455ee 45#include "TArrayD.h"
46#include "TList.h"
47#include "TMinuit.h"
48#include "TVirtualFitter.h"
49#include "TLegend.h"
50#include "TCanvas.h"
51#include "TStyle.h"
dc1455ee 52#include "TLine.h"
ad04a83c 53#include "TMath.h"
4292ca60 54#include "TVirtualFitter.h"
ef12d5a5 55#include "TFitResultPtr.h"
4292ca60 56// aliroot includes
dc1455ee 57#include "AliUnfolding.h"
58#include "AliAnaChargedJetResponseMaker.h"
4292ca60 59// class includes
dc1455ee 60#include "AliJetFlowTools.h"
51e6bc5a 61// roo unfold includes (make sure you have these available on your system)
62#include "RooUnfold.h"
63#include "RooUnfoldResponse.h"
64#include "RooUnfoldSvd.h"
65#include "TSVDUnfold.h"
dc1455ee 66
67using namespace std;
dc1455ee 68//_____________________________________________________________________________
69AliJetFlowTools::AliJetFlowTools() :
4292ca60 70 fResponseMaker (new AliAnaChargedJetResponseMaker()),
ef12d5a5 71 fPower (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
4292ca60 72 fSaveFull (kFALSE),
dc1455ee 73 fActiveString (""),
ad04a83c 74 fActiveDir (0x0),
dc1455ee 75 fInputList (0x0),
4292ca60 76 fRefreshInput (kTRUE),
dc1455ee 77 fOutputFileName ("UnfoldedSpectra.root"),
ad04a83c 78 fOutputFile (0x0),
dc1455ee 79 fCentralityBin (0),
dc1455ee 80 fDetectorResponse (0x0),
4292ca60 81 fBetaIn (.1),
82 fBetaOut (.1),
ef12d5a5 83 fAvoidRoundingError (kFALSE),
51e6bc5a 84 fUnfoldingAlgorithm (kChi2),
85 fPrior (kPriorMeasured),
dc1455ee 86 fBinsTrue (0x0),
87 fBinsRec (0x0),
ef12d5a5 88 fBinsTruePrior (0x0),
89 fBinsRecPrior (0x0),
4292ca60 90 fSVDRegIn (5),
91 fSVDRegOut (5),
51e6bc5a 92 fSVDToy (kTRUE),
93 fJetRadius (0.3),
20abfcc4 94 fEventCount (-1),
95 fNormalizeSpectra (kTRUE),
4292ca60 96 fSmoothenSpectrum (kTRUE),
97 fFitMin (60.),
98 fFitMax (105.),
99 fFitStart (75.),
ef12d5a5 100 fTestMode (kFALSE),
101 fNoDphi (kFALSE),
4292ca60 102 fRawInputProvided (kFALSE),
ef12d5a5 103 fEventPlaneRes (.63),
104 fUseDetectorResponse(kTRUE),
105 fTrainPower (kTRUE),
4292ca60 106 fRMSSpectrumIn (0x0),
107 fRMSSpectrumOut (0x0),
108 fRMSRatio (0x0),
ef12d5a5 109 fRMSV2 (0x0),
ad04a83c 110 fDeltaPtDeltaPhi (0x0),
111 fJetPtDeltaPhi (0x0),
dc1455ee 112 fSpectrumIn (0x0),
113 fSpectrumOut (0x0),
114 fDptInDist (0x0),
115 fDptOutDist (0x0),
116 fDptIn (0x0),
117 fDptOut (0x0),
118 fFullResponseIn (0x0),
119 fFullResponseOut (0x0),
120 fUnfoldedIn (0x0),
4292ca60 121 fUnfoldedOut (0x0) { // class constructor
122 // create response maker weight function
123 fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0 ,200));
124}
dc1455ee 125//_____________________________________________________________________________
126void AliJetFlowTools::Make() {
ad04a83c 127 // core function of the class
4292ca60 128 // 1) rebin the raw output of the jet task to the desired binnings
ad04a83c 129 // 2) calls the unfolding routine
130 // 3) writes output to file
4292ca60 131 // can be repeated multiple times with different configurations
132
ad04a83c 133 // 1) manipulation of input histograms
dc1455ee 134 // check if the input variables are present
4292ca60 135 if(fRefreshInput) {
136 if(!PrepareForUnfolding()) {
137 printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
138 return;
139 }
dc1455ee 140 }
4292ca60 141 // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
142 // parts of the spectrum can end up in over or underflow bins
143 TH1D* resizedJetPtIn = GetUnfoldingTemplate(fSpectrumIn, fBinsRec, TString("resized_in_"));
144 TH1D* resizedJetPtOut = GetUnfoldingTemplate(fSpectrumOut, fBinsRec, TString("resized_out_"));
145
146 // 1b) get the unfolding template
147 // the template will be used as a prior for the chi2 unfolding
148 // it holds thie rec spectrum, but is rebinned to the gen binning scheme
dc1455ee 149 TH1D* unfoldingTemplateIn = GetUnfoldingTemplate(fSpectrumIn, fBinsTrue, TString("in"));
150 TH1D* unfoldingTemplateOut = GetUnfoldingTemplate(fSpectrumOut, fBinsTrue, TString("out"));
151
4292ca60 152 // get the full response matrix from the dpt and the detector response
dc1455ee 153 fDetectorResponse = NormalizeTH2D(fDetectorResponse);
ef12d5a5 154 // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
155 // so that unfolding should return the initial spectrum
156 if(!fTestMode) {
157 fFullResponseIn = (fUseDetectorResponse) ? MatrixMultiplicationTH2D(fDptIn, fDetectorResponse) : fDptIn;
158 fFullResponseOut = (fUseDetectorResponse) ? MatrixMultiplicationTH2D(fDptOut, fDetectorResponse) : fDptOut;
159 } else {
160 fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
161 fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
162 }
4292ca60 163 // normalize each slide of the response to one
dc1455ee 164 NormalizeTH2D(fFullResponseIn);
165 NormalizeTH2D(fFullResponseOut);
4292ca60 166 // resize to desired binning scheme
167 TH2D* resizedResonseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
168 TH2D* resizedResonseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
169 // get the kinematic efficiency
dc1455ee 170 TH1D* kinematicEfficiencyIn = resizedResonseIn->ProjectionX();
4292ca60 171 kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
dc1455ee 172 TH1D* kinematicEfficiencyOut = resizedResonseOut->ProjectionX();
4292ca60 173 kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
174 // suppress the errors
51e6bc5a 175 for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
176 kinematicEfficiencyIn->SetBinError(1+i, 0.);
4292ca60 177 kinematicEfficiencyOut->SetBinError(1+i, 0.);
51e6bc5a 178 }
ad04a83c 179 // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
180 fActiveDir->cd(); // select active dir
181 TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
182 dirIn->cd(); // select inplane subdir
51e6bc5a 183 Bool_t convergedIn(kFALSE), convergedOut(kFALSE);
184 // select the unfolding method
185 switch (fUnfoldingAlgorithm) {
186 case kChi2 : {
187 convergedIn = UnfoldSpectrumChi2( // do the inplane unfolding
188 resizedJetPtIn,
189 resizedResonseIn,
190 kinematicEfficiencyIn,
191 unfoldingTemplateIn,
192 fUnfoldedIn,
193 TString("in"));
20abfcc4 194 printf(" > Spectrum (in plane) unfolded using kChi2 unfolding < \n");
51e6bc5a 195 } break;
196 case kSVD : {
197 convergedIn = UnfoldSpectrumSVD( // do the inplane unfolding
198 resizedJetPtIn,
199 resizedResonseIn,
200 kinematicEfficiencyIn,
201 unfoldingTemplateIn,
4292ca60 202 fUnfoldedIn,
51e6bc5a 203 TString("in"));
20abfcc4 204 printf(" > Spectrum (in plane) unfolded using kSVD unfolding < \n");
51e6bc5a 205 } break;
ef12d5a5 206 case kNone : { // do nothing, just rebin and optionally smooothen the spectrum
207 resizedResonseIn->SetNameTitle("measuredSpectrumIn", "measured spectrum, in plane");
208 if(fSmoothenSpectrum) { // see if we want to smooothen the spectrum
209 TH1D* resizedJetPtInJagged((TH1D*)resizedJetPtIn->Clone("cachedRawJetJaggedIn"));
210 resizedJetPtInJagged->SetNameTitle("measured spectrum before smoothening in", "measured spectrum, before smoothening in");
211 resizedJetPtInJagged = ProtectHeap(resizedJetPtInJagged);
212 resizedJetPtInJagged->Write(); // save the original
213 resizedJetPtIn->Sumw2();
214 TFitResultPtr r = resizedJetPtIn->Fit(fPower, "QWILS", "", fFitMin, fFitMax);
215 if((int)r == 0 ) {
216 for(Int_t i(0); i < resizedJetPtIn->GetNbinsX() + 1; i++) {
217 if(resizedJetPtIn->GetBinCenter(i) > fFitStart) { // from this pt value use extrapolation
218 resizedJetPtIn->SetBinContent(i,fPower->Integral(resizedJetPtIn->GetXaxis()->GetBinLowEdge(i),resizedJetPtIn->GetXaxis()->GetBinUpEdge(i))/resizedJetPtIn->GetXaxis()->GetBinWidth(i));
219 }
220 }
221 } else printf(" > PANIC, SMOOTHENING FAILED < \n");
222 }
223 fUnfoldedIn = ProtectHeap(resizedJetPtIn, kTRUE, TString("in"));
224 convergedIn = kTRUE;
225 } break;
226
51e6bc5a 227 default : {
228 printf(" > Selected unfolding method is not implemented yet ! \n");
229 return;
230 }
231 }
4292ca60 232 resizedResonseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
233 resizedResonseIn->SetXTitle("p_{T}^{true} [GeV/c]");
234 resizedResonseIn->SetYTitle("p_{T}^{rec} [GeV/c]");
235 resizedResonseIn = ProtectHeap(resizedResonseIn);
ad04a83c 236 resizedResonseIn->Write();
4292ca60 237 kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
238 kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
ad04a83c 239 kinematicEfficiencyIn->Write();
240 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
4292ca60 241 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
ad04a83c 242 fDetectorResponse->Write();
4292ca60 243 // optional histograms
244 if(fSaveFull) {
245 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
246 fSpectrumIn->Write();
247 fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
248 fDptInDist->Write();
249 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
250 fDptIn->Write();
251 fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
252 fFullResponseIn->Write();
253 }
ad04a83c 254 fActiveDir->cd();
255 TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
256 dirOut->cd();
51e6bc5a 257 switch (fUnfoldingAlgorithm) {
258 case kChi2 : {
259 convergedOut = UnfoldSpectrumChi2(
260 resizedJetPtOut,
261 resizedResonseOut,
262 kinematicEfficiencyOut,
263 unfoldingTemplateOut,
264 fUnfoldedOut,
265 TString("out"));
20abfcc4 266 printf(" > Spectrum (out of plane) unfolded using kChi2 < \n");
51e6bc5a 267 } break;
268 case kSVD : {
269 convergedOut = UnfoldSpectrumSVD(
270 resizedJetPtOut,
271 resizedResonseOut,
272 kinematicEfficiencyOut,
273 unfoldingTemplateOut,
4292ca60 274 fUnfoldedOut,
51e6bc5a 275 TString("out"));
20abfcc4 276 printf(" > Spectrum (out of plane) unfolded using kSVD < \n");
51e6bc5a 277 } break;
ef12d5a5 278 case kNone : { // do nothing, just rebin and optionally smooothen the spectrum
279 resizedResonseOut->SetNameTitle("measuredSpectrumOut", "measured spectrum, out plane");
280 if(fSmoothenSpectrum) { // see if we want to smooothen the spectrum
281 TH1D* resizedJetPtOutJagged((TH1D*)resizedJetPtOut->Clone("cachedRawJetJaggedOut"));
282 resizedJetPtOutJagged->SetNameTitle("measured spectrum before smoothening Out", "measured spectrum, before smoothening Out");
283 resizedJetPtOutJagged = ProtectHeap(resizedJetPtOutJagged);
284 resizedJetPtOutJagged->Write(); // save the original
285 resizedJetPtOut->Sumw2();
286 TFitResultPtr r = resizedJetPtOut->Fit(fPower, "QWILS", "", fFitMin, fFitMax);
287 if((int)r == 0) {
288 for(Int_t i(0); i < resizedJetPtOut->GetNbinsX() + 1; i++) {
289 if(resizedJetPtOut->GetBinCenter(i) > fFitStart) { // from this pt value use extrapolation
290 resizedJetPtOut->SetBinContent(i,fPower->Integral(resizedJetPtOut->GetXaxis()->GetBinLowEdge(i),resizedJetPtOut->GetXaxis()->GetBinUpEdge(i))/resizedJetPtOut->GetXaxis()->GetBinWidth(i));
291 }
292 }
293 }else printf(" > PANIC, SMOOTHENING FAILED < \n");
294
295 }
296 fUnfoldedOut = ProtectHeap(resizedJetPtOut, kTRUE, TString("out"));
297 convergedOut = kTRUE;
298 } break;
51e6bc5a 299 default : {
300 printf(" > Selected unfolding method is not implemented yet ! \n");
301 return;
302 }
303 }
4292ca60 304 resizedResonseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
305 resizedResonseOut->SetXTitle("p_{T}^{true} [GeV/c]");
306 resizedResonseOut->SetYTitle("p_{T}^{rec} [GeV/c]");
307 resizedResonseOut = ProtectHeap(resizedResonseOut);
ad04a83c 308 resizedResonseOut->Write();
4292ca60 309 kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
310 kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
ad04a83c 311 kinematicEfficiencyOut->Write();
312 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
4292ca60 313 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
ad04a83c 314 fDetectorResponse->Write();
4292ca60 315 // optional histograms
316 if(fSaveFull) {
ef12d5a5 317 fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
4292ca60 318 fSpectrumOut->Write();
319 fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
320 fDptOutDist->Write();
321 fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
322 fDptOut->Write();
323 fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
324 fFullResponseOut->Write();
325 }
51e6bc5a 326 // write general output histograms to file
ad04a83c 327 fActiveDir->cd();
4292ca60 328 if(convergedIn && convergedOut && fUnfoldedIn && fUnfoldedOut) {
ef12d5a5 329 TGraphErrors* ratio(GetRatio((TH1D*)fUnfoldedIn->Clone("unfoldedLocal_in"), (TH1D*)fUnfoldedOut->Clone("unfoldedLocal_out")));
4292ca60 330 if(ratio) {
331 ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
332 ratio->GetXaxis()->SetTitle("p_{T} [GeV/c]");
333 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
ef12d5a5 334 ratio = ProtectHeap(ratio);
4292ca60 335 ratio->Write();
336 // write histo values to RMS files if both routines converged
337 // input values are weighted by their uncertainty
338 for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
339 if(fUnfoldedIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), fUnfoldedIn->GetBinContent(i+1), 1./TMath::Power(fUnfoldedIn->GetBinError(i+1), 2.));
340 if(fUnfoldedOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), fUnfoldedOut->GetBinContent(i+1), 1./TMath::Power(fUnfoldedOut->GetBinError(i+1), 2.));
341 if(fUnfoldedOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), fUnfoldedIn->GetBinContent(i+1) / fUnfoldedOut->GetBinContent(i+1));
ef12d5a5 342 }
343 }
344 TGraphErrors* v2(GetV2((TH1D*)fUnfoldedIn->Clone("unfoldedLocal_inv2"), (TH1D*)fUnfoldedOut->Clone("unfoldedLocal_outv2")));
345 if(v2) {
346 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
347 v2->GetXaxis()->SetTitle("p_{T} [GeV/c]");
348 v2->GetYaxis()->SetTitle("v_{2}");
349 v2 = ProtectHeap(v2);
350 v2->Write();
351 }
352 } else if (fUnfoldedOut && fUnfoldedIn) {
353 TGraphErrors* ratio(GetRatio((TH1D*)fUnfoldedIn->Clone("unfoldedLocal_in"), (TH1D*)fUnfoldedOut->Clone("unfoldedLocal_out"), TString(""), fBinsRec->At(fBinsRec->GetSize()-1)));
354 if(ratio) {
355 ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
356 ratio->GetXaxis()->SetTitle("p_{T} [GeV/c]");
357 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
358 ratio = ProtectHeap(ratio);
359 ratio->Write();
360 }
361 TGraphErrors* v2(GetV2((TH1D*)fUnfoldedIn->Clone("unfoldedLocal_inv2"), (TH1D*)fUnfoldedOut->Clone("unfoldedLocal_outv2")));
362 if(v2) {
363 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
364 v2->GetXaxis()->SetTitle("p_{T} [GeV/c]");
365 v2->GetYaxis()->SetTitle("v_{2}");
366 v2 = ProtectHeap(v2);
367 v2->Write();
4292ca60 368 }
20abfcc4 369 }
ad04a83c 370 fDeltaPtDeltaPhi->Write();
371 fJetPtDeltaPhi->Write();
4292ca60 372 SaveConfiguration(convergedIn, convergedOut);
dc1455ee 373}
374//_____________________________________________________________________________
51e6bc5a 375Bool_t AliJetFlowTools::UnfoldSpectrumChi2(
4292ca60 376 TH1D* resizedJetPt, // truncated raw jets (same binning as pt rec of response)
377 TH2D* resizedResonse, // response matrix
378 TH1D* kinematicEfficiency, // kinematic efficiency
379 TH1D* unfoldingTemplate, // unfolding template: same binning is pt gen of response
380 TH1D *&unfolded, // will point to the unfolded spectrum
381 TString suffix) // suffix (in or out of plane)
dc1455ee 382{
51e6bc5a 383 // unfold the spectrum using chi2 minimization
384
4292ca60 385 // step 0) setup the static members of AliUnfolding
386 ResetAliUnfolding(); // reset from previous iteration
387 // also deletes and re-creates the global TVirtualFitter
388 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
389 if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
390 else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
391 if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
392 else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
393 AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
dc1455ee 394
4292ca60 395 // step 1) clone all input histograms.
396
397 // resizedJetPtLocal holds the spectrum that needs to be unfolded
51e6bc5a 398 TH1D *resizedJetPtLocal = (TH1D*)resizedJetPt->Clone(Form("resizedJetPtLocal_%s", suffix.Data()));
4292ca60 399 if(fSmoothenSpectrum) { // see if we want to smooothen the spectrum
400 TH1D* resizedJetPtLocalJagged((TH1D*)resizedJetPtLocal->Clone(Form("cachedRawJetLocalJagged_%s", suffix.Data())));
401 resizedJetPtLocalJagged->SetNameTitle(Form("measured spectrum before smoothening %s", suffix.Data()), Form("measured spectrum, before smoothening %s", suffix.Data()));
402 resizedJetPtLocalJagged = ProtectHeap(resizedJetPtLocalJagged);
403 resizedJetPtLocalJagged->Write(); // save the original
404 resizedJetPtLocal->Sumw2();
ef12d5a5 405 TFitResultPtr r = resizedJetPtLocal->Fit(fPower, "QWILS", "", fFitMin, fFitMax);
406 if((int)r == 0) {
407 for(Int_t i(0); i < resizedJetPtLocal->GetNbinsX() + 1; i++) {
408 if(resizedJetPtLocal->GetBinCenter(i) > fFitStart) { // from this pt value use extrapolation
409 resizedJetPtLocal->SetBinContent(i,fPower->Integral(resizedJetPtLocal->GetXaxis()->GetBinLowEdge(i),resizedJetPtLocal->GetXaxis()->GetBinUpEdge(i))/resizedJetPtLocal->GetXaxis()->GetBinWidth(i));
410 }
4292ca60 411 }
ef12d5a5 412 }else printf(" > PANIC, SMOOTHENING FAILED < \n");
413
4292ca60 414 }
415 // unfolded local will be filled with the result of the unfolding
416 TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
417
dc1455ee 418 // full response matrix and kinematic efficiency
51e6bc5a 419 TH2D* resizedResponseLocal = (TH2D*)resizedResonse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
420 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
4292ca60 421 // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
422 TH1D *priorLocal = (TH1D*)unfoldingTemplate->Clone(Form("priorLocal_%s", suffix.Data()));
ef12d5a5 423 if(fSmoothenSpectrum) {
424 TH1D* priorLocalJagged((TH1D*)priorLocal->Clone(Form("cachedRawJetLocalJagged_%s", suffix.Data())));
425 priorLocalJagged->SetNameTitle(Form("prior before smoothening %s", suffix.Data()), Form("prior before smoothening %s", suffix.Data()));
426 priorLocalJagged = ProtectHeap(priorLocalJagged);
427 priorLocalJagged->Write(); // save the original
428 priorLocal->Sumw2();
429 TFitResultPtr r = priorLocal->Fit(fPower, "WQILS", "", fFitMin, fFitMax);
430 if((int)r == 0) {
431 for(Int_t i(0); i < priorLocal->GetNbinsX() + 1; i++) {
432 if(priorLocal->GetBinCenter(i) > fFitStart) { // from this pt value use extrapolation
433 priorLocal->SetBinContent(i,fPower->Integral(priorLocal->GetXaxis()->GetBinLowEdge(i),priorLocal->GetXaxis()->GetBinUpEdge(i))/priorLocal->GetXaxis()->GetBinWidth(i));
434 }
435 }
436 }else printf(" > PANIC, SMOOTHENING FAILED < \n");
437
438 }
4292ca60 439 // step 2) start the unfolding
51e6bc5a 440 Int_t status(-1), i(0);
441 while(status < 0 && i < 100) {
4292ca60 442 // i > 0 means that the first iteration didn't converge. in that case, the result of the first
443 // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
444 if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
445 status = AliUnfolding::Unfold(
446 resizedResponseLocal, // response matrix
447 kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
448 resizedJetPtLocal, // measured spectrum
449 priorLocal, // initial conditions (set NULL to use measured spectrum)
450 unfoldedLocal); // results
451 // status holds the minuit fit status (where 0 means convergence)
51e6bc5a 452 i++;
453 }
4292ca60 454 // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
455 if(status == 0 && gMinuit->fISW[1] == 3) {
456 // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
ad04a83c 457 TVirtualFitter *fitter(TVirtualFitter::GetFitter());
4292ca60 458 if(gMinuit) gMinuit->Command("SET COV");
459 TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
460 TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
461 pearson->Print();
ad04a83c 462 TH2D *hPearson(new TH2D(*pearson));
4292ca60 463 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
464 hPearson = ProtectHeap(hPearson);
ad04a83c 465 hPearson->Write();
4292ca60 466 } else status = -1;
467 // step 3) refold the unfolded spectrum and save the ratio measured / refolded
468 TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
469 foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
470 unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
471 TGraphErrors* ratio(GetRatio(foldedLocal, resizedJetPtLocal, kTRUE, fBinsTrue->At(0), fBinsTrue->At(fBinsTrue->GetSize()-1)));
472 if(ratio) {
473 ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio refolded and measured spectrum %s plane", suffix.Data()));
474 ratio = ProtectHeap(ratio);
475 ratio->Write();
476 }
477 resizedJetPtLocal = ProtectHeap(resizedJetPtLocal);
478 resizedJetPtLocal->Write();
479 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
480 resizedResponseLocal->Write();
481 unfoldedLocal = ProtectHeap(unfoldedLocal);
51e6bc5a 482 unfoldedLocal->Write();
483 unfolded = unfoldedLocal;
4292ca60 484 foldedLocal = ProtectHeap(foldedLocal);
485 foldedLocal->Write();
ef12d5a5 486 priorLocal = ProtectHeap(priorLocal);
487 priorLocal->Write();
4292ca60 488 return (status == 0) ? kTRUE : kFALSE;
dc1455ee 489}
490//_____________________________________________________________________________
51e6bc5a 491Bool_t AliJetFlowTools::UnfoldSpectrumSVD(
4292ca60 492 TH1D* resizedJetPt, // jet pt in pt rec bins
493 TH2D* resizedResonse, // full response matrix, normalized in slides of pt true
494 TH1D* kinematicEfficiency, // kinematic efficiency
ef12d5a5 495 TH1D* unfoldingTemplate, // jet pt in pt true bins, also the prior when measured is chosen as prior
496 TH1D *&unfolded, // will point to result. temporarily holds prior when chi2 is chosen as prior
4292ca60 497 TString suffix) // suffix (in, out)
51e6bc5a 498{
499 // use SVD (singular value decomposition) method to unfold spectra
500
4292ca60 501 // 1) get a prior for unfolding.
502 // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
51e6bc5a 503 TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
504 dirOut->cd();
20abfcc4 505 switch (fPrior) { // select the prior for unfolding
506 case kPriorChi2 : {
ef12d5a5 507 if(fBinsTruePrior && fBinsRecPrior) { // if set, use different binning for the prior
508 TArrayD* tempArrayTrue(fBinsTrue); // temporarily cache the original (SVD) binning
509 fBinsTrue = fBinsTruePrior; // switch binning schemes (will be used in UnfoldSpectrumChi2())
510 TArrayD* tempArrayRec(fBinsRec);
511 fBinsRec = fBinsRecPrior;
512 TH1D* resizedJetPtChi2 = GetUnfoldingTemplate((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"));
513 TH1D* unfoldingTemplateChi2 = GetUnfoldingTemplate((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"));
514 TH2D* resizedResonseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
515 TH1D* kinematicEfficiencyChi2(resizedResonseChi2->ProjectionX());
516 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
517 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
518 if(! UnfoldSpectrumChi2(
519 resizedJetPtChi2,
520 resizedResonseChi2,
521 kinematicEfficiencyChi2,
522 unfoldingTemplateChi2, // prior for chi2 unfolding (measured)
523 unfolded, // will hold the result from chi2 (is prior for SVD)
524 TString(Form("prior_%s", suffix.Data()))) ) {
525 printf(" > UnfoldSVD:: panic, couldn't get prior from Chi2 unfolding! \n");
526 printf(" probably Chi2 unfolding did not converge < \n");
527 return kFALSE;
528 }
529 fBinsTrue = tempArrayTrue; // reset bins borders
530 fBinsRec = tempArrayRec;
531 unfolded = GetUnfoldingTemplate(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data()))); // rebin unfolded
532 } else if(! UnfoldSpectrumChi2(
20abfcc4 533 resizedJetPt,
534 resizedResonse,
535 kinematicEfficiency,
ef12d5a5 536 unfoldingTemplate, // prior for chi2 unfolding (measured)
537 unfolded, // will hold the result from chi2 (is prior for SVD)
4292ca60 538 TString(Form("prior_%s", suffix.Data()))) ) {
539 printf(" > UnfoldSVD:: panic, couldn't get prior from Chi2 unfolding! \n");
540 printf(" probably Chi2 unfolding did not converge < \n");
541 return kFALSE;
542 }
543 if(!unfolded) {
544 printf(" > UnfoldSVD:: panic, Chi2 unfolding converged but the prior is NULL ! < " );
545 return kFALSE;
546 }
20abfcc4 547 break;
548 }
4292ca60 549 case kPriorMeasured : {
ef12d5a5 550 unfolded = (TH1D*)unfoldingTemplate->Clone(Form("kPriorMeasured_%s", suffix.Data())); // copy template to unfolded to use as prior
551 if(fSmoothenSpectrum) { // optionally smoothen the measured prior
552 unfolded->Sumw2();
553 TFitResultPtr r = unfolded->Fit(fPower, "QWILS", "", fFitMin, fFitMax);
554 if((int)r == 0) {
555 for(Int_t i(1); i < unfolded->GetNbinsX() + 1; i++) {
556 if(unfolded->GetBinCenter(i) > fFitStart) { // from this pt value use extrapolation
557 unfolded->SetBinContent(i,fPower->Integral(unfolded->GetXaxis()->GetBinLowEdge(i),unfolded->GetXaxis()->GetBinUpEdge(i))/unfolded->GetXaxis()->GetBinWidth(i));
558 }
559 }
560 }else printf(" > PANIC, SMOOTHENING FAILED < \n");
561 }
562// unfoldingTemplate = ProtectHeap(unfolded, kFALSE, TString("template"));
4292ca60 563 }
20abfcc4 564 default : break;
565 }
566 // note: true and measured spectrum must have same binning for SVD unfolding
ef12d5a5 567 // a sane starting point for regularization is nbins / 2 (but the user has to set this ! )
20abfcc4 568 if(unfoldingTemplate->GetXaxis()->GetNbins() != resizedJetPt->GetXaxis()->GetNbins()) {
4292ca60 569 printf(" > UnfoldSpectrumSVD:: PANIC, true and measured spectrum must have same numer of bins ! < \n ");
20abfcc4 570 }
571 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
4292ca60 572 cout << " 1) retrieved prior " << endl;
20abfcc4 573
4292ca60 574 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
575 // prior
20abfcc4 576 TH1D *unfoldedLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
4292ca60 577 // raw jets in pt rec binning
51e6bc5a 578 TH1D *cachedRawJetLocal((TH1D*)resizedJetPt->Clone(Form("jets_%s", suffix.Data())));
ef12d5a5 579 // raw jets in pt true binning
51e6bc5a 580 TH1D *cachedRawJetLocalCoarse((TH1D*)unfoldingTemplate->Clone(Form("unfoldingTemplate_%s", suffix.Data())));
ef12d5a5 581 // copy of raw jets in pt true binning
51e6bc5a 582 TH1D *cachedRawJetLocalCoarseOrig((TH1D*)cachedRawJetLocalCoarse->Clone(Form("cachedRawJetLocalCoarseOrig_%s", suffix.Data())));
ef12d5a5 583 // local copies response matrix
51e6bc5a 584 TH2D *cachedResponseLocal((TH2D*)resizedResonse->Clone(Form("cachedResponseLocal_%s", suffix.Data())));
ef12d5a5 585 // local copy of response matrix, all true slides normalized to 1
51e6bc5a 586 TH2D *cachedResponseLocalNorm((TH2D*)resizedResonse->Clone(Form("cachedResponseLocalNorm_%s", suffix.Data())));
ef12d5a5 587 cachedResponseLocalNorm = NormalizeTH2D(cachedResponseLocalNorm);
588 // kinematic efficiency
51e6bc5a 589 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
ef12d5a5 590 // place holder histos
4292ca60 591 TH1 *unfoldedLocalSVD(0x0);
592 TH1 *foldedLocalSVD(0x0);
593 cout << " 2) setup necessary input " << endl;
51e6bc5a 594 // 3) configure routine
595 RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
20abfcc4 596 // prior: use fit for where the histogram is sparsely filled
ef12d5a5 597 if(fSmoothenSpectrum) { // smoothen raw input jets in true and rec binning, and smoothen the prior
4292ca60 598 cachedRawJetLocalCoarse->Sumw2();
ef12d5a5 599 TFitResultPtr r = cachedRawJetLocalCoarse->Fit(fPower, "QWILS", "", fFitMin, fFitMax);
600 if((int)r == 0) {
601 for(Int_t i(1); i < cachedRawJetLocalCoarse->GetNbinsX() + 1; i++) {
602 if(cachedRawJetLocalCoarse->GetBinCenter(i) > fFitStart) {
603 cachedRawJetLocalCoarse->SetBinContent(i,fPower->Integral(cachedRawJetLocalCoarse->GetXaxis()->GetBinLowEdge(i),cachedRawJetLocalCoarse->GetXaxis()->GetBinUpEdge(i))/cachedRawJetLocalCoarse->GetXaxis()->GetBinWidth(i));
604 }
4292ca60 605 }
ef12d5a5 606 } else printf(" > PANIC, SMOOTHENING FAILED < \n");
4292ca60 607 TH1D* cachedRawJetLocalJagged((TH1D*)cachedRawJetLocal->Clone("cachedRawJetLocalJagged"));
608 cachedRawJetLocalJagged->SetNameTitle("measured spectrum before smoothening", "measured spectrum, before smoothening");
609 cachedRawJetLocalJagged->Write();
610 cachedRawJetLocal->Sumw2();
ef12d5a5 611 r = cachedRawJetLocal->Fit(fPower, "S", "", fFitMin, fFitMax);
612 if((int)r == 0) {
613 for(Int_t i(0); i < cachedRawJetLocal->GetNbinsX() + 1; i++) {
614 if(cachedRawJetLocal->GetBinCenter(i) > fFitStart) { // from this pt value use extrapolation
615 cachedRawJetLocal->SetBinContent(i,fPower->Integral(cachedRawJetLocal->GetXaxis()->GetBinLowEdge(i),cachedRawJetLocal->GetXaxis()->GetBinUpEdge(i))/cachedRawJetLocal->GetXaxis()->GetBinWidth(i));
616 }
617 }
618 }
619 r = unfoldedLocal->Fit(fPower, "S", "", fFitMin, fFitMax);
620 if((int)r == 0) {
621 for(Int_t i(0); i < unfoldedLocal->GetNbinsX() + 1; i++) {
622 if(unfoldedLocal->GetBinCenter(i) > fFitStart) { // from this pt value use extrapolation
623 unfoldedLocal->SetBinContent(i,fPower->Integral(unfoldedLocal->GetXaxis()->GetBinLowEdge(i),unfoldedLocal->GetXaxis()->GetBinUpEdge(i))/unfoldedLocal->GetXaxis()->GetBinWidth(i));
624 }
4292ca60 625 }
51e6bc5a 626 }
4292ca60 627 }
628 cout << " 3) configured routine " << endl;
629 // 4) get transpose matrices
630 // a) get the transpose matrix for the prior
631 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(cachedResponseLocal));
632 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
633 // normalize it with the prior
ef12d5a5 634 responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, unfoldedLocal);
4292ca60 635 cout << " 4a) retrieved first transpose matrix " << endl;
636 // b) prior norm
637 TH2* responseMatrixLocalTransposePriorNorm(fResponseMaker->GetTransposeResponsMatrix(cachedResponseLocalNorm));
638 responseMatrixLocalTransposePriorNorm->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePriorNorm->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePriorNorm->GetName(), suffix.Data()));
639 // normalize with the prior
640 responseMatrixLocalTransposePriorNorm = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePriorNorm, unfoldedLocal);
641 cout << " 4b) retrieved second transpose matrix " << endl;
642
643 // 5) get response for SVD unfolding
644 RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
51e6bc5a 645
4292ca60 646 // change to inplane dir
20abfcc4 647 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) :fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
4292ca60 648
649 cout << " 5) retrieved roo unfold response object " << endl;
650 // 6) actualy unfolding loop
651 RooUnfoldSvd unfoldSVD(&responseSVD, cachedRawJetLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
20abfcc4 652 unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
653 TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
654 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
ef12d5a5 655 cout << " Pearson coeffs" << endl;
4292ca60 656 // create the unfolding qa plots
657 cout << " 6) unfolded spectrum " << endl;
658 if(pearson) {
659 TH2D* hPearson = new TH2D(*pearson);
660 pearson->Print();
661 hPearson->SetNameTitle("PearsonCoefficients", "Pearson coefficients");
662 hPearson->Write();
663 } else return kFALSE; // return if unfolding didn't converge
664 // correct for the efficiency
ef12d5a5 665 unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
4292ca60 666
667 // plot singular values and d_i vector
20abfcc4 668 TSVDUnfold* svdUnfold(unfoldSVD.Impl());
51e6bc5a 669 TH1* hSVal(svdUnfold->GetSV());
670 TH1D* hdi(svdUnfold->GetD());
4292ca60 671 hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
51e6bc5a 672 hSVal->SetXTitle("singular values");
673 hSVal->Write();
4292ca60 674 hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
51e6bc5a 675 hdi->SetXTitle("|d_{i}^{kreg}|");
676 hdi->Write();
4292ca60 677 cout << " plotted singular values and d_i vector " << endl;
678
679 // 7) refold the unfolded spectrum
680 foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, cachedResponseLocalNorm,kinematicEfficiencyLocal);
681 TGraphErrors* ratio(GetRatio(cachedRawJetLocal, foldedLocalSVD, "ratio measured / re-folded"));
682 ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
683 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
684 ratio->Write();
685 cout << " 7) refolded the unfolded spectrum " << endl;
686
687 // write to output
688 cachedRawJetLocal->SetNameTitle("input spectrum","input spectrum (measured)");
689 cachedRawJetLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
690 cachedRawJetLocal->Write(); // input spectrum
691 unfoldedLocalSVD->SetNameTitle("unfolded spectrum","unfolded spectrum");
692 unfoldedLocalSVD->Write(); // unfolded spectrum
693 foldedLocalSVD->SetNameTitle(Form("refoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
694 foldedLocalSVD->Write(); // re-folded spectrum
695 // switch back to active root directory
20abfcc4 696 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) :fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
4292ca60 697 responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix");
698 responseMatrixLocalTransposePrior->SetXTitle("p_{T}^{true} [GeV/c]");
699 responseMatrixLocalTransposePrior->SetYTitle("p_{T}^{rec} [GeV/c]");
700 responseMatrixLocalTransposePrior->Write();
701 cachedRawJetLocal->SetNameTitle("PriorOriginal", "Prior, original");
702 cachedRawJetLocal->SetXTitle("p_{t} [GeV/c]");
703 cachedRawJetLocalCoarse->SetNameTitle("PriorSmoothened", "Prior, smoothened");
704 cachedRawJetLocalCoarse->SetXTitle("p_{t} [GeV/c]");
705 cachedRawJetLocalCoarse->Write();
706 cachedRawJetLocalCoarseOrig->SetNameTitle("Prior", "Prior");
707 cachedRawJetLocalCoarseOrig->SetXTitle("p_{t} [GeV/c]");
708 cachedRawJetLocalCoarseOrig->Write();
709 unfolded = (TH1D*)unfoldedLocalSVD;
ef12d5a5 710 cachedResponseLocalNorm = ProtectHeap(cachedResponseLocalNorm);
711 cachedResponseLocalNorm->Write();
20abfcc4 712 return (unfoldedLocalSVD) ? kTRUE : kFALSE;
51e6bc5a 713}
714//_____________________________________________________________________________
dc1455ee 715Bool_t AliJetFlowTools::PrepareForUnfolding()
716{
717 // prepare for unfolding
4292ca60 718 if(fRawInputProvided) return kTRUE;
dc1455ee 719 if(!fInputList) {
720 printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
721 return kFALSE;
722 }
723 if(!fDetectorResponse) {
724 printf(" AliJetFlowTools::PrepareForUnfolding() fDetectorResponse not found \n - Set detector response using AliJetFlowTools::SetDetectorResponse() \n ");
725 return kFALSE;
726 }
4292ca60 727 // check if the pt bin for true and rec have been set
728 if(!fBinsTrue || !fBinsRec) {
729 printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
730 return kFALSE;
dc1455ee 731 }
4292ca60 732 if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
733 // procedures, these profiles will be nonsensical, user is responsible
734 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
735 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
736 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
dc1455ee 737 }
ef12d5a5 738 if(!fTrainPower) {
739 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
740 }
dc1455ee 741 // extract the spectra
742 TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityBin));
ad04a83c 743 fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
744 if(!fJetPtDeltaPhi) {
dc1455ee 745 printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
746 return kFALSE;
747 }
4292ca60 748 fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
dc1455ee 749 // in plane spectrum
ef12d5a5 750 if(fNoDphi) {
751 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40);
752 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40);
753 } else {
754 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10);
755 fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40));
756 fSpectrumIn = ProtectHeap(fSpectrumIn);
757 // out of plane spectrum
758 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30);
759 fSpectrumOut = ProtectHeap(fSpectrumOut);
760 }
20abfcc4 761 // normalize spectra to event count if requested
762 if(fNormalizeSpectra) {
763 TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityBin)));
ef12d5a5 764 if(!rho) return 0x0;
765 Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
766 if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
20abfcc4 767 if(fEventCount > 0) {
4292ca60 768 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
769 fSpectrumOut->Sumw2();
ef12d5a5 770 Double_t pt(0);
771 // Double_t error(0); // lots of issues with the errors here ...
772 for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
4292ca60 773 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount; // normalized count
ef12d5a5 774// error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
4292ca60 775 fSpectrumIn->SetBinContent(1+i, pt);
ef12d5a5 776 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
777// if(error > 0) fSpectrumIn->SetBinError(1+i, TMath::Sqrt(error));
778 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
4292ca60 779 }
ef12d5a5 780 for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
781 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount; // normalized count
782// error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
783 fSpectrumOut->SetBinContent(1+i, pt);
784 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
785// if(error > 0) fSpectrumOut->SetBinError(1+i, TMath::Sqrt(error));
786 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
4292ca60 787 }
20abfcc4 788 }
ef12d5a5 789 if(normalizeToFullSpectrum) fEventCount = -1;
20abfcc4 790 }
dc1455ee 791 // extract the delta pt matrices
792 TString deltaptName(Form("fHistDeltaPtDeltaPhi2_%i", fCentralityBin));
ad04a83c 793 fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
794 if(!fDeltaPtDeltaPhi) {
dc1455ee 795 printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
796 }
4292ca60 797 fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
798 // in plane delta pt distribution
ef12d5a5 799 if(fNoDphi) {
800 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40);
801 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40);
802 } else {
803 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10);
804 fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40));
805 // out of plane delta pt distribution
806 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30);
807 fDptInDist = ProtectHeap(fDptInDist);
808 fDptOutDist = ProtectHeap(fDptOutDist);
809 // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
810 }
4292ca60 811
dc1455ee 812 // create a rec - true smeared response matrix
813 TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
814 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
ad04a83c 815 Bool_t skip = kFALSE;
dc1455ee 816 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
ad04a83c 817 (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
818 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
dc1455ee 819 }
820 }
821 TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
822 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
ad04a83c 823 Bool_t skip = kFALSE;
dc1455ee 824 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
ad04a83c 825 (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
826 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
dc1455ee 827 }
828 }
829 fDptIn = new TH2D(*rfIn);
830 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
831 fDptIn->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
832 fDptIn->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
4292ca60 833 fDptIn = ProtectHeap(fDptIn);
dc1455ee 834 fDptOut = new TH2D(*rfOut);
835 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
836 fDptOut->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
837 fDptOut->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
4292ca60 838 fDptOut = ProtectHeap(fDptOut);
839
840 fRefreshInput = kTRUE; // force cloning of the input
dc1455ee 841 return kTRUE;
842}
843//_____________________________________________________________________________
844TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
845 // resize the x-axis of a th1d
846 if(!histo) {
847 printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
848 return NULL;
849 }
850 // see how many bins we need to copy
851 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);
852 // low is the bin number of the first new bin
853 Int_t l = histo->GetXaxis()->FindBin(low);
854 // set the values
855 Double_t _x(0), _xx(0);
856 for(Int_t i(0); i < up-low; i++) {
857 _x = histo->GetBinContent(l+i);
858 _xx=histo->GetBinError(l+i);
859 resized->SetBinContent(i+1, _x);
860 resized->SetBinError(i+1, _xx);
861 }
862 return resized;
863}
864//_____________________________________________________________________________
865TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
866 // resize the y-axis of a th2d
867 if(!histo) {
868 printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
869 return NULL;
870 }
871 // see how many bins we need to copy
872 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());
873 // assume only the y-axis has changed
874 // low is the bin number of the first new bin
875 Int_t low = histo->GetYaxis()->FindBin(y->At(0));
876 // set the values
877 Double_t _x(0), _xx(0);
878 for(Int_t i(0); i < x->GetSize(); i++) {
879 for(Int_t j(0); j < y->GetSize(); j++) {
880 _x = histo->GetBinContent(i, low+j);
881 _xx=histo->GetBinError(i, low+1+j);
882 resized->SetBinContent(i, j, _x);
883 resized->SetBinError(i, j, _xx);
884 }
885 }
886 return resized;
887}
888//_____________________________________________________________________________
889TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo) {
890 // general method to normalize all vertical slices of a th2 to unity
891 // i.e. get a probability matrix
892 if(!histo) {
893 printf(" > NormalizeTH2D:: fatal error, NULL pointer passed < \n");
894 return NULL;
895 }
896 Int_t binsX = histo->GetXaxis()->GetNbins();
897 Int_t binsY = histo->GetYaxis()->GetNbins();
898
899 // normalize all slices in x
900 for(Int_t i(0); i < binsX; i++) { // for each vertical slice
901 Double_t weight = 0;
902 for(Int_t j(0); j < binsY; j++) { // loop over all the horizontal components
903 weight+=histo->GetBinContent(i+1, j+1);
904 } // now we know the total weight
905 for(Int_t j(0); j < binsY; j++) {
906 if (weight <= 0 ) continue;
907 histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
908 histo->SetBinError( 1+i, j+1, histo->GetBinError( 1+i, j+1)/weight);
909 }
910 }
911 return histo;
912}
913//_____________________________________________________________________________
914TH1D* AliJetFlowTools::GetUnfoldingTemplate(TH1D* histo, TArrayD* bins, TString suffix) {
915 // return a TH1D with the supplied histogram rebinned to the supplied bins
916 // this histogram will be used as a startng point for the chi^2 minimization
4292ca60 917 // the returned histogram is new
dc1455ee 918 if(!histo || !bins) {
919 printf(" > RebinTH2D:: fatal error, NULL pointer passed < \n");
920 return NULL;
921 }
922 // create the output histo
923 TString name = histo->GetName();
924 name+="_template";
925 name+=suffix;
926 TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
927 for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
4292ca60 928 // loop over the bins of the old histo and fill the new one with its data
929 rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
dc1455ee 930 }
931 return rebinned;
932}
933//_____________________________________________________________________________
4292ca60 934TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
935 if(!fResponseMaker || !binsTrue || !binsRec) {
936 printf(" > RebinTH2D:: function called with NULL arguments < \n");
937 return 0x0;
dc1455ee 938 }
4292ca60 939 TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
940 return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
dc1455ee 941}
942//_____________________________________________________________________________
943TH2D* AliJetFlowTools::MatrixMultiplicationTH2D(TH2D* A, TH2D* B, TString name) {
944 // general matrix multiplication
945 if(!A || !B) {
946 printf(" > MatrixMultiplicationTH2D:: fatal error, NULL pointer passed < \n");
947 return NULL;
948 }
949 // step A -> see if the matrices are both square and of equal dimensions
950 Int_t aX(A->GetXaxis()->GetNbins());
951 Int_t aY(A->GetYaxis()->GetNbins());
952 Int_t bX(B->GetXaxis()->GetNbins());
953 Int_t bY(B->GetYaxis()->GetNbins());
954 if(!(aX == aY && aX == bX && aX == bY)) {
955 printf(" > MatrixMultiplicationTH2D:: Error, matrices with incompatible dimensions passed ! < \n");
956 printf(" - matrix A (%i, %i) \n", aX, aY);
957 printf(" - matrix B (%i, %i) \n", bX, bY);
958 return NULL;
959 }
960
961 // step B -> setup a new matrix template to store the results
962 TH2D* C = (TH2D*)A->Clone();
963 C->SetNameTitle(name.Data(), name.Data());
964 C->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
965 C->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
966
967 // step C -> nested loop multiplication
968 printf ( " >> parsing matrix data, may take some time (depending on matrix dimensions) << \n");
969 for(Int_t i(0); i < aX; i++) { // x coordinate of target matrix
970 for(Int_t j(0); j < aY; j++) { // y coordinate of target matrix
971
51e6bc5a 972 // the coordinate of the target matrix (x, y) is the dot
973 // product of row x of matrix A with column y of matrix B
dc1455ee 974 // so in this nexted loop, we need to fill point i, j of the target matrix
975 // with the dot product of
976 // matrix A) row i \cdotp matrix B) column j
977 // looping through a row means that y remains constant
978 // looping through a column means that x remains constant
979
980 Double_t value(0); // placeholder for target value
981 for(Int_t x(0); x < aX; x++) value += A->GetBinContent(x, i) * B->GetBinContent(j, x);
982 C->SetBinContent(i, j);
983 }
984 }
985 return C;
986}
987//_____________________________________________________________________________
988TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale) {
989 // normalize a th1d to a certain scale
4292ca60 990 histo->Sumw2();
991 Double_t integral = histo->Integral()*scale;
992 if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
993 else if (scale != 1.) histo->Scale(1./scale, "width");
994 else printf(" > Histogram integral < 0, cannot normalize \n");
995 return histo;
dc1455ee 996}
997//_____________________________________________________________________________
51e6bc5a 998TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
dc1455ee 999{
1000 // Calculate pearson coefficients from covariance matrix
51e6bc5a 1001 TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1002 Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
ad04a83c 1003 Double_t pearson(0.);
dc1455ee 1004 if(nrows==0 && ncols==0) return 0x0;
ef12d5a5 1005 for(Int_t row = 0; row < nrows; row++) {
1006 for(Int_t col = 0; col<ncols; col++) {
51e6bc5a 1007 if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1008 (*pearsonCoefficients)(row,col) = pearson;
dc1455ee 1009 }
1010 }
51e6bc5a 1011 return pearsonCoefficients;
dc1455ee 1012}
1013//_____________________________________________________________________________
4292ca60 1014Bool_t AliJetFlowTools::SetRawInput (
1015 TH2D* detectorResponse, // detector response matrix
1016 TH1D* jetPtIn, // in plane jet spectrum
1017 TH1D* jetPtOut, // out of plane jet spectrum
1018 TH1D* dptIn, // in plane delta pt distribution
1019 TH1D* dptOut, // out of plane delta pt distribution
1020 Int_t eventCount) {
1021 // set input histograms manually
1022 fDetectorResponse = detectorResponse;
1023 fSpectrumIn = jetPtIn;
1024 fSpectrumOut = jetPtOut;
1025 fDptInDist = dptIn;
1026 fDptOutDist = dptOut;
1027 fRawInputProvided = kTRUE;
1028 // check if all data is provided
1029 if(!fDetectorResponse) {
1030 printf(" fDetectorResponse not found \n ");
1031 return kFALSE;
1032 }
1033 // check if the pt bin for true and rec have been set
1034 if(!fBinsTrue || !fBinsRec) {
1035 printf(" No true or rec bins set, please set binning ! \n");
1036 return kFALSE;
1037 }
1038 if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
1039 // procedures, these profiles will be nonsensical, user is responsible
1040 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1041 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1042 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1043 }
1044 // normalize spectra to event count if requested
1045 if(fNormalizeSpectra) {
1046 fEventCount = eventCount;
1047 if(fEventCount > 0) {
1048 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
1049 fSpectrumOut->Sumw2();
1050 fSpectrumIn->Scale(1./((double)fEventCount));
1051 fSpectrumOut->Scale(1./((double)fEventCount));
1052 }
1053 }
1054 if(!fNormalizeSpectra && fEventCount > 0) {
1055 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
1056 fSpectrumOut->Sumw2();
1057 fSpectrumIn->Scale(1./((double)fEventCount));
1058 fSpectrumOut->Scale(1./((double)fEventCount));
1059 }
1060 fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
1061 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityBin), Form("dpt_response_INPLANE_%i", fCentralityBin));
1062 fDptIn->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1063 fDptIn->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1064 fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
1065 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityBin), Form("dpt_response_OUTOFPLANE_%i", fCentralityBin));
1066 fDptOut->GetXaxis()->SetTitle("p_{T}^{gen} [GeV/c]");
1067 fDptOut->GetYaxis()->SetTitle("p_{T}^{rec} [GeV/c]");
1068
1069 return kTRUE;
1070}
1071//_____________________________________________________________________________
1072TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
dc1455ee 1073{
ef12d5a5 1074 // return ratio of h1 / h2
1075 // histograms can have different binning. errors are propagated as uncorrelated
20abfcc4 1076 if(!(h1 && h2) ) {
1077 printf(" GetRatio called with NULL argument(s) \n ");
1078 return 0x0;
1079 }
ad04a83c 1080 Int_t j(0);
4292ca60 1081 TGraphErrors *gr = new TGraphErrors();
ef12d5a5 1082 Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
4292ca60 1083 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
dc1455ee 1084 binCent = h1->GetXaxis()->GetBinCenter(i);
4292ca60 1085 if(xmax > 0. && binCent > xmax) continue;
dc1455ee 1086 j = h2->FindBin(binCent);
1087 binWidth = h1->GetXaxis()->GetBinWidth(i);
ad04a83c 1088 if(h2->GetBinContent(j) > 0.) {
dc1455ee 1089 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
1090 Double_t A = 1./h2->GetBinContent(j)*h1->GetBinError(i);
4292ca60 1091 Double_t B = 0.;
1092 if(h2->GetBinError(j)>0.) {
1093 B = -1.*h1->GetBinContent(i)/(h2->GetBinContent(j)*h2->GetBinContent(j))*h2->GetBinError(j);
1094 error2 = A*A + B*B;
ef12d5a5 1095 } else error2 = A*A;
1096 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
dc1455ee 1097 gr->SetPoint(gr->GetN(),binCent,ratio);
ef12d5a5 1098 gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
dc1455ee 1099 }
1100 }
ad04a83c 1101 if(appendFit) {
4292ca60 1102 TF1* fit(new TF1("lin", "pol0", 10, 100));
ad04a83c 1103 gr->Fit(fit);
1104 }
4292ca60 1105 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
dc1455ee 1106 return gr;
1107}
1108//_____________________________________________________________________________
ef12d5a5 1109TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
1110{
1111 // get v2 from difference of in plane, out of plane yield
1112 // h1 must hold the in-plane yield, h2 holds the out of plane yield
1113 // different binning is allowed
1114 // r is the event plane resolution for the chosen centrality
1115 if(!(h1 && h2) ) {
1116 printf(" GetV2 called with NULL argument(s) \n ");
1117 return 0x0;
1118 }
1119 Int_t j(0);
1120 TGraphErrors *gr = new TGraphErrors();
1121 Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
1122 Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
1123 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
1124 binCent = h1->GetXaxis()->GetBinCenter(i);
1125 j = h2->FindBin(binCent);
1126 binWidth = h1->GetXaxis()->GetBinWidth(i);
1127 if(h2->GetBinContent(j) > 0.) {
1128 in = h1->GetBinContent(i);
1129 ein = h1->GetBinError(i);
1130 out = h2->GetBinContent(j);
1131 eout = h2->GetBinError(j);
1132 ratio = pre*((in-out)/(in+out));
1133 error2 = (r*4.)/(TMath::Pi())*((out*out/(TMath::Power(in+out, 4)))*ein*ein+(in*in/(TMath::Power(in+out, 4)))*eout*eout);
1134 if(error2 > 0) error2 = TMath::Sqrt(error2);
1135 gr->SetPoint(gr->GetN(),binCent,ratio);
1136 gr->SetPointError(gr->GetN()-1,0.5*binWidth,error2);
1137 }
1138 }
1139 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
1140 return gr;
1141}
1142//_____________________________________________________________________________
4292ca60 1143void AliJetFlowTools::WriteObject(TObject* object) {
1144 // write object with unique identifier to active TDirectoryFile
1145 if(!object) {
1146 printf(" > WriteObject:: called with NULL arguments \n ");
1147 return;
1148 } else object->Write();
1149 // FIXME to be implememnted
1150}
1151//_____________________________________________________________________________
1152TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
1153 // construt a delta pt response matrix from supplied dpt distribution
1154 // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
1155 // do a weighted rebinning to a (coarser) dpt distribution
1156 // be careful with the binning of the dpt response: it should be equal to that
1157 // of the response matrix, otherwise dpt and response matrices cannot be multiplied
1158 //
1159 // the response matrix will be square and have the same binning
1160 // (min, max and granularity) of the input histogram
1161 Int_t bins(dpt->GetXaxis()->GetNbins()); // number of bins, will also be no of rows, columns
1162 Double_t _bins[bins+1]; // prepare array with bin borders
1163 for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
1164 _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1); // get upper edge
1165 TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
1166 for(Int_t j(0); j < bins+1; j++) { // loop on pt true slices j
1167 Bool_t skip = kFALSE;
1168 for(Int_t k(0); k < bins+1; k++) { // loop on pt gen slices k
1169 (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
1170 if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
1171 }
1172 }
1173 return res;
1174}
ef12d5a5 1175//_____________________________________________________________________________
1176TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1177 if(!binsTrue || !binsRec) {
1178 printf(" > GetUnityResponse:: function called with NULL arguments < \n");
1179 return 0x0;
1180 }
1181 TString name(Form("unityResponse_%s", suffix.Data()));
1182 TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
1183 for(Int_t i(0); i < binsTrue->GetSize(); i++) {
1184 for(Int_t j(0); j < binsRec->GetSize(); j++) {
1185 if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
1186 }
1187 }
1188 return unity;
1189}
4292ca60 1190//_____________________________________________________________________________
1191void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) {
1192 // save configuration parameters to histogram
ef12d5a5 1193 TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 16, -.5, 16.5);
1194 summary->SetBinContent(1, fBetaIn);
1195 summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
1196 summary->SetBinContent(2, fBetaOut);
1197 summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
1198 summary->SetBinContent(3, fCentralityBin);
1199 summary->GetXaxis()->SetBinLabel(3, "fCentralityBin");
1200 summary->SetBinContent(4, (int)convergedIn);
1201 summary->GetXaxis()->SetBinLabel(4, "convergedIn");
1202 summary->SetBinContent(5, (int)convergedOut);
1203 summary->GetXaxis()->SetBinLabel(5, "convergedOut");
1204 summary->SetBinContent(6, (int)fAvoidRoundingError);
1205 summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
1206 summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
1207 summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
1208 summary->SetBinContent(8, (int)fPrior);
1209 summary->GetXaxis()->SetBinLabel(8, "fPrior");
1210 summary->SetBinContent(9, fSVDRegIn);
1211 summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
1212 summary->SetBinContent(10, fSVDRegOut);
1213 summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
1214 summary->SetBinContent(11, (int)fSVDToy);
1215 summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
1216 summary->SetBinContent(12, fJetRadius);
1217 summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
1218 summary->SetBinContent(13, (int)fNormalizeSpectra);
1219 summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
1220 summary->SetBinContent(14, (int)fSmoothenSpectrum);
1221 summary->GetXaxis()->SetBinLabel(14, "fSmoothenSpectrum");
1222 summary->SetBinContent(15, (int)fTestMode);
1223 summary->GetXaxis()->SetBinLabel(15, "fTestMode");
1224 summary->SetBinContent(16, (int)fUseDetectorResponse);
1225 summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
1226 summary->Write();
4292ca60 1227}
1228//_____________________________________________________________________________
1229void AliJetFlowTools::ResetAliUnfolding() {
1230 // ugly function: reset all unfolding parameters
1231 TVirtualFitter* fitter(TVirtualFitter::GetFitter());
1232 if(fitter) {
1233 printf(" > Found fitter, will delete it < \n");
1234 delete fitter;
1235 }
ef12d5a5 1236// if(gMinuit) {
1237// printf(" > Found gMinuit, will re-create it < \n");
1238// delete gMinuit;
1239// gMinuit = new TMinuit;
1240// }
4292ca60 1241 AliUnfolding::fgCorrelationMatrix = 0;
1242 AliUnfolding::fgCorrelationMatrixSquared = 0;
1243 AliUnfolding::fgCorrelationCovarianceMatrix = 0;
1244 AliUnfolding::fgCurrentESDVector = 0;
1245 AliUnfolding::fgEntropyAPriori = 0;
1246 AliUnfolding::fgEfficiency = 0;
1247 AliUnfolding::fgUnfoldedAxis = 0;
1248 AliUnfolding::fgMeasuredAxis = 0;
1249 AliUnfolding::fgFitFunction = 0;
1250 AliUnfolding::fgMaxInput = -1;
1251 AliUnfolding::fgMaxParams = -1;
1252 AliUnfolding::fgOverflowBinLimit = -1;
1253 AliUnfolding::fgRegularizationWeight = 10000;
1254 AliUnfolding::fgSkipBinsBegin = 0;
1255 AliUnfolding::fgMinuitStepSize = 0.1;
1256 AliUnfolding::fgMinuitPrecision = 1e-6;
1257 AliUnfolding::fgMinuitMaxIterations = 1000000;
1258 AliUnfolding::fgMinuitStrategy = 1.;
1259 AliUnfolding::fgMinimumInitialValue = kFALSE;
1260 AliUnfolding::fgMinimumInitialValueFix = -1;
1261 AliUnfolding::fgNormalizeInput = kFALSE;
1262 AliUnfolding::fgNotFoundEvents = 0;
1263 AliUnfolding::fgSkipBin0InChi2 = kFALSE;
1264 AliUnfolding::fgBayesianSmoothing = 1;
1265 AliUnfolding::fgBayesianIterations = 10;
1266 AliUnfolding::fgDebug = kFALSE;
1267 AliUnfolding::fgCallCount = 0;
1268 AliUnfolding::fgPowern = 5;
1269 AliUnfolding::fChi2FromFit = 0.;
1270 AliUnfolding::fPenaltyVal = 0.;
1271 AliUnfolding::fAvgResidual = 0.;
1272 AliUnfolding::fgPrintChi2Details = 0;
1273 AliUnfolding::fgCanvas = 0;
1274 AliUnfolding::fghUnfolded = 0;
1275 AliUnfolding::fghCorrelation = 0;
1276 AliUnfolding::fghEfficiency = 0;
1277 AliUnfolding::fghMeasured = 0;
1278 AliUnfolding::SetMinuitStepSize(1.);
1279 AliUnfolding::SetMinuitPrecision(1e-6);
1280 AliUnfolding::SetMinuitMaxIterations(100000);
1281 AliUnfolding::SetMinuitStrategy(2.);
1282 AliUnfolding::SetDebug(1);
1283}
1284//_____________________________________________________________________________
ef12d5a5 1285TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) {
4292ca60 1286 // protect heap by adding unique qualifier to name
1287 if(!protect) return 0x0;
1288 TH1D* p = (TH1D*)protect->Clone();
ef12d5a5 1289 TString tempString(fActiveString);
1290 tempString+=suffix;
1291 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1292 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 1293 if(kill) delete protect;
1294 return p;
1295}
1296//_____________________________________________________________________________
ef12d5a5 1297TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) {
4292ca60 1298 // protect heap by adding unique qualifier to name
1299 if(!protect) return 0x0;
1300 TH2D* p = (TH2D*)protect->Clone();
ef12d5a5 1301 TString tempString(fActiveString);
1302 tempString+=suffix;
1303 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1304 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 1305 if(kill) delete protect;
1306 return p;
1307}
1308//_____________________________________________________________________________
ef12d5a5 1309TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) {
4292ca60 1310 // protect heap by adding unique qualifier to name
1311 if(!protect) return 0x0;
1312 TGraphErrors* p = (TGraphErrors*)protect->Clone();
ef12d5a5 1313 TString tempString(fActiveString);
1314 tempString+=suffix;
1315 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
1316 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4292ca60 1317 if(kill) delete protect;
1318 return p;
1319}
1320//_____________________________________________________________________________