1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Author: Redmer Alexander Bertens, Utrecht University, Utrecht, Netherlands
17 // (rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl)
19 // Tools class for Jet Flow Analysis, replaces 'extractJetFlow.C' macro
21 // The task uses input from two analysis tasks:
22 // - $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetV2.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
28 // libraries must be present on the system
29 // ( see http://hepunx.rl.ac.uk/~adye/software/unfold/RooUnfold.html ).
31 // The weak spot of this class is the function PrepareForUnfolding, which will read
32 // output from two output files and expects histograms with certain names and binning.
33 // Unfolding methods itself are general and should be able to handle any input, therefore one
34 // can forgo the PrepareForUnfolding method, and supply necessary input information via the
35 // SetRawInput() method
37 // to see an example of how to use this class, see $ALICE_ROOT/PWGCF/FLOW/macros/jetFlowTools.C
43 #include "TGraphErrors.h"
44 #include "TGraphAsymmErrors.h"
51 #include "TVirtualFitter.h"
57 #include "TVirtualFitter.h"
58 #include "TFitResultPtr.h"
59 #include "Minuit2/Minuit2Minimizer.h"
60 #include "Math/Functor.h"
62 #include "AliUnfolding.h"
63 #include "AliAnaChargedJetResponseMaker.h"
65 #include "AliJetFlowTools.h"
66 // roo unfold includes (make sure you have these available on your system)
67 #include "RooUnfold.h"
68 #include "RooUnfoldResponse.h"
69 #include "RooUnfoldSvd.h"
70 #include "RooUnfoldBayes.h"
71 #include "TSVDUnfold.h"
74 //_____________________________________________________________________________
75 AliJetFlowTools::AliJetFlowTools() :
76 fResponseMaker (new AliAnaChargedJetResponseMaker()),
81 fPower (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
86 fRefreshInput (kTRUE),
87 fOutputFileName ("UnfoldedSpectra.root"),
89 fCentralityArray (0x0),
90 fMergeBinsArray (0x0),
91 fCentralityWeights (0x0),
92 fDetectorResponse (0x0),
98 fBayesianSmoothIn (0.),
99 fBayesianSmoothOut (0.),
100 fAvoidRoundingError (kFALSE),
101 fUnfoldingAlgorithm (kChi2),
102 fPrior (kPriorMeasured),
106 fBinsTruePrior (0x0),
113 fNormalizeSpectra (kFALSE),
114 fSmoothenPrior (kFALSE),
118 fSmoothenCounts (kTRUE),
120 fRawInputProvided (kFALSE),
121 fEventPlaneRes (.63),
122 fUseDetectorResponse(kTRUE),
123 fUseDptResponse (kTRUE),
125 fDphiUnfolding (kTRUE),
126 fDphiDptUnfolding (kFALSE),
128 fTitleFontSize (-999.),
129 fRMSSpectrumIn (0x0),
130 fRMSSpectrumOut (0x0),
133 fDeltaPtDeltaPhi (0x0),
134 fJetPtDeltaPhi (0x0),
141 fFullResponseIn (0x0),
142 fFullResponseOut (0x0) { // class constructor
143 // create response maker weight function (tuned to PYTHIA spectrum)
144 fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0, 200));
145 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
147 //_____________________________________________________________________________
148 void AliJetFlowTools::Make() {
149 // core function of the class
150 if(fDphiDptUnfolding) {
151 // to extract the yield as function of Dphi, Dpt - experimental
155 // 1) rebin the raw output of the jet task to the desired binnings
156 // 2) calls the unfolding routine
157 // 3) writes output to file
158 // can be repeated multiple times with different configurations
160 // 1) manipulation of input histograms
161 // check if the input variables are present
163 if(!PrepareForUnfolding()) {
164 printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
168 // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
169 // parts of the spectrum can end up in over or underflow bins
171 // if bootstrap mode is kTRUE, resample the underlying distributions
172 // FIXME think about resampling the rebinned results or raw results, could lead to difference
173 // in smoothness of tail of spectrum (which is probably not used in any case, but still ... )
176 // resample but leave original spectra intact for the next unfolding round
177 fSpectrumIn = reinterpret_cast<TH1D*>(Bootstrap(fSpectrumIn, kFALSE));
178 fSpectrumOut = reinterpret_cast<TH1D*>(Bootstrap(fSpectrumOut, kFALSE));
181 TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
182 TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec, TString("resized_out_"), kFALSE);
185 measuredJetSpectrumIn = reinterpret_cast<TH1D*>(Bootstrap(measuredJetSpectrumIn, kFALSE));
186 measuredJetSpectrumOut = reinterpret_cast<TH1D*>(Bootstrap(measuredJetSpectrumOut, kFALSE));
188 // for now do it BEFORE as after gives an issue in Rebin function (counts are wrong)
192 // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
193 // the template will be used as a prior for the chi2 unfolding
194 TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
195 TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
196 // get the full response matrix from the dpt and the detector response
197 fDetectorResponse = NormalizeTH2D(fDetectorResponse);
198 // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
199 // so that unfolding should return the initial spectrum
201 if(fUseDptResponse && fUseDetectorResponse) {
202 fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
203 fFullResponseOut = MatrixMultiplication(fDptOut, fDetectorResponse);
204 } else if (fUseDptResponse && !fUseDetectorResponse) {
205 fFullResponseIn = fDptIn;
206 fFullResponseOut = fDptOut;
207 } else if (!fUseDptResponse && fUseDetectorResponse) {
208 fFullResponseIn = fDetectorResponse;
209 fFullResponseOut = fDetectorResponse;
210 } else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) {
211 printf(" > No response, exiting ! < \n" );
215 fFullResponseIn = GetUnityResponse(fBinsTrue, fBinsRec, TString("in"));
216 fFullResponseOut = GetUnityResponse(fBinsTrue, fBinsRec, TString("out"));
218 // normalize each slide of the response to one
219 NormalizeTH2D(fFullResponseIn);
220 NormalizeTH2D(fFullResponseOut);
221 // resize to desired binning scheme
222 TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
223 TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
224 // get the kinematic efficiency
225 TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
226 kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
227 TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
228 kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
229 // suppress the errors
230 for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
231 kinematicEfficiencyIn->SetBinError(1+i, 0.);
232 kinematicEfficiencyOut->SetBinError(1+i, 0.);
234 TH1D* jetFindingEfficiency(0x0);
236 jetFindingEfficiency = ProtectHeap(fJetFindingEff);
237 jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
238 jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
240 // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
241 TH1D* unfoldedJetSpectrumIn(0x0);
242 TH1D* unfoldedJetSpectrumOut(0x0);
243 fActiveDir->cd(); // select active dir
244 TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
245 dirIn->cd(); // select inplane subdir
246 // do the inplane unfolding
247 unfoldedJetSpectrumIn = UnfoldWrapper(
248 measuredJetSpectrumIn,
250 kinematicEfficiencyIn,
251 measuredJetSpectrumTrueBinsIn,
253 jetFindingEfficiency);
254 resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
255 resizedResponseIn->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
256 resizedResponseIn->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
257 resizedResponseIn = ProtectHeap(resizedResponseIn);
258 resizedResponseIn->Write();
259 kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
260 kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
261 kinematicEfficiencyIn->Write();
262 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
263 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
264 fDetectorResponse->Write();
265 // optional histograms
267 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
268 fSpectrumIn->Write();
269 fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
271 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
273 fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
274 fFullResponseIn->Write();
278 TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
280 // do the out of plane unfolding
281 unfoldedJetSpectrumOut = UnfoldWrapper(
282 measuredJetSpectrumOut,
284 kinematicEfficiencyOut,
285 measuredJetSpectrumTrueBinsOut,
287 jetFindingEfficiency);
288 resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
289 resizedResponseOut->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
290 resizedResponseOut->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
291 resizedResponseOut = ProtectHeap(resizedResponseOut);
292 resizedResponseOut->Write();
293 kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
294 kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
295 kinematicEfficiencyOut->Write();
296 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
297 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
298 fDetectorResponse->Write();
299 if(jetFindingEfficiency) jetFindingEfficiency->Write();
300 // optional histograms
302 fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
303 fSpectrumOut->Write();
304 fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
305 fDptOutDist->Write();
306 fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
308 fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
309 fFullResponseOut->Write();
312 // write general output histograms to file
314 if(unfoldedJetSpectrumIn && unfoldedJetSpectrumOut && unfoldedJetSpectrumIn && unfoldedJetSpectrumOut) {
315 TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out")));
317 ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
318 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
319 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
320 ratio = ProtectHeap(ratio);
322 // write histo values to RMS files if both routines converged
323 // input values are weighted by their uncertainty
324 for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
325 if(unfoldedJetSpectrumIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumIn->GetBinError(i+1), 2.));
326 if(unfoldedJetSpectrumOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), unfoldedJetSpectrumOut->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumOut->GetBinError(i+1), 2.));
327 if(unfoldedJetSpectrumOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1) / unfoldedJetSpectrumOut->GetBinContent(i+1));
330 TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
332 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
333 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
334 v2->GetYaxis()->SetTitle("v_{2}");
335 v2 = ProtectHeap(v2);
338 } else if (unfoldedJetSpectrumOut && unfoldedJetSpectrumIn) {
339 TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
341 ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
342 ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
343 ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
344 ratio = ProtectHeap(ratio);
347 TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
349 v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
350 v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
351 v2->GetYaxis()->SetTitle("v_{2}");
352 v2 = ProtectHeap(v2);
356 } // end of if(fDphiUnfolding)
357 fDeltaPtDeltaPhi->Write();
358 unfoldedJetSpectrumIn->Sumw2();
359 ProtectHeap(unfoldedJetSpectrumIn, kFALSE);
360 unfoldedJetSpectrumIn->Write();
361 unfoldedJetSpectrumOut->Sumw2();
362 ProtectHeap(unfoldedJetSpectrumOut, kFALSE);
363 unfoldedJetSpectrumOut->Write();
364 fJetPtDeltaPhi->Write();
365 // save the current state of the unfolding object
366 SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
367 TH1D* unfoldedJetSpectrumInForSub((TH1D*)unfoldedJetSpectrumIn->Clone("forSubIn"));
368 TH1D* unfoldedJetSpectrumOutForSub((TH1D*)unfoldedJetSpectrumOut->Clone("forSubOut"));
369 unfoldedJetSpectrumInForSub->Add(unfoldedJetSpectrumOutForSub, -1.);
370 unfoldedJetSpectrumInForSub= ProtectHeap(unfoldedJetSpectrumInForSub);
371 unfoldedJetSpectrumInForSub->Write();
374 //_____________________________________________________________________________
375 TH1D* AliJetFlowTools::UnfoldWrapper(
376 const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
377 const TH2D* resizedResponse, // response matrix
378 const TH1D* kinematicEfficiency, // kinematic efficiency
379 const TH1D* measuredJetSpectrumTrueBins, // unfolding template: same binning is pt gen of response
380 const TString suffix, // suffix (in or out of plane)
381 const TH1D* jetFindingEfficiency) // jet finding efficiency
383 // wrapper function to call specific unfolding routine
384 TH1D* (AliJetFlowTools::*myFunction)(const TH1D*, const TH2D*, const TH1D*, const TH1D*, const TString, const TH1D*);
385 // initialize functon pointer
386 if(fUnfoldingAlgorithm == kChi2) myFunction = &AliJetFlowTools::UnfoldSpectrumChi2;
387 else if(fUnfoldingAlgorithm == kBayesian) myFunction = &AliJetFlowTools::UnfoldSpectrumBayesian;
388 else if(fUnfoldingAlgorithm == kBayesianAli) myFunction = &AliJetFlowTools::UnfoldSpectrumBayesianAli;
389 else if(fUnfoldingAlgorithm == kSVD) myFunction = &AliJetFlowTools::UnfoldSpectrumSVD;
390 else if(fUnfoldingAlgorithm == kNone) {
391 TH1D* clone((TH1D*)measuredJetSpectrum->Clone("clone"));
392 clone->SetNameTitle(Form("MeasuredJetSpectrum%s", suffix.Data()), Form("measuredJetSpectrum %s", suffix.Data()));
393 return clone;//RebinTH1D(clone, fBinsTrue, clone->GetName(), kFALSE);
396 // do the actual unfolding with the selected function
397 return (this->*myFunction)( measuredJetSpectrum, resizedResponse, kinematicEfficiency, measuredJetSpectrumTrueBins, suffix, jetFindingEfficiency);
399 //_____________________________________________________________________________
400 TH1D* AliJetFlowTools::UnfoldSpectrumChi2(
401 const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
402 const TH2D* resizedResponse, // response matrix
403 const TH1D* kinematicEfficiency, // kinematic efficiency
404 const TH1D* measuredJetSpectrumTrueBins, // unfolding template: same binning is pt gen of response
405 const TString suffix, // suffix (in or out of plane)
406 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
408 // unfold the spectrum using chi2 minimization
410 // step 0) setup the static members of AliUnfolding
411 ResetAliUnfolding(); // reset from previous iteration
412 // also deletes and re-creates the global TVirtualFitter
413 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
414 if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
415 else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
416 if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
417 else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
418 AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
420 // step 1) clone all input histograms. the histograms are cloned to make sure that the original histograms
421 // stay intact. a local copy of a histogram (which only exists in the scope of this function) is
422 // denoted by the suffix 'Local'
424 // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
425 TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
426 // unfolded local will be filled with the result of the unfolding
427 TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
429 // full response matrix and kinematic efficiency
430 TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
431 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
433 // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
434 TH1D *priorLocal = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("priorLocal_%s", suffix.Data()));
435 // optionally, the prior can be smoothened by extrapolating the spectrum using a power law fit
436 if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
438 // step 2) start the unfolding
439 Int_t status(-1), i(0);
440 while(status < 0 && i < 100) {
441 // i > 0 means that the first iteration didn't converge. in that case, the result of the first
442 // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
443 if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
444 status = AliUnfolding::Unfold(
445 resizedResponseLocal, // response matrix
446 kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
447 measuredJetSpectrumLocal, // measured spectrum
448 priorLocal, // initial conditions (set NULL to use measured spectrum)
449 unfoldedLocal); // results
450 // status holds the minuit fit status (where 0 means convergence)
453 // 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
457 TVirtualFitter *fitter(TVirtualFitter::GetFitter());
458 if(gMinuit) gMinuit->Command("SET COV");
459 TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
460 TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
462 hPearson = new TH2D(*pearson);
463 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
464 hPearson = ProtectHeap(hPearson);
466 if(fMergeBinsArray) unfoldedLocal = MergeSpectrumBins(fMergeBinsArray, unfoldedLocal, hPearson);
469 // step 3) refold the unfolded spectrum and save the ratio measured / refolded
470 TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
471 foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
472 unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
473 TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
475 ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
476 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
477 ratio = ProtectHeap(ratio);
481 // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
482 // objects are cloned using 'ProtectHeap()'
483 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
484 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
485 measuredJetSpectrumLocal->Write();
487 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
488 resizedResponseLocal->Write();
490 unfoldedLocal = ProtectHeap(unfoldedLocal);
491 if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
492 unfoldedLocal->Write();
494 foldedLocal = ProtectHeap(foldedLocal);
495 foldedLocal->Write();
497 priorLocal = ProtectHeap(priorLocal);
500 // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
501 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));
502 fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
503 fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
504 fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
505 fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
506 fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
507 fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
508 fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
509 fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
512 return unfoldedLocal;
514 //_____________________________________________________________________________
515 TH1D* AliJetFlowTools::UnfoldSpectrumSVD(
516 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
517 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
518 const TH1D* kinematicEfficiency, // kinematic efficiency
519 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
520 const TString suffix, // suffix (in, out)
521 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
524 TH1D* priorLocal( GetPrior(
525 measuredJetSpectrum, // jet pt in pt rec bins
526 resizedResponse, // full response matrix, normalized in slides of pt true
527 kinematicEfficiency, // kinematic efficiency
528 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
529 suffix, // suffix (in, out)
530 jetFindingEfficiency)); // jet finding efficiency (optional)
532 printf(" > couldn't find prior ! < \n");
534 } else printf(" 1) retrieved prior \n");
536 // go back to the 'root' directory of this instance of the SVD unfolding routine
537 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
539 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
540 // measured jets in pt rec binning
541 TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
542 // local copie of the response matrix
543 TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
544 // local copy of response matrix, all true slides normalized to 1
545 // this response matrix will eventually be used in the re-folding routine
546 TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
547 resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
548 // kinematic efficiency
549 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
550 // place holder histos
551 TH1D *unfoldedLocalSVD(0x0);
552 TH1D *foldedLocalSVD(0x0);
553 cout << " 2) setup necessary input " << endl;
554 // 3) configure routine
555 RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
556 cout << " step 3) configured routine " << endl;
558 // 4) get transpose matrices
559 // a) get the transpose of the full response matrix
560 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
561 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
562 // normalize it with the prior. this will ensure that high statistics bins will constrain the
563 // end result most strenuously than bins with limited number of counts
564 responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
565 cout << " 4) retrieved first transpose matrix " << endl;
567 // 5) get response for SVD unfolding
568 RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
569 cout << " 5) retrieved roo unfold response object " << endl;
571 // 6) actualy unfolding loop
572 RooUnfoldSvd unfoldSVD(&responseSVD, measuredJetSpectrumLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
573 unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
574 // correct the spectrum for the kinematic efficiency
575 unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
577 // get the pearson coefficients from the covariance matrix
578 TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
579 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
582 hPearson = new TH2D(*pearson);
584 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
585 hPearson = ProtectHeap(hPearson);
587 if(fMergeBinsArray) unfoldedLocalSVD = MergeSpectrumBins(fMergeBinsArray, unfoldedLocalSVD, hPearson);
588 } else return 0x0; // return if unfolding didn't converge
590 // plot singular values and d_i vector
591 TSVDUnfold* svdUnfold(unfoldSVD.Impl());
592 TH1* hSVal(svdUnfold->GetSV());
593 TH1D* hdi(svdUnfold->GetD());
594 hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
595 hSVal->SetXTitle("singular values");
597 hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
598 hdi->SetXTitle("|d_{i}^{kreg}|");
600 cout << " plotted singular values and d_i vector " << endl;
602 // 7) refold the unfolded spectrum
603 foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, resizedResponseLocalNorm, kinematicEfficiencyLocal);
604 TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalSVD, "ratio measured / re-folded", kTRUE));
605 ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
606 ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
607 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
609 cout << " 7) refolded the unfolded spectrum " << endl;
611 // write the measured, unfolded and re-folded spectra to the output directory
612 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
613 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
614 measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
615 measuredJetSpectrumLocal->Write(); // input spectrum
616 unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
617 unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
618 if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
619 unfoldedLocalSVD->Write(); // unfolded spectrum
620 foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
621 foldedLocalSVD = ProtectHeap(foldedLocalSVD);
622 foldedLocalSVD->Write(); // re-folded spectrum
624 // save more general bookkeeeping histograms to the output directory
625 responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
626 responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
627 responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
628 responseMatrixLocalTransposePrior->Write();
629 priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
630 priorLocal->SetXTitle("p_{t} [GeV/c]");
631 priorLocal = ProtectHeap(priorLocal);
633 resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
634 resizedResponseLocalNorm->Write();
637 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));
638 fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
639 fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
642 return unfoldedLocalSVD;
644 //_____________________________________________________________________________
645 TH1D* AliJetFlowTools::UnfoldSpectrumBayesianAli(
646 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
647 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
648 const TH1D* kinematicEfficiency, // kinematic efficiency
649 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
650 const TString suffix, // suffix (in, out)
651 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
653 // unfold the spectrum using the bayesian unfolding impelmented in AliUnfolding
654 // FIXME careful, not tested yet ! (06122013) FIXME
656 // step 0) setup the static members of AliUnfolding
657 ResetAliUnfolding(); // reset from previous iteration
658 // also deletes and re-creates the global TVirtualFitter
659 AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
660 if(!strcmp("in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
661 else if(!strcmp("out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
662 else if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
663 else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
664 AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
666 // 1) get a prior for unfolding and clone all the input histograms
667 TH1D* priorLocal( GetPrior(
668 measuredJetSpectrum, // jet pt in pt rec bins
669 resizedResponse, // full response matrix, normalized in slides of pt true
670 kinematicEfficiency, // kinematic efficiency
671 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
672 suffix, // suffix (in, out)
673 jetFindingEfficiency)); // jet finding efficiency (optional)
675 printf(" > couldn't find prior ! < \n");
677 } else printf(" 1) retrieved prior \n");
678 // switch back to root dir of this unfolding procedure
679 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
681 // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
682 TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
683 // unfolded local will be filled with the result of the unfolding
684 TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
686 // full response matrix and kinematic efficiency
687 TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
688 TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
690 // step 2) start the unfolding
691 Int_t status(-1), i(0);
692 while(status < 0 && i < 100) {
693 // i > 0 means that the first iteration didn't converge. in that case, the result of the first
694 // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
695 if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
696 status = AliUnfolding::Unfold(
697 resizedResponseLocal, // response matrix
698 kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
699 measuredJetSpectrumLocal, // measured spectrum
700 priorLocal, // initial conditions (set NULL to use measured spectrum)
701 unfoldedLocal); // results
702 // status holds the minuit fit status (where 0 means convergence)
705 // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
707 if(status == 0 && gMinuit->fISW[1] == 3) {
708 // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
709 TVirtualFitter *fitter(TVirtualFitter::GetFitter());
710 if(gMinuit) gMinuit->Command("SET COV");
711 TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
712 TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
714 hPearson= new TH2D(*pearson);
715 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
716 hPearson = ProtectHeap(hPearson);
719 if(fMergeBinsArray) unfoldedLocal = MergeSpectrumBins(fMergeBinsArray, unfoldedLocal, hPearson);
721 // step 3) refold the unfolded spectrum and save the ratio measured / refolded
722 TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
723 foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
724 unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
725 TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
727 ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
728 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
729 ratio = ProtectHeap(ratio);
733 // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
734 // objects are cloned using 'ProtectHeap()'
735 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
736 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
737 measuredJetSpectrumLocal->Write();
739 resizedResponseLocal = ProtectHeap(resizedResponseLocal);
740 resizedResponseLocal->Write();
742 unfoldedLocal = ProtectHeap(unfoldedLocal);
743 if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
744 unfoldedLocal->Write();
746 foldedLocal = ProtectHeap(foldedLocal);
747 foldedLocal->Write();
749 priorLocal = ProtectHeap(priorLocal);
752 // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
753 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));
754 fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
755 fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
756 fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
757 fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
758 fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
759 fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
760 fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
761 fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
764 return unfoldedLocal;
766 //_____________________________________________________________________________
767 TH1D* AliJetFlowTools::UnfoldSpectrumBayesian(
768 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
769 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
770 const TH1D* kinematicEfficiency, // kinematic efficiency
771 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
772 const TString suffix, // suffix (in, out)
773 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
775 // use bayesian unfolding from the RooUnfold package to unfold jet spectra
777 // 1) get a prior for unfolding.
778 TH1D* priorLocal( GetPrior(
779 measuredJetSpectrum, // jet pt in pt rec bins
780 resizedResponse, // full response matrix, normalized in slides of pt true
781 kinematicEfficiency, // kinematic efficiency
782 measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
783 suffix, // suffix (in, out)
784 jetFindingEfficiency)); // jet finding efficiency (optional)
786 printf(" > couldn't find prior ! < \n");
788 } else printf(" 1) retrieved prior \n");
789 (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
791 // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
792 // measured jets in pt rec binning
793 TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
794 // local copie of the response matrix
795 TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
796 // local copy of response matrix, all true slides normalized to 1
797 // this response matrix will eventually be used in the re-folding routine
798 TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
799 resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
800 // kinematic efficiency
801 TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
802 // place holder histos
803 TH1D *unfoldedLocalBayes(0x0);
804 TH1D *foldedLocalBayes(0x0);
805 cout << " 2) setup necessary input " << endl;
806 // 4) get transpose matrices
807 // a) get the transpose of the full response matrix
808 TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
809 responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
810 // normalize it with the prior. this will ensure that high statistics bins will constrain the
811 // end result most strenuously than bins with limited number of counts
812 responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
813 // 3) get response for Bayesian unfolding
814 RooUnfoldResponse responseBayes(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedBayes_%s", suffix.Data()), Form("respCombinedBayes_%s", suffix.Data()));
816 // 4) actualy unfolding loop
817 RooUnfoldBayes unfoldBayes(&responseBayes, measuredJetSpectrumLocal, (!strcmp("in", suffix.Data())) ? fBayesianIterIn : fBayesianIterOut);
818 RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
819 unfoldedLocalBayes = (TH1D*)unfoldBayes.Hreco(errorTreatment);
820 // correct the spectrum for the kinematic efficiency
821 unfoldedLocalBayes->Divide(kinematicEfficiencyLocal);
822 // get the pearson coefficients from the covariance matrix
823 TMatrixD covarianceMatrix = unfoldBayes.Ereco(errorTreatment);
824 TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
827 hPearson = new TH2D(*pearson);
829 hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
830 hPearson = ProtectHeap(hPearson);
832 if(fMergeBinsArray) unfoldedLocalBayes = MergeSpectrumBins(fMergeBinsArray, unfoldedLocalBayes, hPearson);
833 } else return 0x0; // return if unfolding didn't converge
835 // 5) refold the unfolded spectrum
836 foldedLocalBayes = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalBayes, resizedResponseLocalNorm, kinematicEfficiencyLocal);
837 TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalBayes, "ratio measured / re-folded", kTRUE));
838 ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
839 ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
840 ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
842 cout << " 7) refolded the unfolded spectrum " << endl;
844 // write the measured, unfolded and re-folded spectra to the output directory
845 measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
846 measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
847 measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
848 measuredJetSpectrumLocal->Write(); // input spectrum
849 unfoldedLocalBayes->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
850 unfoldedLocalBayes = ProtectHeap(unfoldedLocalBayes);
851 if(jetFindingEfficiency) unfoldedLocalBayes->Divide(jetFindingEfficiency);
852 unfoldedLocalBayes->Write(); // unfolded spectrum
853 foldedLocalBayes->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
854 foldedLocalBayes = ProtectHeap(foldedLocalBayes);
855 foldedLocalBayes->Write(); // re-folded spectrum
857 // save more general bookkeeeping histograms to the output directory
858 responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
859 responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
860 responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
861 responseMatrixLocalTransposePrior->Write();
862 priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
863 priorLocal->SetXTitle("p_{t} [GeV/c]");
864 priorLocal = ProtectHeap(priorLocal);
866 resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
867 resizedResponseLocalNorm->Write();
870 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));
871 fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fBayesianIterIn : fBayesianIterOut);
872 fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fBayesianIterIn" : "fBayesianIterOut");
875 return unfoldedLocalBayes;
877 //_____________________________________________________________________________
878 Bool_t AliJetFlowTools::PrepareForUnfolding()
880 // prepare for unfolding
881 if(fRawInputProvided) return kTRUE;
883 printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
886 if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
887 // check if the pt bin for true and rec have been set
888 if(!fBinsTrue || !fBinsRec) {
889 printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
892 if(!fRMSSpectrumIn && fDphiUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
893 // procedures, these profiles will be nonsensical, user is responsible
894 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
895 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
896 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
899 // clear minuit state to avoid constraining the fit with the results of the previous iteration
900 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
902 // extract the spectra
903 TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
904 if(fRho0) spectrumName = Form("fHistJetPsi2PtRho0_%i", fCentralityArray->At(0));
905 if(!fInputList->FindObject(spectrumName.Data())) {
906 printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
910 // get the first scaled spectrum
911 fJetPtDeltaPhi = (TH2D*)fInputList->FindObject(spectrumName.Data());
912 // clone the spectrum on the heap. this is necessary since scale or add change the
913 // contents of the original histogram
914 fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
915 fJetPtDeltaPhi->Scale(fCentralityWeights->At(0));
916 printf("Extracted %s wight weight %.2f \n", spectrumName.Data(), fCentralityWeights->At(0));
917 // merge subsequent bins (if any)
918 for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
919 spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i));
920 printf( " Merging with %s with weight %.4f \n", spectrumName.Data(), fCentralityWeights->At(i));
921 fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())), fCentralityWeights->At(i));
925 if(!fDphiUnfolding) {
926 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e");
927 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40, "e");
929 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10, "e");
930 fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40, "e"));
931 fSpectrumIn = ProtectHeap(fSpectrumIn);
932 // out of plane spectrum
933 fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30, "e");
934 fSpectrumOut = ProtectHeap(fSpectrumOut);
936 // normalize spectra to event count if requested
937 if(fNormalizeSpectra) {
938 TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(0))));
939 rho->Scale(fCentralityWeights->At(0));
940 for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
941 rho->Add((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(i))), fCentralityWeights->At(i));
944 Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
945 if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
946 if(fEventCount > 0) {
947 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
948 fSpectrumOut->Sumw2();
950 Double_t error(0); // lots of issues with the errors here ...
951 for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
952 pt = fSpectrumIn->GetBinContent(1+i)/fEventCount; // normalized count
953 error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
954 fSpectrumIn->SetBinContent(1+i, pt);
955 if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
956 if(error > 0) fSpectrumIn->SetBinError(1+i, error);
957 else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
959 for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
960 pt = fSpectrumOut->GetBinContent(1+i)/fEventCount; // normalized count
961 error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
962 fSpectrumOut->SetBinContent(1+i, pt);
963 if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
964 if(error > 0) fSpectrumOut->SetBinError(1+i, error);
965 else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
968 if(normalizeToFullSpectrum) fEventCount = -1;
970 // extract the delta pt matrices
971 TString deltaptName("");
973 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
975 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJRho0_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2Rho0_%i", fCentralityArray->At(0));
977 fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
978 if(!fDeltaPtDeltaPhi) {
979 printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
980 printf(" > may be ok, depending no what you want to do < \n");
981 fRefreshInput = kTRUE;
985 // clone the distribution on the heap and if requested merge with other centralities
986 fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
987 fDeltaPtDeltaPhi->Scale(fCentralityWeights->At(0));
988 printf("Extracted %s with weight %.2f \n", deltaptName.Data(), fCentralityWeights->At(0));
989 for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
990 deltaptName = (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i));
991 printf(" Merging with %s with weight %.4f \n", deltaptName.Data(), fCentralityWeights->At(i));
992 fDeltaPtDeltaPhi->Add((TH2D*)fInputList->FindObject(deltaptName.Data()), fCentralityWeights->At(i));
995 // in plane delta pt distribution
996 if(!fDphiUnfolding) {
997 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e");
998 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40, "e");
1000 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10, "e");
1001 fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40, "e"));
1002 // out of plane delta pt distribution
1003 fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30, "e");
1004 fDptInDist = ProtectHeap(fDptInDist);
1005 fDptOutDist = ProtectHeap(fDptOutDist);
1006 // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
1009 // create a rec - true smeared response matrix
1010 TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1011 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
1012 Bool_t skip = kFALSE;
1013 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
1014 (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1015 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1018 TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
1019 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
1020 Bool_t skip = kFALSE;
1021 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
1022 (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
1023 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1026 fDptIn = new TH2D(*rfIn);
1027 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
1028 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1029 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1030 fDptIn = ProtectHeap(fDptIn);
1031 fDptOut = new TH2D(*rfOut);
1032 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
1033 fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1034 fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1035 fDptOut = ProtectHeap(fDptOut);
1037 fRefreshInput = kTRUE; // force cloning of the input
1040 //_____________________________________________________________________________
1041 Bool_t AliJetFlowTools::PrepareForUnfolding(Int_t low, Int_t up) {
1042 // prepare for unfoldingUA - more robust method to extract input spectra from file
1043 // will replace PrepareForUnfolding eventually (09012014)
1045 printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
1048 if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
1049 // check if the pt bin for true and rec have been set
1050 if(!fBinsTrue || !fBinsRec) {
1051 printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
1055 // clear minuit state to avoid constraining the fit with the results of the previous iteration
1056 for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
1058 // extract the spectra
1059 TString spectrumName(Form("fHistJetPsi2Pt_%i", fCentralityArray->At(0)));
1060 if(fRho0) spectrumName = Form("fHistJetPsi2PtRho0_%i", fCentralityArray->At(0));
1061 fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
1062 if(!fJetPtDeltaPhi) {
1063 printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
1066 if(fCentralityArray) {
1067 for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1068 spectrumName = Form("fHistJetPsi2Pt_%i", fCentralityArray->At(i));
1069 fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())));
1072 fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
1073 // in plane spectrum
1074 fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e");
1075 // extract the delta pt matrices
1076 TString deltaptName("");
1078 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(0));
1080 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJRho0_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi2Rho0_%i", fCentralityArray->At(0));
1082 fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1083 if(!fDeltaPtDeltaPhi) {
1084 printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
1085 printf(" > may be ok, depending no what you want to do < \n");
1086 fRefreshInput = kTRUE;
1089 if(fCentralityArray) {
1090 for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1091 deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi2ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi2_%i", fCentralityArray->At(i));
1092 fDeltaPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(deltaptName.Data())));
1096 fDeltaPtDeltaPhi = ProtectHeap(fDeltaPtDeltaPhi, kFALSE);
1097 // in plane delta pt distribution
1098 fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), low, up, "e");
1099 // create a rec - true smeared response matrix
1100 TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1101 for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
1102 Bool_t skip = kFALSE;
1103 for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
1104 (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1105 if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1108 fDptIn = new TH2D(*rfIn);
1109 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
1110 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1111 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1112 fDptIn = ProtectHeap(fDptIn);
1116 //_____________________________________________________________________________
1117 TH1D* AliJetFlowTools::GetPrior(
1118 const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
1119 const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
1120 const TH1D* kinematicEfficiency, // kinematic efficiency
1121 const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
1122 const TString suffix, // suffix (in, out)
1123 const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
1125 // 1) get a prior for unfolding.
1126 // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
1127 TH1D* unfolded(0x0);
1128 TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
1130 switch (fPrior) { // select the prior for unfolding
1133 printf("> GetPrior:: FATAL ERROR! TF1 prior requested but prior has not been set ! < \n");
1138 case kPriorPythia : {
1140 printf("> GetPrior:: FATAL ERROR! pythia prior requested but prior has not been set ! < \n");
1143 // rebin the given prior to the true spectrum (creates a new histo)
1144 return RebinTH1D(fPriorUser, fBinsTrue, Form("kPriorPythia_%s", suffix.Data()), kFALSE);
1147 TArrayI* placeHolder(0x0);
1148 if(fMergeBinsArray) {
1149 placeHolder = fMergeBinsArray;
1150 fMergeBinsArray = 0x0;
1152 if(fBinsTruePrior && fBinsRecPrior) { // if set, use different binning for the prior
1153 TArrayD* tempArrayTrue(fBinsTrue); // temporarily cache the original binning
1154 fBinsTrue = fBinsTruePrior; // switch binning schemes (will be used in UnfoldSpectrumChi2())
1155 TArrayD* tempArrayRec(fBinsRec);
1156 fBinsRec = fBinsRecPrior;
1157 // for the prior, do not re-bin the output
1158 TH1D* measuredJetSpectrumChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
1159 TH1D* measuredJetSpectrumTrueBinsChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
1160 TH2D* resizedResponseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
1161 TH1D* kinematicEfficiencyChi2(resizedResponseChi2->ProjectionX());
1162 kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
1163 for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
1164 unfolded= UnfoldSpectrumChi2(
1165 measuredJetSpectrumChi2,
1166 resizedResponseChi2,
1167 kinematicEfficiencyChi2,
1168 measuredJetSpectrumTrueBinsChi2, // prior for chi2 unfolding (measured)
1169 TString(Form("prior_%s", suffix.Data())));
1171 printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1172 printf(" probably Chi2 unfolding did not converge < \n");
1173 if(placeHolder) fMergeBinsArray = placeHolder;
1176 fBinsTrue = tempArrayTrue; // reset bins borders
1177 fBinsRec = tempArrayRec;
1178 // if the chi2 prior has a different binning, rebin to the true binning for the unfolding
1179 unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data()))); // rebin unfolded
1181 unfolded = UnfoldSpectrumChi2(
1182 measuredJetSpectrum,
1184 kinematicEfficiency,
1185 measuredJetSpectrumTrueBins, // prior for chi2 unfolding (measured)
1186 TString(Form("prior_%s", suffix.Data())));
1188 printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1189 printf(" probably Chi2 unfolding did not converge < \n");
1190 if(placeHolder) fMergeBinsArray = placeHolder;
1193 if(placeHolder) fMergeBinsArray = placeHolder;
1198 case kPriorMeasured : {
1199 unfolded = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("kPriorMeasured_%s", suffix.Data())); // copy template to unfolded to use as prior
1203 // it can be important that the prior is smooth, this can be achieved by
1204 // extrapolating the spectrum with a fitted power law when bins are sparsely filed
1205 if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
1206 TH1D *priorLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
1207 if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
1210 //_____________________________________________________________________________
1211 TH1D* AliJetFlowTools::ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix) {
1212 // resize the x-axis of a th1d
1214 printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1217 // see how many bins we need to copy
1218 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);
1219 // low is the bin number of the first new bin
1220 Int_t l = histo->GetXaxis()->FindBin(low);
1222 Double_t _x(0), _xx(0);
1223 for(Int_t i(0); i < up-low; i++) {
1224 _x = histo->GetBinContent(l+i);
1225 _xx=histo->GetBinError(l+i);
1226 resized->SetBinContent(i+1, _x);
1227 resized->SetBinError(i+1, _xx);
1231 //_____________________________________________________________________________
1232 TH2D* AliJetFlowTools::ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix) {
1233 // resize the y-axis of a th2d
1235 printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1238 // see how many bins we need to copy
1239 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());
1240 // assume only the y-axis has changed
1241 // low is the bin number of the first new bin
1242 Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1244 Double_t _x(0), _xx(0);
1245 for(Int_t i(0); i < x->GetSize(); i++) {
1246 for(Int_t j(0); j < y->GetSize(); j++) {
1247 _x = histo->GetBinContent(i, low+j);
1248 _xx=histo->GetBinError(i, low+1+j);
1249 resized->SetBinContent(i, j, _x);
1250 resized->SetBinError(i, j, _xx);
1255 //_____________________________________________________________________________
1256 TH2D* AliJetFlowTools::NormalizeTH2D(TH2D* histo, Bool_t noError) {
1257 // general method to normalize all vertical slices of a th2 to unity
1258 // i.e. get a probability matrix
1260 printf(" > NormalizeTH2D:: NULL pointer passed, returning NULL < \n");
1263 Int_t binsX = histo->GetXaxis()->GetNbins();
1264 Int_t binsY = histo->GetYaxis()->GetNbins();
1266 // normalize all slices in x
1267 for(Int_t i(0); i < binsX; i++) { // for each vertical slice
1268 Double_t weight = 0;
1269 for(Int_t j(0); j < binsY; j++) { // loop over all the horizontal components
1270 weight+=histo->GetBinContent(i+1, j+1);
1271 } // now we know the total weight
1272 for(Int_t j(0); j < binsY; j++) {
1273 if (weight <= 0 ) continue;
1274 histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
1275 if(noError) histo->SetBinError( 1+i, j+1, 0.);
1276 else histo->SetBinError( 1+i, j+1, histo->GetBinError( 1+i, j+1)/weight);
1281 //_____________________________________________________________________________
1282 TH1* AliJetFlowTools::Bootstrap(TH1* hist, Bool_t kill) {
1284 // the returned histogram is new, the original is deleted from the heap if kill is true
1286 printf(" > Bootstrap:: fatal error,NULL pointer passed! \n");
1289 else printf(" > Bootstrap:: resampling, may take some time \n");
1290 // clone input histo
1291 TH1* bootstrapped((TH1*)(hist->Clone(Form("%s_bootstrapped", hist->GetName()))));
1292 bootstrapped->Reset();
1294 /* OLD method - slightly underestimates fluctuations
1295 // reset the content
1296 bootstrapped->Reset();
1297 // resample the input histogram
1298 for(Int_t i(0); i < hist->GetEntries(); i++) bootstrapped->Fill(hist->GetRandom()); */
1301 Double_t mean(0), sigma(0);
1302 Int_t sampledMean(0), entries(0);
1304 for(Int_t i(0); i < hist->GetXaxis()->GetNbins(); i++) {
1305 // for each bin, get the value
1306 mean = hist->GetBinContent(i+1);
1307 sigma = hist->GetBinError(i+1);
1309 sampledMean = TMath::Nint(gRandom->Gaus(mean, sigma));
1310 printf(" sampled %i from original number %.2f \n", sampledMean, mean);
1311 // set the new bin content
1312 bootstrapped->SetBinContent(i+1, sampledMean);
1313 if(sampledMean > 0) bootstrapped->SetBinError(i+1, TMath::Sqrt(sampledMean));
1314 entries += sampledMean;
1316 printf(" Done bootstrapping, setting number of entries to %i \n", entries);
1317 bootstrapped->SetEntries((double)entries);
1320 // if requested kill input histo
1321 if(kill) delete hist;
1322 // return resampled histogram
1323 return bootstrapped;
1325 //_____________________________________________________________________________
1326 TH1D* AliJetFlowTools::RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix, Bool_t kill) {
1327 // return a TH1D with the supplied histogram rebinned to the supplied bins
1328 // the returned histogram is new, the original is deleted from the heap if kill is true
1329 if(!histo || !bins) {
1330 printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
1333 // create the output histo
1334 TString name = histo->GetName();
1337 TH1D* rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1338 for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
1339 // loop over the bins of the old histo and fill the new one with its data
1340 rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
1342 if(kill) delete histo;
1345 //_____________________________________________________________________________
1346 TH2D* AliJetFlowTools::RebinTH2D(TH2D* rebinMe, TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
1347 // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
1348 // not static as it is just a wrapper for the response maker object
1349 if(!fResponseMaker || !binsTrue || !binsRec) {
1350 printf(" > RebinTH2D:: function called with NULL arguments < \n");
1353 TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1354 return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
1356 //_____________________________________________________________________________
1357 TH2D* AliJetFlowTools::MatrixMultiplication(TH2D* a, TH2D* b, TString name)
1359 // multiply two matrices
1360 if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1361 TH2D* c = (TH2D*)a->Clone("c");
1362 for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1363 for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1365 for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1367 val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1369 c->SetBinContent(x2, y1, val);
1370 c->SetBinError(x2, y1, 0.);
1373 if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1376 //_____________________________________________________________________________
1377 TH1D* AliJetFlowTools::NormalizeTH1D(TH1D* histo, Double_t scale)
1379 // normalize a th1d to a certain scale
1381 Double_t integral = histo->Integral()*scale;
1382 if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1383 else if (scale != 1.) histo->Scale(1./scale, "width");
1384 else printf(" > Histogram integral < 0, cannot normalize \n");
1387 //_____________________________________________________________________________
1388 TH1D* AliJetFlowTools::MergeSpectrumBins(TArrayI* bins, TH1D* spectrum, TH2D* corr)
1390 // merge a spectrum histogram taking into account the correlation terms
1391 if(!(bins&&spectrum)) {
1392 printf(" > NULL pointer passed as argument in MergeSpectrumBins ! < \n");
1395 Double_t sum(0), error(0), pearson(0);
1396 // take the sum of the bin content
1397 sum += spectrum->GetBinContent(bins->At(0));
1398 sum += spectrum->GetBinContent(bins->At(1));
1399 // quadratically sum the errors
1400 error += TMath::Power(spectrum->GetBinError(bins->At(0)), 2);
1401 error += TMath::Power(spectrum->GetBinError(bins->At(1)), 2);
1402 // add the covariance term
1403 pearson = corr->GetBinContent(bins->At(0), bins->At(1));
1405 printf(" > PANIC ! something is wrong with the covariance matrix, assuming full correlation ! < \n ");
1408 error += 2.*spectrum->GetBinError(bins->At(0))*spectrum->GetBinError(bins->At(1))*pearson;
1409 spectrum->SetBinContent(1, sum);
1410 if(error <= 0) return spectrum;
1411 spectrum->SetBinError(1, TMath::Sqrt(error));
1412 printf(" > sum is %.2f \t error is %.8f < \n", sum, error);
1413 printf(" > REPLACING BIN CONTENT OF FIRST BIN (%i) WITH MERGED BINS (%i) and (%i) < \n", 1, bins->At(0), bins->At(1));
1416 //_____________________________________________________________________________
1417 TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
1419 // Calculate pearson coefficients from covariance matrix
1420 TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1421 Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
1422 Double_t pearson(0.);
1423 if(nrows==0 && ncols==0) return 0x0;
1424 for(Int_t row = 0; row < nrows; row++) {
1425 for(Int_t col = 0; col<ncols; col++) {
1426 if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1427 (*pearsonCoefficients)(row,col) = pearson;
1430 return pearsonCoefficients;
1432 //_____________________________________________________________________________
1433 TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
1434 // smoothen the spectrum using a user defined function
1435 // returns a clone of the original spectrum if fitting failed
1436 // if kill is kTRUE the input spectrum will be deleted from the heap
1437 // if 'count' is selected, bins are filled with integers (necessary if the
1438 // histogram is interpreted in a routine which accepts only counts)
1439 if(!spectrum || !function) return 0x0;
1440 if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
1441 TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1442 temp->Sumw2(); // if already called on the original, this will give off a warning but do nothing
1443 TFitResultPtr r = temp->Fit(function, "", "", min, max);
1444 if((int)r == 0) { // MINUIT status
1445 for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1446 if(temp->GetBinCenter(i) > start) { // from this pt value use extrapolation
1447 (counts) ? temp->SetBinContent(i, (int)(function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i))) : temp->SetBinContent(i, function->Integral(temp->GetXaxis()->GetBinLowEdge(i),temp->GetXaxis()->GetBinUpEdge(i))/temp->GetXaxis()->GetBinWidth(i));
1448 if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1452 if(kill) delete spectrum;
1455 //_____________________________________________________________________________
1456 void AliJetFlowTools::Style(Bool_t legacy)
1458 // set global style for your current aliroot session
1460 // legacy style is pleasing to the eye, default is the formal ALICE style
1462 gStyle->SetCanvasColor(-1);
1463 gStyle->SetPadColor(-1);
1464 gStyle->SetFrameFillColor(-1);
1465 gStyle->SetHistFillColor(-1);
1466 gStyle->SetTitleFillColor(-1);
1467 gStyle->SetFillColor(-1);
1468 gStyle->SetFillStyle(4000);
1469 gStyle->SetStatStyle(0);
1470 gStyle->SetTitleStyle(0);
1471 gStyle->SetCanvasBorderSize(0);
1472 gStyle->SetFrameBorderSize(0);
1473 gStyle->SetLegendBorderSize(0);
1474 gStyle->SetStatBorderSize(0);
1475 gStyle->SetTitleBorderSize(0);
1477 gStyle->Reset("Plain");
1478 gStyle->SetOptTitle(0);
1479 gStyle->SetOptStat(0);
1480 gStyle->SetPalette(1);
1481 gStyle->SetCanvasColor(10);
1482 gStyle->SetCanvasBorderMode(0);
1483 gStyle->SetFrameLineWidth(1);
1484 gStyle->SetFrameFillColor(kWhite);
1485 gStyle->SetPadColor(10);
1486 gStyle->SetPadTickX(1);
1487 gStyle->SetPadTickY(1);
1488 gStyle->SetPadBottomMargin(0.15);
1489 gStyle->SetPadLeftMargin(0.15);
1490 gStyle->SetHistLineWidth(1);
1491 gStyle->SetHistLineColor(kRed);
1492 gStyle->SetFuncWidth(2);
1493 gStyle->SetFuncColor(kGreen);
1494 gStyle->SetLineWidth(2);
1495 gStyle->SetLabelSize(0.045,"xyz");
1496 gStyle->SetLabelOffset(0.01,"y");
1497 gStyle->SetLabelOffset(0.01,"x");
1498 gStyle->SetLabelColor(kBlack,"xyz");
1499 gStyle->SetTitleSize(0.05,"xyz");
1500 gStyle->SetTitleOffset(1.25,"y");
1501 gStyle->SetTitleOffset(1.2,"x");
1502 gStyle->SetTitleFillColor(kWhite);
1503 gStyle->SetTextSizePixels(26);
1504 gStyle->SetTextFont(42);
1505 gStyle->SetLegendBorderSize(0);
1506 gStyle->SetLegendFillColor(kWhite);
1507 gStyle->SetLegendFont(42);
1510 //_____________________________________________________________________________
1511 void AliJetFlowTools::Style(TCanvas* c, TString style)
1513 // set a default style for a canvas
1514 if(!strcmp(style.Data(), "PEARSON")) {
1515 printf(" > style PEARSON canvas < \n");
1516 gStyle->SetOptStat(0);
1521 } else if(!strcmp(style.Data(), "SPECTRUM")) {
1522 printf(" > style SPECTRUM canvas < \n");
1523 gStyle->SetOptStat(0);
1529 } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1531 //_____________________________________________________________________________
1532 void AliJetFlowTools::Style(TVirtualPad* c, TString style, Bool_t legacy)
1534 // set a default style for a canva
1537 c->SetLeftMargin(.25);
1538 c->SetBottomMargin(.25);
1541 if(!strcmp(style.Data(), "PEARSON")) {
1542 printf(" > style PEARSON pad < \n");
1543 gStyle->SetOptStat(0);
1548 } else if(!strcmp(style.Data(), "SPECTRUM")) {
1549 printf(" > style SPECTRUM pad < \n");
1550 gStyle->SetOptStat(0);
1556 } else if (!strcmp(style.Data(), "GRID")) {
1557 printf(" > style GRID pad < \n");
1558 gStyle->SetOptStat(0);
1562 } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1564 //_____________________________________________________________________________
1565 void AliJetFlowTools::Style(TLegend* l)
1567 // set a default style for a legend
1569 l->SetBorderSize(0);
1570 if(gStyle) l->SetTextSize(gStyle->GetTextSize()*.08);
1572 //_____________________________________________________________________________
1573 void AliJetFlowTools::Style(TH1* h, EColor col, histoType type, Bool_t legacy)
1576 h->GetYaxis()->SetNdivisions(505);
1577 h->SetLineColor(col);
1578 h->SetMarkerColor(col);
1580 h->SetMarkerSize(1);
1583 h->GetYaxis()->SetLabelSize(0.05);
1584 h->GetXaxis()->SetLabelSize(0.05);
1585 h->GetYaxis()->SetTitleOffset(1.5);
1586 h->GetXaxis()->SetTitleOffset(1.5);
1587 h->GetYaxis()->SetTitleSize(.05);
1588 h->GetXaxis()->SetTitleSize(.05);
1591 case kInPlaneSpectrum : {
1592 h->SetTitle("IN PLANE");
1593 h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1594 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1596 case kOutPlaneSpectrum : {
1597 h->SetTitle("OUT OF PLANE");
1598 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1599 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1601 case kUnfoldedSpectrum : {
1602 h->SetTitle("UNFOLDED");
1603 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1604 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1606 case kFoldedSpectrum : {
1607 h->SetTitle("FOLDED");
1608 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1609 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1611 case kMeasuredSpectrum : {
1612 h->SetTitle("MEASURED");
1613 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1614 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{cht} (GeV/#it{c})");
1617 h->SetFillColor(col);
1619 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1620 h->SetBarOffset(0.2);
1623 h->SetMarkerStyle(8);
1624 h->SetMarkerSize(1);
1627 h->GetYaxis()->SetTitle("[counts]");
1628 h->GetXaxis()->SetTitle("#Delta #phi");
1633 //_____________________________________________________________________________
1634 void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type, Bool_t legacy)
1637 h->GetYaxis()->SetNdivisions(505);
1638 h->SetLineColor(col);
1639 h->SetMarkerColor(col);
1641 h->SetMarkerSize(1);
1643 h->SetFillColor(kCyan);
1645 h->GetYaxis()->SetLabelSize(0.05);
1646 h->GetXaxis()->SetLabelSize(0.05);
1647 h->GetYaxis()->SetTitleOffset(1.6);
1648 h->GetXaxis()->SetTitleOffset(1.6);
1649 h->GetYaxis()->SetTitleSize(.05);
1650 h->GetXaxis()->SetTitleSize(.05);
1652 h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1654 case kInPlaneSpectrum : {
1655 h->SetTitle("IN PLANE");
1656 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1658 case kOutPlaneSpectrum : {
1659 h->SetTitle("OUT OF PLANE");
1660 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1662 case kUnfoldedSpectrum : {
1663 h->SetTitle("UNFOLDED");
1664 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1666 case kFoldedSpectrum : {
1667 h->SetTitle("FOLDED");
1668 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1670 case kMeasuredSpectrum : {
1671 h->SetTitle("MEASURED");
1672 h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1675 // h->GetYaxis()->SetTitle("#frac{d#it{N_{in plane}^{jet}}}{d#it{p}_{T}} / #frac{d#it{N_{out of plane}^{jet}}}{d#it{p}_{T}}");
1676 h->GetYaxis()->SetTitle("(d#it{N}^{ch, jet}_{in plane}/(d#it{p}_{T}d#eta))/(d#it{N}^{ch,jet}_{out of plane}/(d#it{p}_{T}d#eta))");
1679 // h->GetYaxis()->SetTitle("#it{v}_{2} = #frac{1}{#it{R}} #frac{#pi}{4} #frac{#it{N_{in plane}} - #it{N_{out of plane}}}{#it{N_{in plane}} + #it{N_{out of plane}}}");
1680 h->GetYaxis()->SetTitle("#it{v}_{2}^{ch, jet} \{EP, |#Delta#eta|>0.9 \} ");
1681 h->GetYaxis()->SetRangeUser(-.5, 1.);
1682 h->SetMarkerStyle(8);
1683 h->SetMarkerSize(1);
1688 //_____________________________________________________________________________
1689 void AliJetFlowTools::GetNominalValues(
1690 TH1D*& ratio, // pointer reference, output of this function
1691 TGraphErrors*& v2, // pointer reference, as output of this function
1695 TString outFile) const
1697 // pass clones of the nominal points and nominal v2 values
1698 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close(); // if for some weird reason the unfolding output is still mutable
1699 TFile* readMe(new TFile(inFile.Data(), "READ")); // open unfolding output read-only
1700 if(readMe->IsZombie()) {
1701 printf(" > Fatal error, couldn't read %s for post processing ! < \n", inFile.Data());
1704 printf("\n\n\n\t\t GetNominalValues \n > Recovered the following file structure : \n <");
1706 TFile* output(new TFile(outFile.Data(), "RECREATE")); // create new output file
1707 // get some placeholders, have to be initialized but will be deleted
1708 ratio = new TH1D("nominal", "nominal", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1709 TH1D* nominalIn(new TH1D("nominal in", "nominal in", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1710 TH1D* nominalOut(new TH1D("nominal out", "nominal out", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1711 Int_t iIn[] = {in->At(0), in->At(0)}; // first value is the nominal point
1712 Int_t iOut[] = {out->At(0), out->At(0)};
1714 // call the functions and set the relevant pointer references
1716 DoIntermediateSystematics(
1717 new TArrayI(2, iIn),
1718 new TArrayI(2, iOut),
1719 dud, dud, dud, dud, dud, dud,
1720 ratio, // pointer reference, output of this function
1725 fBinsTrue->At(fBinsTrue->GetSize()-1),
1728 v2 = GetV2(nominalIn, nominalOut, fEventPlaneRes, "nominal v_{2}");
1730 // close the open files, reclaim ownership of histograms which are necessary outside of the file scope
1731 ratio->SetDirectory(0); // disassociate from current gDirectory
1734 //_____________________________________________________________________________
1735 void AliJetFlowTools::GetCorrelatedUncertainty(
1736 TGraphAsymmErrors*& corrRatio, // correlated uncertainty pointer reference
1737 TGraphAsymmErrors*& corrV2, // correlated uncertainty pointer reference
1738 TArrayI* variationsIn, // variantions in plane
1739 TArrayI* variationsOut, // variantions out of plane
1740 Bool_t sym, // treat as symmmetric
1741 TArrayI* variations2ndIn, // second source of variations
1742 TArrayI* variations2ndOut, // second source of variations
1743 Bool_t sym2nd, // treat as symmetric
1744 TString type, // type of uncertaitny
1746 Int_t columns, // divide the output canvasses in this many columns
1747 Float_t rangeLow, // lower pt range
1748 Float_t rangeUp, // upper pt range
1749 Float_t corr, // correlation strength
1750 TString in, // input file name (created by this unfolding class)
1751 TString out // output file name (which will hold results of the systematic test)
1754 // do full systematics
1755 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close(); // if for some weird reason the unfolding output is still mutable
1756 TFile* readMe(new TFile(in.Data(), "READ")); // open unfolding output read-only
1757 if(readMe->IsZombie()) {
1758 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1761 printf("\n\n\n\t\t GetCorrelatedUncertainty \n > Recovered the following file structure : \n <");
1763 TFile* output(new TFile(out.Data(), "RECREATE")); // create new output file
1765 // create some null placeholder pointers
1766 TH1D* relativeErrorVariationInLow(0x0);
1767 TH1D* relativeErrorVariationInUp(0x0);
1768 TH1D* relativeErrorVariationOutLow(0x0);
1769 TH1D* relativeErrorVariationOutUp(0x0);
1770 TH1D* relativeError2ndVariationInLow(0x0);
1771 TH1D* relativeError2ndVariationInUp(0x0);
1772 TH1D* relativeError2ndVariationOutLow(0x0);
1773 TH1D* relativeError2ndVariationOutUp(0x0);
1774 TH1D* relativeStatisticalErrorIn(0x0);
1775 TH1D* relativeStatisticalErrorOut(0x0);
1776 // histo for the nominal ratio in / out
1777 TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1778 TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1779 TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1781 // call the functions
1782 if(variationsIn && variationsOut) {
1783 DoIntermediateSystematics(
1786 relativeErrorVariationInUp, // pointer reference
1787 relativeErrorVariationInLow, // pointer reference
1788 relativeErrorVariationOutUp, // pointer reference
1789 relativeErrorVariationOutLow, // pointer reference
1790 relativeStatisticalErrorIn, // pointer reference
1791 relativeStatisticalErrorOut, // pointer reference
1800 if(relativeErrorVariationInUp) {
1801 // canvas with the error from variation strength
1802 TCanvas* relativeErrorVariation(new TCanvas(Form("relativeError_%s", type.Data()), Form("relativeError_%s", type.Data())));
1803 relativeErrorVariation->Divide(2);
1804 relativeErrorVariation->cd(1);
1805 Style(gPad, "GRID");
1806 relativeErrorVariationInUp->DrawCopy("b");
1807 relativeErrorVariationInLow->DrawCopy("same b");
1808 Style(AddLegend(gPad));
1809 relativeErrorVariation->cd(2);
1810 Style(gPad, "GRID");
1811 relativeErrorVariationOutUp->DrawCopy("b");
1812 relativeErrorVariationOutLow->DrawCopy("same b");
1813 SavePadToPDF(relativeErrorVariation);
1814 Style(AddLegend(gPad));
1815 relativeErrorVariation->Write();
1817 // now smoothen the diced response error (as it is expected to be flat)
1818 // this is done by fitting a constant to the diced resonse error histo
1820 TF1* lin = new TF1("lin", "[0]", rangeLow, rangeUp);
1821 relativeErrorVariationInUp->Fit(lin, "L", "", rangeLow, rangeUp);
1822 if(!gMinuit->fISW[1] == 3) printf(" fit is NOT ok ! " );
1823 for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1824 relativeErrorVariationInUp->SetBinContent(i+1, lin->GetParameter(0));
1826 relativeErrorVariationInLow->Fit(lin, "L", "", rangeLow, rangeUp);
1827 printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
1828 for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1829 relativeErrorVariationInLow->SetBinContent(i+1, 0);//lin->GetParameter(0));
1831 relativeErrorVariationOutUp->Fit(lin, "L", "", rangeLow, rangeUp);
1832 printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
1833 for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1834 relativeErrorVariationOutUp->SetBinContent(i+1, lin->GetParameter(0));
1836 relativeErrorVariationOutLow->Fit(lin, "L", "", rangeLow, rangeUp);
1837 printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
1838 for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1839 relativeErrorVariationOutLow->SetBinContent(i+1, 0);//lin->GetParameter(0));
1843 // call the functions for a second set of systematic sources
1844 if(variations2ndIn && variations2ndOut) {
1845 DoIntermediateSystematics(
1848 relativeError2ndVariationInUp, // pointer reference
1849 relativeError2ndVariationInLow, // pointer reference
1850 relativeError2ndVariationOutUp, // pointer reference
1851 relativeError2ndVariationOutLow, // pointer reference
1852 relativeStatisticalErrorIn, // pointer reference
1853 relativeStatisticalErrorOut, // pointer reference
1862 if(relativeError2ndVariationInUp) {
1863 // canvas with the error from variation strength
1864 TCanvas* relativeError2ndVariation(new TCanvas(Form("relativeError2nd_%s", type2.Data()), Form("relativeError2nd_%s", type2.Data())));
1865 relativeError2ndVariation->Divide(2);
1866 relativeError2ndVariation->cd(1);
1867 Style(gPad, "GRID");
1868 relativeError2ndVariationInUp->DrawCopy("b");
1869 relativeError2ndVariationInLow->DrawCopy("same b");
1870 Style(AddLegend(gPad));
1871 relativeError2ndVariation->cd(2);
1872 Style(gPad, "GRID");
1873 relativeError2ndVariationOutUp->DrawCopy("b");
1874 relativeError2ndVariationOutLow->DrawCopy("same b");
1875 SavePadToPDF(relativeError2ndVariation);
1876 Style(AddLegend(gPad));
1877 relativeError2ndVariation->Write();
1882 // and the placeholder for the final systematic
1883 TH1D* relativeErrorInUp(new TH1D("max correlated uncertainty in plane", "max correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1884 TH1D* relativeErrorOutUp(new TH1D("max correlated uncertainty out of plane", "max correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1885 TH1D* relativeErrorInLow(new TH1D("min correlated uncertainty in plane", "min correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1886 TH1D* relativeErrorOutLow(new TH1D("min correlated uncertainty out of plane", "min correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1887 relativeErrorInUp->SetYTitle("relative uncertainty");
1888 relativeErrorOutUp->SetYTitle("relative uncertainty");
1889 relativeErrorInLow->SetYTitle("relative uncertainty");
1890 relativeErrorOutLow->SetYTitle("relative uncertainty");
1892 // merge the correlated errors (FIXME) trivial for one set
1893 Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.);
1894 Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.);
1895 Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.);
1896 Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.);
1898 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
1899 // for the upper bound
1900 if(relativeErrorVariationInUp) aInUp = relativeErrorVariationInUp->GetBinContent(b+1);
1901 if(relativeErrorVariationOutUp) aOutUp = relativeErrorVariationOutUp->GetBinContent(b+1);
1902 if(relativeError2ndVariationInUp) bInUp = relativeError2ndVariationInUp->GetBinContent(b+1);
1903 if(relativeError2ndVariationOutUp) bInLow = relativeError2ndVariationOutUp->GetBinContent(b+1);
1904 dInUp = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp;
1905 // for a symmetric first variation
1906 if(sym) dInUp += aInLow*aInLow;
1907 if(dInUp > 0) relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(dInUp));
1908 dOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp;
1909 if(sym) dOutUp += aOutLow*aOutLow;
1910 if(dOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(dOutUp));
1911 // for the lower bound
1912 if(relativeErrorVariationInLow) aInLow = relativeErrorVariationInLow->GetBinContent(b+1);
1913 if(relativeErrorVariationOutLow) aOutLow = relativeErrorVariationOutLow->GetBinContent(b+1);
1914 if(relativeError2ndVariationInLow) bInLow = relativeError2ndVariationInLow->GetBinContent(b+1);
1915 if(relativeError2ndVariationOutLow) bOutLow = relativeError2ndVariationOutLow->GetBinContent(b+1);
1916 dInLow = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow;
1917 if(sym) dInLow += aInUp*aInUp;
1918 if(dInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1*TMath::Sqrt(dInLow));
1919 dOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow;
1920 if(sym) dOutLow += aOutUp*aOutUp;
1921 if(dOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(dOutLow));
1923 // project the estimated errors on the nominal ratio
1925 Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
1926 Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
1927 Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
1928 Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
1929 Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
1930 Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
1931 Double_t lowErr(0.), upErr(0.);
1932 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
1933 // add the in and out of plane errors in quadrature
1934 // propagation is tricky because we should consider two cases:
1935 // [1] the error is uncorrelated, and we propagate upper 1 with lower 2 and vice versa
1936 // [2] the error is correlated - in this case upper 1 should be propagated with upper 2
1937 // as the fluctuations are bound to always to of same sign
1939 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(1+i)*relativeErrorOutUp->GetBinContent(i+1);
1940 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
1942 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(1+i)*relativeErrorOutLow->GetBinContent(i+1) -2.*corr*relativeErrorInLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
1943 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1) - 2.*corr*relativeErrorInUp->GetBinContent(i+1)*relativeErrorOutUp->GetBinContent(i+1);
1945 ayl[i] = TMath::Sqrt(TMath::Abs(lowErr))*nominal->GetBinContent(i+1);
1946 ayh[i] = TMath::Sqrt(TMath::Abs(upErr))*nominal->GetBinContent(i+1);
1947 // get the bin width (which is the 'error' on x
1948 Double_t binWidth(nominal->GetBinWidth(i+1));
1949 axl[i] = binWidth/2.;
1950 axh[i] = binWidth/2.;
1951 // now get the coordinate for the point
1952 ax[i] = nominal->GetBinCenter(i+1);
1953 ay[i] = nominal->GetBinContent(i+1);
1955 // save the nominal ratio
1956 TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
1957 corrRatio = (TGraphAsymmErrors*)nominalError->Clone();
1958 nominalError->SetName("correlated uncertainty");
1959 TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
1960 nominalCanvas->Divide(2);
1961 nominalCanvas->cd(1);
1962 Style(nominal, kBlack);
1963 Style(nominalError, kCyan, kRatio);
1964 nominalError->SetLineColor(kCyan);
1965 nominalError->SetMarkerColor(kCyan);
1966 nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1967 nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
1968 nominalError->DrawClone("a2");
1969 nominal->DrawCopy("same E1");
1970 Style(AddLegend(gPad));
1971 Style(gPad, "GRID");
1972 Style(nominalCanvas);
1973 // save nominal v2 and systematic errors
1974 TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
1978 "v_{2} with shape uncertainty",
1982 relativeErrorOutLow,
1984 // pass the nominal values to the pointer references
1985 corrV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
1986 TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
1987 nominalCanvas->cd(2);
1988 Style(nominalV2, kBlack);
1989 Style(nominalV2Error, kCyan, kV2);
1990 nominalV2Error->SetLineColor(kCyan);
1991 nominalV2Error->SetMarkerColor(kCyan);
1992 nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
1993 nominalV2Error->DrawClone("a2");
1994 nominalV2->DrawClone("same E1");
1995 Style(AddLegend(gPad));
1996 Style(gPad, "GRID");
1997 Style(nominalCanvas);
1998 SavePadToPDF(nominalCanvas);
1999 nominalCanvas->Write();
2002 TCanvas* relativeError(new TCanvas("relativeCorrelatedError"," relativeCorrelatedError"));
2003 relativeError->Divide(2);
2004 relativeError->cd(1);
2005 Style(gPad, "GRID");
2006 relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2007 Style(relativeErrorInUp, kBlue, kBar);
2008 Style(relativeErrorInLow, kGreen, kBar);
2009 relativeErrorInUp->DrawCopy("b");
2010 relativeErrorInLow->DrawCopy("same b");
2011 Style(relativeStatisticalErrorIn, kRed);
2012 relativeStatisticalErrorIn->DrawCopy("same");
2013 Style(AddLegend(gPad));
2014 relativeError->cd(2);
2015 Style(gPad, "GRID");
2016 relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2017 Style(relativeErrorOutUp, kBlue, kBar);
2018 Style(relativeErrorOutLow, kGreen, kBar);
2019 relativeErrorOutUp->DrawCopy("b");
2020 relativeErrorOutLow->DrawCopy("same b");
2021 Style(relativeStatisticalErrorOut, kRed);
2022 relativeStatisticalErrorOut->DrawCopy("same");
2023 Style(AddLegend(gPad));
2025 // write the buffered file to disk and close the file
2026 SavePadToPDF(relativeError);
2027 relativeError->Write();
2031 //_____________________________________________________________________________
2032 void AliJetFlowTools::GetShapeUncertainty(
2033 TGraphAsymmErrors*& shapeRatio, // pointer reference to final shape uncertainty of ratio
2034 TGraphAsymmErrors*& shapeV2, // pointer reference to final shape uncertainty of v2
2035 TArrayI* regularizationIn, // regularization strength systematics, in plane
2036 TArrayI* regularizationOut, // regularization strenght systematics, out of plane
2037 TArrayI* trueBinIn, // true bin ranges, in plane
2038 TArrayI* trueBinOut, // true bin ranges, out of plane
2039 TArrayI* recBinIn, // rec bin ranges, in plane
2040 TArrayI* recBinOut, // rec bin ranges, out of plane
2041 TArrayI* methodIn, // method in
2042 TArrayI* methodOut, // method out
2043 Int_t columns, // divide the output canvasses in this many columns
2044 Float_t rangeLow, // lower pt range
2045 Float_t rangeUp, // upper pt range
2046 TString in, // input file name (created by this unfolding class)
2047 TString out // output file name (which will hold results of the systematic test)
2050 // do full systematics
2051 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close(); // if for some weird reason the unfolding output is still mutable
2052 TFile* readMe(new TFile(in.Data(), "READ")); // open unfolding output read-only
2053 if(readMe->IsZombie()) {
2054 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
2057 printf("\n\n\n\t\t DOSYSTEMATICS \n > Recovered the following file structure : \n <");
2059 TFile* output(new TFile(out.Data(), "RECREATE")); // create new output file
2061 // create some null placeholder pointers
2062 TH1D* relativeErrorRegularizationInLow(0x0);
2063 TH1D* relativeErrorRegularizationInUp(0x0);
2064 TH1D* relativeErrorTrueBinInLow(0x0);
2065 TH1D* relativeErrorTrueBinInUp(0x0);
2066 TH1D* relativeErrorRecBinInLow(0x0);
2067 TH1D* relativeErrorRecBinInUp(0x0);
2068 TH1D* relativeErrorMethodInLow(0x0);
2069 TH1D* relativeErrorMethodInUp(0x0);
2070 TH1D* relativeErrorRegularizationOutLow(0x0);
2071 TH1D* relativeErrorRegularizationOutUp(0x0);
2072 TH1D* relativeErrorTrueBinOutLow(0x0);
2073 TH1D* relativeErrorTrueBinOutUp(0x0);
2074 TH1D* relativeErrorRecBinOutLow(0x0);
2075 TH1D* relativeErrorRecBinOutUp(0x0);
2076 TH1D* relativeStatisticalErrorIn(0x0);
2077 TH1D* relativeStatisticalErrorOut(0x0);
2078 TH1D* relativeErrorMethodOutLow(0x0);
2079 TH1D* relativeErrorMethodOutUp(0x0);
2080 // histo for the nominal ratio in / out
2081 TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2082 TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2083 TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2085 // call the functions
2086 if(regularizationIn && regularizationOut) {
2087 DoIntermediateSystematics(
2090 relativeErrorRegularizationInUp, // pointer reference
2091 relativeErrorRegularizationInLow, // pointer reference
2092 relativeErrorRegularizationOutUp, // pointer reference
2093 relativeErrorRegularizationOutLow, // pointer reference
2094 relativeStatisticalErrorIn, // pointer reference
2095 relativeStatisticalErrorOut, // pointer reference
2105 if(relativeErrorRegularizationInUp) {
2106 // canvas with the error from regularization strength
2107 TCanvas* relativeErrorRegularization(new TCanvas("relativeErrorRegularization", "relativeErrorRegularization"));
2108 relativeErrorRegularization->Divide(2);
2109 relativeErrorRegularization->cd(1);
2110 Style(gPad, "GRID");
2111 relativeErrorRegularizationInUp->DrawCopy("b");
2112 relativeErrorRegularizationInLow->DrawCopy("same b");
2113 Style(AddLegend(gPad));
2114 relativeErrorRegularization->cd(2);
2115 Style(gPad, "GRID");
2116 relativeErrorRegularizationOutUp->DrawCopy("b");
2117 relativeErrorRegularizationOutLow->DrawCopy("same b");
2118 SavePadToPDF(relativeErrorRegularization);
2119 Style(AddLegend(gPad));
2120 relativeErrorRegularization->Write();
2123 if(trueBinIn && trueBinOut) {
2124 DoIntermediateSystematics(
2127 relativeErrorTrueBinInUp, // pointer reference
2128 relativeErrorTrueBinInLow, // pointer reference
2129 relativeErrorTrueBinOutUp, // pointer reference
2130 relativeErrorTrueBinOutLow, // pointer reference
2131 relativeStatisticalErrorIn,
2132 relativeStatisticalErrorOut,
2141 if(relativeErrorTrueBinInUp) {
2142 TCanvas* relativeErrorTrueBin(new TCanvas("relativeErrorTrueBin", "relativeErrorTrueBin"));
2143 relativeErrorTrueBin->Divide(2);
2144 relativeErrorTrueBin->cd(1);
2145 Style(gPad, "GRID");
2146 relativeErrorTrueBinInUp->DrawCopy("b");
2147 relativeErrorTrueBinInLow->DrawCopy("same b");
2148 Style(AddLegend(gPad));
2149 relativeErrorTrueBin->cd(2);
2150 Style(gPad, "GRID");
2151 relativeErrorTrueBinOutUp->DrawCopy("b");
2152 relativeErrorTrueBinOutLow->DrawCopy("same b");
2153 SavePadToPDF(relativeErrorTrueBin);
2154 Style(AddLegend(gPad));
2155 relativeErrorTrueBin->Write();
2158 if(recBinIn && recBinOut) {
2159 DoIntermediateSystematics(
2162 relativeErrorRecBinInUp, // pointer reference
2163 relativeErrorRecBinInLow, // pointer reference
2164 relativeErrorRecBinOutUp, // pointer reference
2165 relativeErrorRecBinOutLow, // pointer reference
2166 relativeStatisticalErrorIn,
2167 relativeStatisticalErrorOut,
2177 if(relativeErrorRecBinOutUp) {
2178 // canvas with the error from regularization strength
2179 TCanvas* relativeErrorRecBin(new TCanvas("relativeErrorRecBin"," relativeErrorRecBin"));
2180 relativeErrorRecBin->Divide(2);
2181 relativeErrorRecBin->cd(1);
2182 Style(gPad, "GRID");
2183 relativeErrorRecBinInUp->DrawCopy("b");
2184 relativeErrorRecBinInLow->DrawCopy("same b");
2185 Style(AddLegend(gPad));
2186 relativeErrorRecBin->cd(2);
2187 Style(gPad, "GRID");
2188 relativeErrorRecBinOutUp->DrawCopy("b");
2189 relativeErrorRecBinOutLow->DrawCopy("same b");
2190 SavePadToPDF(relativeErrorRecBin);
2191 Style(AddLegend(gPad));
2192 relativeErrorRecBin->Write();
2195 if(methodIn && methodOut) {
2196 DoIntermediateSystematics(
2199 relativeErrorMethodInUp, // pointer reference
2200 relativeErrorMethodInLow, // pointer reference
2201 relativeErrorMethodOutUp, // pointer reference
2202 relativeErrorMethodOutLow, // pointer reference
2203 relativeStatisticalErrorIn,
2204 relativeStatisticalErrorOut,
2214 if(relativeErrorMethodInUp) {
2215 TCanvas* relativeErrorMethod(new TCanvas("relativeErrorMethod", "relativeErrorMethod"));
2216 relativeErrorMethod->Divide(2);
2217 relativeErrorMethod->cd(1);
2218 Style(gPad, "GRID");
2219 relativeErrorMethodInUp->DrawCopy("b");
2220 relativeErrorMethodInLow->DrawCopy("same b");
2221 Style(AddLegend(gPad));
2222 relativeErrorMethod->cd(2);
2223 Style(gPad, "GRID");
2224 relativeErrorMethodOutUp->DrawCopy("b");
2225 relativeErrorMethodOutLow->DrawCopy("same b");
2226 SavePadToPDF(relativeErrorMethod);
2227 Style(AddLegend(gPad));
2228 relativeErrorMethod->Write();
2232 // and the placeholder for the final systematic
2233 TH1D* relativeErrorInUp(new TH1D("max shape uncertainty in plane", "max shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2234 TH1D* relativeErrorOutUp(new TH1D("max shape uncertainty out of plane", "max shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2235 TH1D* relativeErrorInLow(new TH1D("min shape uncertainty in plane", "min shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2236 TH1D* relativeErrorOutLow(new TH1D("min shape uncertainty out of plane", "min shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2237 relativeErrorInUp->SetYTitle("relative uncertainty");
2238 relativeErrorOutUp->SetYTitle("relative uncertainty");
2239 relativeErrorInLow->SetYTitle("relative uncertainty");
2240 relativeErrorOutLow->SetYTitle("relative uncertainty");
2242 // sum of squares for the total systematic error
2243 Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.), eInUp(0.);
2244 Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.), eOutUp(0.);
2245 Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.), eInLow(0.);
2246 Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.), eOutLow(0.);
2248 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2249 // for the upper bound
2250 if(relativeErrorRegularizationInUp) aInUp = relativeErrorRegularizationInUp->GetBinContent(b+1);
2251 if(relativeErrorRegularizationOutUp) aOutUp = relativeErrorRegularizationOutUp->GetBinContent(b+1);
2252 if(relativeErrorTrueBinInUp) bInUp = relativeErrorTrueBinInUp->GetBinContent(b+1);
2253 if(relativeErrorTrueBinOutUp) bOutUp = relativeErrorTrueBinOutUp->GetBinContent(b+1);
2254 if(relativeErrorRecBinInUp) cInUp = relativeErrorRecBinInUp->GetBinContent(b+1);
2255 if(relativeErrorRecBinOutUp) cOutUp = relativeErrorRecBinOutUp->GetBinContent(b+1);
2256 if(relativeErrorMethodInUp) dInUp = relativeErrorMethodInUp->GetBinContent(b+1);
2257 if(relativeErrorMethodOutUp) dOutUp = relativeErrorMethodOutUp->GetBinContent(b+1);
2258 eInUp = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp + dInUp*dInUp;
2259 if(eInUp > 0) relativeErrorInUp->SetBinContent(b+1, 1.*TMath::Sqrt(eInUp));
2260 eOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp + dOutUp*dOutUp;
2261 if(eOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, 1.*TMath::Sqrt(eOutUp));
2262 // for the lower bound
2263 if(relativeErrorRegularizationInLow) aInLow = relativeErrorRegularizationInLow->GetBinContent(b+1);
2264 if(relativeErrorRegularizationOutLow) aOutLow = relativeErrorRegularizationOutLow->GetBinContent(b+1);
2265 if(relativeErrorTrueBinInLow) bInLow = relativeErrorTrueBinInLow->GetBinContent(b+1);
2266 if(relativeErrorTrueBinOutLow) bOutLow = relativeErrorTrueBinOutLow->GetBinContent(b+1);
2267 if(relativeErrorRecBinInLow) cInLow = relativeErrorRecBinInLow->GetBinContent(b+1);
2268 if(relativeErrorRecBinOutLow) cOutLow = relativeErrorRecBinOutLow->GetBinContent(b+1);
2269 if(relativeErrorMethodInLow) dInLow = relativeErrorMethodInLow->GetBinContent(b+1);
2270 if(relativeErrorMethodOutLow) dOutLow = relativeErrorMethodOutLow->GetBinContent(b+1);
2271 if(fSymmRMS) { // take first category as symmetric
2275 cOutLow = cOutUp; // other sources
2276 if(dInLow < dInUp) dInLow = dInUp;
2277 if(dOutLow < dOutUp) dOutLow = dOutUp;
2280 eInLow = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow + dInLow*dInLow;
2281 if(eInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(eInLow));
2282 eOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow + dOutLow*dOutLow;
2283 if(eOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(eOutLow));
2285 // project the estimated errors on the nominal ratio
2287 Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
2288 Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
2289 Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
2290 Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
2291 Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
2292 Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
2293 Double_t lowErr(0.), upErr(0.);
2294 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
2295 // add the in and out of plane errors in quadrature
2296 // take special care here: to propagate the assymetric error, we need to correlate the
2297 // InLow with OutUp (minimum value of ratio) and InUp with OutLow (maximum value of ratio)
2298 lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(1+i)*relativeErrorOutUp->GetBinContent(i+1);
2299 upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
2301 ayl[i] = TMath::Sqrt(lowErr)*nominal->GetBinContent(i+1);
2302 ayh[i] = TMath::Sqrt(upErr)*nominal->GetBinContent(i+1);
2303 // get the bin width (which is the 'error' on x
2304 Double_t binWidth(nominal->GetBinWidth(i+1));
2305 axl[i] = binWidth/2.;
2306 axh[i] = binWidth/2.;
2307 // now get the coordinate for the point
2308 ax[i] = nominal->GetBinCenter(i+1);
2309 ay[i] = nominal->GetBinContent(i+1);
2311 // save the nominal ratio
2312 TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
2313 shapeRatio = (TGraphAsymmErrors*)nominalError->Clone();
2314 nominalError->SetName("shape uncertainty");
2315 TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
2316 nominalCanvas->Divide(2);
2317 nominalCanvas->cd(1);
2318 Style(nominal, kBlack);
2319 Style(nominalError, kCyan, kRatio);
2320 nominalError->SetLineColor(kCyan);
2321 nominalError->SetMarkerColor(kCyan);
2322 nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2323 nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
2324 nominalError->DrawClone("a2");
2325 nominal->DrawCopy("same E1");
2326 Style(AddLegend(gPad));
2327 Style(gPad, "GRID");
2328 Style(nominalCanvas);
2329 // save nominal v2 and systematic errors
2330 TGraphAsymmErrors* nominalV2Error(GetV2WithSystematicErrors(
2334 "v_{2} with shape uncertainty",
2338 relativeErrorOutLow));
2339 shapeV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
2340 TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
2341 nominalCanvas->cd(2);
2342 Style(nominalV2, kBlack);
2343 Style(nominalV2Error, kCyan, kV2);
2344 nominalV2Error->SetLineColor(kCyan);
2345 nominalV2Error->SetMarkerColor(kCyan);
2346 nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2347 nominalV2Error->DrawClone("a2");
2348 nominalV2->DrawClone("same E1");
2349 Style(AddLegend(gPad));
2350 Style(gPad, "GRID");
2351 Style(nominalCanvas);
2352 SavePadToPDF(nominalCanvas);
2353 nominalCanvas->Write();
2356 TCanvas* relativeError(new TCanvas("relativeError"," relativeError"));
2357 relativeError->Divide(2);
2358 relativeError->cd(1);
2359 Style(gPad, "GRID");
2360 relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2361 Style(relativeErrorInUp, kBlue, kBar);
2362 Style(relativeErrorInLow, kGreen, kBar);
2363 relativeErrorInUp->DrawCopy("b");
2364 relativeErrorInLow->DrawCopy("same b");
2365 Style(relativeStatisticalErrorIn, kRed);
2366 relativeStatisticalErrorIn->DrawCopy("same");
2367 Style(AddLegend(gPad));
2368 relativeError->cd(2);
2369 Style(gPad, "GRID");
2370 relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2371 Style(relativeErrorOutUp, kBlue, kBar);
2372 Style(relativeErrorOutLow, kGreen, kBar);
2373 relativeErrorOutUp->DrawCopy("b");
2374 relativeErrorOutLow->DrawCopy("same b");
2375 Style(relativeStatisticalErrorOut, kRed);
2376 relativeStatisticalErrorOut->DrawCopy("same");
2377 Style(AddLegend(gPad));
2379 // write the buffered file to disk and close the file
2380 SavePadToPDF(relativeError);
2381 relativeError->Write();
2385 //_____________________________________________________________________________
2386 void AliJetFlowTools::DoIntermediateSystematics(
2387 TArrayI* variationsIn, // variantions in plane
2388 TArrayI* variationsOut, // variantions out of plane
2389 TH1D*& relativeErrorInUp, // pointer reference to minimum relative error histogram in plane
2390 TH1D*& relativeErrorInLow, // pointer reference to maximum relative error histogram in plane
2391 TH1D*& relativeErrorOutUp, // pointer reference to minimum relative error histogram out of plane
2392 TH1D*& relativeErrorOutLow, // pointer reference to maximum relative error histogram out of plane
2393 TH1D*& relativeStatisticalErrorIn, // relative systematic error on ratio
2394 TH1D*& relativeStatisticalErrorOut, // relative systematic error on ratio
2395 TH1D*& nominal, // clone of the nominal ratio in / out of plane
2396 TH1D*& nominalIn, // clone of the nominal in plane yield
2397 TH1D*& nominalOut, // clone of the nominal out of plane yield
2398 Int_t columns, // divide the output canvasses in this many columns
2399 Float_t rangeLow, // lower pt range
2400 Float_t rangeUp, // upper pt range
2401 TFile* readMe, // input file name (created by this unfolding class)
2402 TString source, // source of the variation
2403 Bool_t RMS // return RMS of distribution of variations as error
2406 // intermediate systematic check function. first index of supplied array is nominal value
2407 TList* listOfKeys((TList*)readMe->GetListOfKeys());
2409 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2412 // check input params
2413 if(variationsIn->GetSize() != variationsOut->GetSize()) {
2414 printf(" > DoSystematics: fatal error, input arrays have different sizes ! < \n ");
2417 TDirectoryFile* defRootDirIn(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(0))->GetName())));
2418 TDirectoryFile* defRootDirOut(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(0))->GetName())));
2419 if(!(defRootDirIn && defRootDirOut)) {
2420 printf(" > DoSystematics: fatal error, couldn't retrieve nominal values ! < \n ");
2423 TString defIn(defRootDirIn->GetName());
2424 TString defOut(defRootDirOut->GetName());
2426 // define lines to make the output prettier
2427 TLine* lineLow(new TLine(rangeLow, 0., rangeLow, 2.));
2428 TLine* lineUp(new TLine(rangeUp, 0., rangeUp, 2.));
2429 lineLow->SetLineColor(11);
2430 lineUp->SetLineColor(11);
2431 lineLow->SetLineWidth(3);
2432 lineUp->SetLineWidth(3);
2434 // define an output histogram with the maximum relative error from this systematic constribution
2435 // if the option RMS is set to false, sigma is not really a standard deviation but holds the maximum (or minimum) relative value that the data has
2436 // reached in this function call.
2437 // if the option RMS is set to true, sigma holds the RMS value (equal to sigma as the mean is zero for relative errors) of the distribution of variations
2438 // which should correspond to a 68% confidence level
2439 relativeErrorInUp = new TH1D(Form("relative error (up) from %s", source.Data()), Form("relative error (up) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2440 relativeErrorInLow = new TH1D(Form("relative error (low) from %s", source.Data()), Form("relative error (low) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2441 relativeErrorOutUp = new TH1D(Form("relative error (up) from %s", source.Data()), Form("relative error (up) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2442 relativeErrorOutLow = new TH1D(Form("relative error (low) from %s", source.Data()), Form("relative error (low) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2443 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2445 relativeErrorInUp->SetBinContent(b+1, 1.);
2446 relativeErrorInUp->SetBinError(b+1, 0.);
2447 relativeErrorOutUp->SetBinContent(b+1, 1.);
2448 relativeErrorOutUp->SetBinError(b+1, .0);
2449 relativeErrorInLow->SetBinContent(b+1, 1.);
2450 relativeErrorInLow->SetBinError(b+1, 0.);
2451 relativeErrorOutLow->SetBinContent(b+1, 1.);
2452 relativeErrorOutLow->SetBinError(b+1, .0);
2454 relativeErrorInUp->SetBinContent(b+1, 0.);
2455 relativeErrorInUp->SetBinError(b+1, 0.);
2456 relativeErrorOutUp->SetBinContent(b+1, 0.);
2457 relativeErrorOutUp->SetBinError(b+1, 0.);
2458 relativeErrorInLow->SetBinContent(b+1, 0.);
2459 relativeErrorInLow->SetBinError(b+1, 0.);
2460 relativeErrorOutLow->SetBinContent(b+1, 0.);
2461 relativeErrorOutLow->SetBinError(b+1, 0.);
2464 Int_t relativeErrorInUpN[100] = {0};
2465 Int_t relativeErrorOutUpN[100] = {0};
2466 Int_t relativeErrorInLowN[100] = {0};
2467 Int_t relativeErrorOutLowN[100] = {0};
2469 // define an output histogram with the systematic error from this systematic constribution
2470 if(!relativeStatisticalErrorIn && !relativeStatisticalErrorOut) {
2471 relativeStatisticalErrorIn = new TH1D("relative statistical error, in plane", "relative statistical error, in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2472 relativeStatisticalErrorOut = new TH1D("relative statistical error, out of plane", "relative statistital error, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2475 // prepare necessary canvasses
2476 TCanvas* canvasIn(new TCanvas(Form("SYST_%s_PearsonIn", source.Data()), Form("SYST_%s_PearsonIn", source.Data())));
2477 TCanvas* canvasOut(0x0);
2478 if(fDphiUnfolding) canvasOut = new TCanvas(Form("SYST_%s_PearsonOut", source.Data()), Form("SYST_%s_PearsonOut", source.Data()));
2479 TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas(Form("SYST_%s_RefoldedIn", source.Data()), Form("SYST_%s_RefoldedIn", source.Data())));
2480 TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
2481 if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas(Form("SYST_%s_RefoldedOut", source.Data()), Form("SYST_%s_RefoldedOut", source.Data()));
2482 TCanvas* canvasSpectraIn(new TCanvas(Form("SYST_%s_SpectraIn", source.Data()), Form("SYST_%s_SpectraIn", source.Data())));
2483 TCanvas* canvasSpectraOut(0x0);
2484 if(fDphiUnfolding) canvasSpectraOut = new TCanvas(Form("SYST_%s_SpectraOut", source.Data()), Form("SYST_%s_SpectraOut", source.Data()));
2485 TCanvas* canvasRatio(0x0);
2486 if(fDphiUnfolding) canvasRatio = new TCanvas(Form("SYST_%s_Ratio", source.Data()), Form("SYST_%s_Ratio", source.Data()));
2487 TCanvas* canvasV2(0x0);
2488 if(fDphiUnfolding) canvasV2 = new TCanvas(Form("SYST_%s_V2", source.Data()), Form("SYST_%s_V2", source.Data()));
2489 TCanvas* canvasMISC(new TCanvas(Form("SYST_%s_MISC", source.Data()), Form("SYST_%s_MISC", source.Data())));
2490 TCanvas* canvasMasterIn(new TCanvas(Form("SYST_%s_defaultIn", source.Data()), Form("SYST_%s_defaultIn", source.Data())));
2491 TCanvas* canvasMasterOut(0x0);
2492 if(fDphiUnfolding) canvasMasterOut = new TCanvas(Form("SYST_%s_defaultOut", source.Data()), Form("SYST_%s_defaultOut", source.Data()));
2493 (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
2495 TCanvas* canvasProfiles(new TCanvas(Form("SYST_%s_canvasProfiles", source.Data()), Form("SYST_%s_canvasProfiles", source.Data())));
2496 canvasProfiles->Divide(2);
2497 TProfile* ratioProfile(new TProfile(Form("SYST_%s_ratioProfile", source.Data()), Form("SYST_%s_ratioProfile", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2498 TProfile* v2Profile(new TProfile(Form("SYST_%s_v2Profile", source.Data()), Form("SYST_%s_v2Profile", source.Data()),fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2499 // get an estimate of the number of outputs and find the default set
2502 columns = variationsIn->GetSize()-1;
2503 (TMath::Floor(variationsIn->GetSize()/(float)columns)+((variationsIn->GetSize()%columns)>0));
2504 canvasIn->Divide(columns, rows);
2505 if(canvasOut) canvasOut->Divide(columns, rows);
2506 canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
2507 if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
2508 canvasSpectraIn->Divide(columns, rows);
2509 if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
2510 if(canvasRatio) canvasRatio->Divide(columns, rows);
2511 if(canvasV2) canvasV2->Divide(columns, rows);
2512 canvasMasterIn->Divide(columns, rows);
2513 if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
2515 // prepare a separate set of canvases to hold the nominal points
2516 TCanvas* canvasNominalIn(new TCanvas(Form("Nominal_%s_PearsonIn", source.Data()), Form("Nominal_%s_PearsonIn", source.Data())));
2517 TCanvas* canvasNominalOut(0x0);
2518 if(fDphiUnfolding) canvasNominalOut = new TCanvas(Form("Nominal_%s_PearsonOut", source.Data()), Form("Nominal_%s_PearsonOut", source.Data()));
2519 TCanvas* canvasNominalRatioMeasuredRefoldedIn(new TCanvas(Form("Nominal_%s_RefoldedIn", source.Data()), Form("Nominal_%s_RefoldedIn", source.Data())));
2520 TCanvas* canvasNominalRatioMeasuredRefoldedOut(0x0);
2521 if(fDphiUnfolding) canvasNominalRatioMeasuredRefoldedOut = new TCanvas(Form("Nominal_%s_RefoldedOut", source.Data()), Form("Nominal_%s_RefoldedOut", source.Data()));
2522 TCanvas* canvasNominalSpectraIn(new TCanvas(Form("Nominal_%s_SpectraIn", source.Data()), Form("Nominal_%s_SpectraIn", source.Data())));
2523 TCanvas* canvasNominalSpectraOut(0x0);
2524 if(fDphiUnfolding) canvasNominalSpectraOut = new TCanvas(Form("Nominal_%s_SpectraOut", source.Data()), Form("Nominal_%s_SpectraOut", source.Data()));
2525 TCanvas* canvasNominalRatio(0x0);
2526 if(fDphiUnfolding) canvasNominalRatio = new TCanvas(Form("Nominal_%s_Ratio", source.Data()), Form("Nominal_%s_Ratio", source.Data()));
2527 TCanvas* canvasNominalV2(0x0);
2528 if(fDphiUnfolding) canvasNominalV2 = new TCanvas(Form("Nominal_%s_V2", source.Data()), Form("Nominal_%s_V2", source.Data()));
2529 TCanvas* canvasNominalMISC(new TCanvas(Form("Nominal_%s_MISC", source.Data()), Form("Nominal_%s_MISC", source.Data())));
2530 TCanvas* canvasNominalMasterIn(new TCanvas(Form("Nominal_%s_defaultIn", source.Data()), Form("Nominal_%s_defaultIn", source.Data())));
2531 TCanvas* canvasNominalMasterOut(0x0);
2532 if(fDphiUnfolding) canvasNominalMasterOut = new TCanvas(Form("Nominal_%s_defaultOut", source.Data()), Form("Nominal_%s_defaultOut", source.Data()));
2533 (fDphiUnfolding) ? canvasNominalMISC->Divide(4, 2) : canvasNominalMISC->Divide(4, 1);
2535 canvasNominalSpectraIn->Divide(2);
2536 if(canvasNominalSpectraOut) canvasNominalSpectraOut->Divide(2);
2538 canvasNominalMasterIn->Divide(2);
2539 if(canvasNominalMasterOut) canvasNominalMasterOut->Divide(2);
2541 // extract the default output
2542 TH1D* defaultUnfoldedJetSpectrumIn(0x0);
2543 TH1D* defaultUnfoldedJetSpectrumOut(0x0);
2544 TDirectoryFile* defDirIn = (TDirectoryFile*)defRootDirIn->Get(Form("InPlane___%s", defIn.Data()));
2545 TDirectoryFile* defDirOut = (TDirectoryFile*)defRootDirOut->Get(Form("OutOfPlane___%s", defOut.Data()));
2546 if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", defIn.Data()));
2547 if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", defOut.Data()));
2548 printf(" > succesfully extracted default results < \n");
2550 // loop through the directories, only plot the graphs if the deconvolution converged
2551 TDirectoryFile* tempDirIn(0x0);
2552 TDirectoryFile* tempDirOut(0x0);
2553 TDirectoryFile* tempIn(0x0);
2554 TDirectoryFile* tempOut(0x0);
2555 TH1D* unfoldedSpectrumInForRatio(0x0);
2556 TH1D* unfoldedSpectrumOutForRatio(0x0);
2557 for(Int_t i(0), j(-1); i < variationsIn->GetSize(); i++) {
2558 tempDirIn = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(i))->GetName())));
2559 tempDirOut = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(i))->GetName())));
2560 if(!(tempDirIn && tempDirOut)) {
2561 printf(" > DoSystematics: couldn't get a set of variations < \n");
2564 TString dirNameIn(tempDirIn->GetName());
2565 TString dirNameOut(tempDirOut->GetName());
2566 // try to read the in- and out of plane subdirs
2567 tempIn = (TDirectoryFile*)tempDirIn->Get(Form("InPlane___%s", dirNameIn.Data()));
2568 tempOut = (TDirectoryFile*)tempDirOut->Get(Form("OutOfPlane___%s", dirNameOut.Data()));
2571 // to see if the unfolding converged try to extract the pearson coefficients
2572 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirNameIn.Data())));
2574 printf(" - %s in plane converged \n", dirNameIn.Data());
2576 if(i==0) canvasNominalIn->cd(j);
2577 Style(gPad, "PEARSON");
2578 pIn->DrawCopy("colz");
2579 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirNameIn.Data())));
2582 printf(" > found RatioRefoldedMeasured < \n");
2583 canvasRatioMeasuredRefoldedIn->cd(j);
2584 if(i==0) canvasNominalRatioMeasuredRefoldedIn->cd(j);
2585 Style(gPad, "GRID");
2586 rIn->SetFillColor(kRed);
2589 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2590 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2591 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirNameIn.Data())));
2592 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirNameIn.Data())));
2593 if(dvector && avalue && rm && eff) {
2599 if(i==0) canvasNominalMISC->cd(1);
2600 Style(gPad, "SPECTRUM");
2601 dvector->DrawCopy();
2603 if(i==0) canvasNominalMISC->cd(2);
2604 Style(gPad, "SPECTRUM");
2607 if(i==0) canvasNominalMISC->cd(3);
2608 Style(gPad, "PEARSON");
2609 rm->DrawCopy("colz");
2611 if(i==0) canvasNominalMISC->cd(4);
2612 Style(gPad, "GRID");
2614 } else if(rm && eff) {
2618 if(i==0) canvasNominalMISC->cd(3);
2619 Style(gPad, "PEARSON");
2620 rm->DrawCopy("colz");
2622 if(i==0) canvasNominalMISC->cd(4);
2623 Style(gPad, "GRID");
2627 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirNameIn.Data())));
2628 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirNameIn.Data())));
2629 unfoldedSpectrumInForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2630 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirNameIn.Data())));
2631 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2632 if(defaultUnfoldedJetSpectrumIn) {
2633 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2634 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirNameIn.Data())));
2635 temp->Divide(unfoldedSpectrum);
2636 // get the absolute relative error
2637 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2638 if(!RMS) { // save the maximum deviation that a variation can cause
2639 // the variation is HIGHER than the nominal point, so the bar goes UP
2640 if( temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorInUp->GetBinContent(b+1)) {
2641 relativeErrorInUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2642 relativeErrorInUp->SetBinError(b+1, 0.);
2644 // the variation is LOWER than the nominal point, so the bar goes DOWN
2645 else if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorInLow->GetBinContent(b+1)) {
2646 relativeErrorInLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2647 relativeErrorInLow->SetBinError(b+1, 0.);
2649 } else if (RMS && !fSymmRMS) { // save info necessary for evaluating the RMS of a distribution of variations
2650 printf(" oops shouldnt be here \n " );
2651 if(temp->GetBinContent(b+1) < 1) {
2652 relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2653 relativeErrorInUpN[b]++;
2655 // the variation is LOWER than the nominal point, so the bar goes DOWN
2656 else if(temp->GetBinContent(b+1) > 1) {
2657 relativeErrorInLow->SetBinContent(b+1, relativeErrorInLow->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2658 relativeErrorInLowN[b]++;
2660 } else if (fSymmRMS) {
2661 // save symmetric sum of square to get a symmetric rms
2662 relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(temp->GetBinContent(b+1)-1., 2));
2663 relativeErrorInUpN[b]++;
2665 if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorIn->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2667 temp->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2668 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2669 temp->GetYaxis()->SetTitle("ratio");
2670 canvasMasterIn->cd(j);
2671 temp->GetYaxis()->SetRangeUser(0., 2);
2672 Style(gPad, "GRID");
2674 canvasNominalMasterIn->cd(1);
2675 Style(gPad, "GRID");
2677 TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2678 tempSyst->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2679 Style(tempSyst, (EColor)(i+2));
2680 if(i==1) tempSyst->DrawCopy();
2681 else tempSyst->DrawCopy("same");
2684 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirNameIn.Data())));
2685 canvasSpectraIn->cd(j);
2686 if(i==0) canvasNominalSpectraIn->cd(1);
2688 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2689 unfoldedSpectrum->DrawCopy();
2690 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2691 inputSpectrum->DrawCopy("same");
2692 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2693 refoldedSpectrum->DrawCopy("same");
2694 TLegend* l(AddLegend(gPad));
2696 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2697 Float_t chi(fitStatus->GetBinContent(1));
2698 Float_t pen(fitStatus->GetBinContent(2));
2699 Int_t dof((int)fitStatus->GetBinContent(3));
2700 Float_t beta(fitStatus->GetBinContent(4));
2701 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2702 } else if (fitStatus) { // only available in SVD fit
2703 Int_t reg((int)fitStatus->GetBinContent(1));
2704 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2706 canvasNominalSpectraIn->cd(2);
2707 TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2708 tempSyst->SetTitle(Form("[%s]", dirNameIn.Data()));
2709 Style(tempSyst, (EColor)(i+2));
2710 Style(gPad, "SPECTRUM");
2711 if(i==0) tempSyst->DrawCopy();
2712 else tempSyst->DrawCopy("same");
2716 TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirNameOut.Data())));
2718 printf(" - %s out of plane converged \n", dirNameOut.Data());
2720 if(i==0) canvasNominalOut->cd(j);
2721 Style(gPad, "PEARSON");
2722 pOut->DrawCopy("colz");
2723 TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirNameOut.Data())));
2726 printf(" > found RatioRefoldedMeasured < \n");
2727 canvasRatioMeasuredRefoldedOut->cd(j);
2728 if(i==0) canvasNominalRatioMeasuredRefoldedOut->cd(j);
2729 Style(gPad, "GRID");
2730 rOut->SetFillColor(kRed);
2733 TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2734 TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2735 TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirNameOut.Data())));
2736 TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirNameOut.Data())));
2737 if(dvector && avalue && rm && eff) {
2743 if(i==0) canvasNominalMISC->cd(5);
2744 Style(gPad, "SPECTRUM");
2745 dvector->DrawCopy();
2747 if(i==0) canvasNominalMISC->cd(6);
2748 Style(gPad, "SPECTRUM");
2751 if(i==0) canvasNominalMISC->cd(7);
2752 Style(gPad, "PEARSON");
2753 rm->DrawCopy("colz");
2755 if(i==0) canvasNominalMISC->cd(8);
2756 Style(gPad, "GRID");
2758 } else if(rm && eff) {
2762 if(i==0) canvasNominalMISC->cd(7);
2763 Style(gPad, "PEARSON");
2764 rm->DrawCopy("colz");
2766 if(i==0) canvasNominalMISC->cd(8);
2767 Style(gPad, "GRID");
2771 TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirNameOut.Data())));
2772 TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirNameOut.Data())));
2773 unfoldedSpectrumOutForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2774 TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirNameOut.Data())));
2775 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2776 if(defaultUnfoldedJetSpectrumOut) {
2777 Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
2778 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirNameOut.Data())));
2779 temp->Divide(unfoldedSpectrum);
2780 // get the absolute relative error
2781 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2783 // check if the error is larger than the current maximum
2784 if(temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorOutUp->GetBinContent(b+1)) {
2785 relativeErrorOutUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2786 relativeErrorOutUp->SetBinError(b+1, 0.);
2788 // check if the error is smaller than the current minimum
2789 else if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorOutLow->GetBinContent(b+1)) {
2790 relativeErrorOutLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2791 relativeErrorOutLow->SetBinError(b+1, 0.);
2793 } else if (RMS && !fSymmRMS) {
2794 printf(" OOps \n ");
2795 if(temp->GetBinContent(b+1) < 1) {
2796 relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2797 relativeErrorOutUpN[b]++;
2799 else if(temp->GetBinContent(b+1) > 1) {
2800 relativeErrorOutLow->SetBinContent(b+1, relativeErrorOutLow->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2801 relativeErrorOutLowN[b]++;
2803 } else if (fSymmRMS) {
2804 // save symmetric rms value
2805 relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(temp->GetBinContent(b+1)-1., 2));
2806 relativeErrorOutUpN[b]++;
2808 if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorOut->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2810 temp->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2811 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2812 temp->GetYaxis()->SetTitle("ratio");
2813 canvasMasterOut->cd(j);
2814 temp->GetYaxis()->SetRangeUser(0., 2);
2815 Style(gPad, "GRID");
2817 canvasNominalMasterOut->cd(1);
2818 Style(gPad, "GRID");
2820 TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2821 tempSyst->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
2822 Style(tempSyst, (EColor)(i+2));
2823 if(i==1) tempSyst->DrawCopy();
2824 else tempSyst->DrawCopy("same");
2827 TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirNameOut.Data())));
2828 canvasSpectraOut->cd(j);
2829 if(i==0) canvasNominalSpectraOut->cd(1);
2831 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2832 unfoldedSpectrum->DrawCopy();
2833 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2834 inputSpectrum->DrawCopy("same");
2835 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2836 refoldedSpectrum->DrawCopy("same");
2837 TLegend* l(AddLegend(gPad));
2839 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2840 Float_t chi(fitStatus->GetBinContent(1));
2841 Float_t pen(fitStatus->GetBinContent(2));
2842 Int_t dof((int)fitStatus->GetBinContent(3));
2843 Float_t beta(fitStatus->GetBinContent(4));
2844 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2845 } else if (fitStatus) { // only available in SVD fit
2846 Int_t reg((int)fitStatus->GetBinContent(1));
2847 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2849 canvasNominalSpectraOut->cd(2);
2850 TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2851 tempSyst->SetTitle(Form("[%s]", dirNameOut.Data()));
2852 Style(tempSyst, (EColor)(i+2));
2853 Style(gPad, "SPECTRUM");
2854 if(i==0) tempSyst->DrawCopy();
2855 else tempSyst->DrawCopy("same");
2858 if(canvasRatio && canvasV2) {
2861 canvasNominalRatio->cd(j);
2862 if(nominal && nominalIn && nominalOut) {
2863 // if a nominal ratio is requested, delete the placeholder and update the nominal point
2867 nominalIn = (TH1D*)unfoldedSpectrumInForRatio->Clone("in plane jet yield");
2868 nominalOut = (TH1D*)unfoldedSpectrumOutForRatio->Clone("out of plane jet yield");
2869 nominal = (TH1D*)unfoldedSpectrumInForRatio->Clone("ratio in plane out of plane");
2870 nominal->Divide(nominal, unfoldedSpectrumOutForRatio); // standard root divide for uncorrelated histos
2873 TGraphErrors* ratioYield(GetRatio(unfoldedSpectrumInForRatio, unfoldedSpectrumOutForRatio, TString(Form("ratio [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2874 Double_t _tempx(0), _tempy(0);
2877 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2878 ratioYield->GetPoint(b, _tempx, _tempy);
2879 ratioProfile->Fill(_tempx, _tempy);
2881 ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
2882 ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2883 ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
2884 ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2885 ratioYield->SetFillColor(kRed);
2886 ratioYield->Draw("ap");
2889 if(i==0) canvasNominalV2->cd(j);
2890 TGraphErrors* ratioV2(GetV2(unfoldedSpectrumInForRatio,unfoldedSpectrumOutForRatio, fEventPlaneRes, TString(Form("v_{2} [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
2893 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
2894 ratioV2->GetPoint(b, _tempx, _tempy);
2895 v2Profile->Fill(_tempx, _tempy);
2898 v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
2899 v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2900 ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
2901 ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2902 ratioV2->SetFillColor(kRed);
2903 ratioV2->Draw("ap");
2906 delete unfoldedSpectrumInForRatio;
2907 delete unfoldedSpectrumOutForRatio;
2909 // save the canvasses
2910 canvasProfiles->cd(1);
2911 Style(ratioProfile);
2912 ratioProfile->DrawCopy();
2913 canvasProfiles->cd(2);
2915 v2Profile->DrawCopy();
2916 SavePadToPDF(canvasProfiles);
2917 canvasProfiles->Write();
2918 SavePadToPDF(canvasIn);
2921 SavePadToPDF(canvasOut);
2924 SavePadToPDF(canvasRatioMeasuredRefoldedIn);
2925 canvasRatioMeasuredRefoldedIn->Write();
2926 if(canvasRatioMeasuredRefoldedOut) {
2927 SavePadToPDF(canvasRatioMeasuredRefoldedOut);
2928 canvasRatioMeasuredRefoldedOut->Write();
2930 SavePadToPDF(canvasSpectraIn);
2931 canvasSpectraIn->Write();
2932 if(canvasSpectraOut) {
2933 SavePadToPDF(canvasSpectraOut);
2934 canvasSpectraOut->Write();
2935 SavePadToPDF(canvasRatio);
2936 canvasRatio->Write();
2937 SavePadToPDF(canvasV2);
2940 SavePadToPDF(canvasMasterIn);
2941 canvasMasterIn->Write();
2942 if(canvasMasterOut) {
2943 SavePadToPDF(canvasMasterOut);
2944 canvasMasterOut->Write();
2946 SavePadToPDF(canvasMISC);
2947 canvasMISC->Write();
2948 // save the nomial canvasses
2949 SavePadToPDF(canvasNominalIn);
2950 canvasNominalIn->Write();
2951 if(canvasNominalOut) {
2952 SavePadToPDF(canvasNominalOut);
2953 canvasNominalOut->Write();
2955 SavePadToPDF(canvasNominalRatioMeasuredRefoldedIn);
2956 canvasNominalRatioMeasuredRefoldedIn->Write();
2957 if(canvasNominalRatioMeasuredRefoldedOut) {
2958 SavePadToPDF(canvasNominalRatioMeasuredRefoldedOut);
2959 canvasNominalRatioMeasuredRefoldedOut->Write();
2961 canvasNominalSpectraIn->cd(2);
2962 Style(AddLegend(gPad));
2963 SavePadToPDF(canvasNominalSpectraIn);
2964 canvasNominalSpectraIn->Write();
2965 if(canvasNominalSpectraOut) {
2966 canvasNominalSpectraOut->cd(2);
2967 Style(AddLegend(gPad));
2968 SavePadToPDF(canvasNominalSpectraOut);
2969 canvasNominalSpectraOut->Write();
2970 SavePadToPDF(canvasNominalRatio);
2971 canvasNominalRatio->Write();
2972 SavePadToPDF(canvasNominalV2);
2973 canvasNominalV2->Write();
2975 canvasNominalMasterIn->cd(1);
2976 Style(AddLegend(gPad));
2977 lineUp->DrawClone("same");
2978 lineLow->DrawClone("same");
2979 SavePadToPDF(canvasNominalMasterIn);
2980 canvasNominalMasterIn->Write();
2981 if(canvasNominalMasterOut) {
2982 canvasNominalMasterOut->cd(1);
2983 Style(AddLegend(gPad));
2984 lineUp->DrawClone("same");
2985 lineLow->DrawClone("same");
2986 SavePadToPDF(canvasNominalMasterOut);
2987 canvasNominalMasterOut->Write();
2989 SavePadToPDF(canvasNominalMISC);
2990 canvasNominalMISC->Write();
2992 // save the relative errors
2993 for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2994 // to arrive at a min and max from here, combine in up and out low
2996 relativeErrorInUp->SetBinContent(b+1, -1.*(relativeErrorInUp->GetBinContent(b+1)-1));
2997 relativeErrorInUp->SetBinError(b+1, 0.);
2998 relativeErrorOutUp->SetBinContent(b+1, -1.*(relativeErrorOutUp->GetBinContent(b+1)-1));
2999 relativeErrorOutUp->SetBinError(b+1, .0);
3000 relativeErrorInLow->SetBinContent(b+1, -1.*(relativeErrorInLow->GetBinContent(b+1)-1));
3001 relativeErrorInLow->SetBinError(b+1, 0.);
3002 relativeErrorOutLow->SetBinContent(b+1, -1.*(relativeErrorOutLow->GetBinContent(b+1)-1));
3003 relativeErrorOutLow->SetBinError(b+1, .0);
3005 // these guys are already stored as percentages, so no need to remove the offset of 1
3006 // RMS is defined as sqrt(sum(squared))/N
3007 // min is defined as negative, max is defined as positive
3009 if(relativeErrorInUpN[b] < 1) relativeErrorInUpN[b] = 1;
3010 if(relativeErrorInLowN[b] < 1) relativeErrorInLowN[b] = 1;
3011 if(relativeErrorOutUpN[b] < 1) relativeErrorOutUpN[b] = 1;
3012 if(relativeErrorOutLowN[b] < 1) relativeErrorOutLowN[b] = 1;
3013 relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(b+1)/relativeErrorInUpN[b]));
3014 relativeErrorInUp->SetBinError(b+1, 0.);
3015 relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorOutUp->GetBinContent(b+1)/relativeErrorOutUpN[b]));
3016 relativeErrorOutUp->SetBinError(b+1, .0);
3017 relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(relativeErrorInLow->GetBinContent(b+1)/relativeErrorInLowN[b]));
3018 relativeErrorInLow->SetBinError(b+1, 0.);
3019 relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(relativeErrorOutLow->GetBinContent(b+1)/relativeErrorOutLowN[b]));
3020 relativeErrorOutLow->SetBinError(b+1, .0);
3021 } else if (fSymmRMS) {
3022 if(relativeErrorInUpN[b] < 1) relativeErrorInUpN[b] = 1;
3023 if(relativeErrorOutUpN[b] < 1) relativeErrorOutUpN[b] = 1;
3024 relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(b+1)/relativeErrorInUpN[b]));
3025 relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorOutUp->GetBinContent(b+1)/relativeErrorOutUpN[b]));
3029 relativeErrorInUp->SetYTitle("relative uncertainty");
3030 relativeErrorOutUp->SetYTitle("relative uncertainty");
3031 relativeErrorInLow->SetYTitle("relative uncertainty");
3032 relativeErrorOutLow->SetYTitle("relative uncertainty");
3033 relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
3034 relativeErrorInLow->GetYaxis()->SetRangeUser(-1.5, 3.);
3035 relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
3036 relativeErrorOutLow->GetYaxis()->SetRangeUser(-1.5, 3.);
3038 canvasNominalMasterIn->cd(2);
3039 Style(gPad, "GRID");
3040 Style(relativeErrorInUp, kBlue, kBar);
3041 Style(relativeErrorInLow, kGreen, kBar);
3042 relativeErrorInUp->DrawCopy("b");
3043 relativeErrorInLow->DrawCopy("same b");
3044 Style(AddLegend(gPad));
3045 SavePadToPDF(canvasNominalMasterIn);
3046 canvasNominalMasterIn->Write();
3047 canvasNominalMasterOut->cd(2);
3048 Style(gPad, "GRID");
3049 Style(relativeErrorOutUp, kBlue, kBar);
3050 Style(relativeErrorOutLow, kGreen, kBar);
3051 relativeErrorOutUp->DrawCopy("b");
3052 relativeErrorOutLow->DrawCopy("same b");
3053 Style(AddLegend(gPad));
3054 SavePadToPDF(canvasNominalMasterOut);
3055 canvasNominalMasterOut->Write();
3057 //_____________________________________________________________________________
3058 void AliJetFlowTools::PostProcess(TString def, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
3060 // go through the output file and perform post processing routines
3061 // can either be performed in one go with the unfolding, or at a later stage
3062 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
3063 TFile readMe(in.Data(), "READ"); // open file read-only
3064 if(readMe.IsZombie()) {
3065 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
3068 printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
3070 TList* listOfKeys((TList*)readMe.GetListOfKeys());
3072 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
3075 // prepare necessary canvasses
3076 TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
3077 TCanvas* canvasOut(0x0);
3078 if(fDphiUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
3079 TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
3080 TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
3081 if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
3082 TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
3083 TCanvas* canvasSpectraOut(0x0);
3084 if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
3085 TCanvas* canvasRatio(0x0);
3086 if(fDphiUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
3087 TCanvas* canvasV2(0x0);
3088 if(fDphiUnfolding) canvasV2 = new TCanvas("V2", "V2");
3089 TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
3090 TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
3091 TCanvas* canvasMasterOut(0x0);
3092 if(fDphiUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
3093 (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
3094 TDirectoryFile* defDir(0x0);
3095 TCanvas* canvasProfiles(new TCanvas("canvasProfiles", "canvasProfiles"));
3096 canvasProfiles->Divide(2);
3097 TProfile* ratioProfile(new TProfile("ratioProfile", "ratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3098 TProfile* v2Profile(new TProfile("v2Profile", "v2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3099 // get an estimate of the number of outputs and find the default set
3101 for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
3102 if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
3103 if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3107 Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%columns)>0));
3108 canvasIn->Divide(columns, rows);
3109 if(canvasOut) canvasOut->Divide(columns, rows);
3110 canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
3111 if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
3112 canvasSpectraIn->Divide(columns, rows);
3113 if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
3114 if(canvasRatio) canvasRatio->Divide(columns, rows);
3115 if(canvasV2) canvasV2->Divide(columns, rows);
3117 canvasMasterIn->Divide(columns, rows);
3118 if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
3119 // extract the default output
3120 TH1D* defaultUnfoldedJetSpectrumIn(0x0);
3121 TH1D* defaultUnfoldedJetSpectrumOut(0x0);
3123 TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
3124 TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
3125 if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
3126 if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
3127 printf(" > succesfully extracted default results < \n");
3130 // loop through the directories, only plot the graphs if the deconvolution converged
3131 TDirectoryFile* tempDir(0x0);
3132 TDirectoryFile* tempIn(0x0);
3133 TDirectoryFile* tempOut(0x0);
3134 for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
3135 // read the directory by its name
3136 tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3137 if(!tempDir) continue;
3138 TString dirName(tempDir->GetName());
3139 // try to read the in- and out of plane subdirs
3140 tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
3141 tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
3143 if(!tempIn) { // attempt to get MakeAU output
3144 TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
3145 TCanvas* tempPearson(new TCanvas(Form("pearson_%s", dirName.Data()), Form("pearson_%s", dirName.Data())));
3146 tempPearson->Divide(4,2);
3147 TCanvas* tempRatio(new TCanvas(Form("ratio_%s", dirName.Data()), Form("ratio_%s", dirName.Data())));
3148 tempRatio->Divide(4,2);
3149 TCanvas* tempResult(new TCanvas(Form("result_%s", dirName.Data()), Form("result_%s", dirName.Data())));
3150 tempResult->Divide(4,2);
3151 for(Int_t q(0); q < 8; q++) {
3152 tempIn = (TDirectoryFile*)tempDir->Get(Form("%s___%s", stringArray[q].Data(), dirName.Data()));
3154 // to see if the unfolding converged try to extract the pearson coefficients
3155 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
3157 printf(" - %s in plane converged \n", dirName.Data());
3158 tempPearson->cd(1+q);
3159 Style(gPad, "PEARSON");
3160 pIn->DrawCopy("colz");
3161 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3164 printf(" > found RatioRefoldedMeasured < \n");
3166 rIn->SetFillColor(kRed);
3169 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
3170 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
3171 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
3172 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
3173 if(dvector && avalue && rm && eff) {
3179 Style(gPad, "SPECTRUM");
3180 dvector->DrawCopy();
3182 Style(gPad, "SPECTRUM");
3185 Style(gPad, "PEARSON");
3186 rm->DrawCopy("colz");
3188 Style(gPad, "GRID");
3190 } else if(rm && eff) {
3194 Style(gPad, "PEARSON");
3195 rm->DrawCopy("colz");
3197 Style(gPad, "GRID");
3201 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
3202 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
3203 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
3204 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3205 if(defaultUnfoldedJetSpectrumIn) {
3206 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
3207 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
3208 temp->Divide(unfoldedSpectrum);
3209 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3210 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3211 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3212 canvasMasterIn->cd(j);
3213 temp->GetYaxis()->SetRangeUser(0., 2);
3216 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
3217 tempResult->cd(q+1);
3219 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3220 unfoldedSpectrum->DrawCopy();
3221 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3222 inputSpectrum->DrawCopy("same");
3223 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3224 refoldedSpectrum->DrawCopy("same");
3225 TLegend* l(AddLegend(gPad));
3227 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3228 Float_t chi(fitStatus->GetBinContent(1));
3229 Float_t pen(fitStatus->GetBinContent(2));
3230 Int_t dof((int)fitStatus->GetBinContent(3));
3231 Float_t beta(fitStatus->GetBinContent(4));
3232 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3233 } else if (fitStatus) { // only available in SVD fit
3234 Int_t reg((int)fitStatus->GetBinContent(1));
3235 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3242 // to see if the unfolding converged try to extract the pearson coefficients
3243 TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
3245 printf(" - %s in plane converged \n", dirName.Data());
3247 Style(gPad, "PEARSON");
3248 pIn->DrawCopy("colz");
3249 TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3252 printf(" > found RatioRefoldedMeasured < \n");
3253 canvasRatioMeasuredRefoldedIn->cd(j);
3254 rIn->SetFillColor(kRed);
3257 TH1D* dvector((TH1D*)tempIn->Get("dVector"));
3258 TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
3259 TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
3260 TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
3261 if(dvector && avalue && rm && eff) {
3267 Style(gPad, "SPECTRUM");
3268 dvector->DrawCopy();
3270 Style(gPad, "SPECTRUM");
3273 Style(gPad, "PEARSON");
3274 rm->DrawCopy("colz");
3276 Style(gPad, "GRID");
3278 } else if(rm && eff) {
3282 Style(gPad, "PEARSON");
3283 rm->DrawCopy("colz");
3285 Style(gPad, "GRID");
3289 TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
3290 TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
3291 TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
3292 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3293 if(defaultUnfoldedJetSpectrumIn) {
3294 Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
3295 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
3296 temp->Divide(unfoldedSpectrum);
3297 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3298 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3299 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3300 canvasMasterIn->cd(j);
3301 temp->GetYaxis()->SetRangeUser(0., 2);
3304 TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
3305 canvasSpectraIn->cd(j);
3307 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3308 unfoldedSpectrum->DrawCopy();
3309 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3310 inputSpectrum->DrawCopy("same");
3311 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3312 refoldedSpectrum->DrawCopy("same");
3313 TLegend* l(AddLegend(gPad));
3315 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3316 Float_t chi(fitStatus->GetBinContent(1));
3317 Float_t pen(fitStatus->GetBinContent(2));
3318 Int_t dof((int)fitStatus->GetBinContent(3));
3319 Float_t beta(fitStatus->GetBinContent(4));
3320 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3321 } else if (fitStatus) { // only available in SVD fit
3322 Int_t reg((int)fitStatus->GetBinContent(1));
3323 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3328 TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
3330 printf(" - %s out of plane converged \n", dirName.Data());
3332 Style(gPad, "PEARSON");
3333 pOut->DrawCopy("colz");
3334 TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3337 printf(" > found RatioRefoldedMeasured < \n");
3338 canvasRatioMeasuredRefoldedOut->cd(j);
3339 rOut->SetFillColor(kRed);
3342 TH1D* dvector((TH1D*)tempOut->Get("dVector"));
3343 TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
3344 TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
3345 TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
3346 if(dvector && avalue && rm && eff) {
3352 Style(gPad, "SPECTRUM");
3353 dvector->DrawCopy();
3355 Style(gPad, "SPECTRUM");
3358 Style(gPad, "PEARSON");
3359 rm->DrawCopy("colz");
3361 Style(gPad, "GRID");
3363 } else if(rm && eff) {
3367 Style(gPad, "PEARSON");
3368 rm->DrawCopy("colz");
3370 Style(gPad, "GRID");
3374 TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
3375 TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
3376 TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
3377 if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3378 if(defaultUnfoldedJetSpectrumOut) {
3379 Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
3380 TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirName.Data())));
3381 temp->Divide(unfoldedSpectrum);
3382 temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3383 temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3384 temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3385 canvasMasterOut->cd(j);
3386 temp->GetYaxis()->SetRangeUser(0., 2.);
3389 TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
3390 canvasSpectraOut->cd(j);
3392 Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3393 unfoldedSpectrum->DrawCopy();
3394 Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3395 inputSpectrum->DrawCopy("same");
3396 Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3397 refoldedSpectrum->DrawCopy("same");
3398 TLegend* l(AddLegend(gPad));
3400 if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3401 Float_t chi(fitStatus->GetBinContent(1));
3402 Float_t pen(fitStatus->GetBinContent(2));
3403 Int_t dof((int)fitStatus->GetBinContent(3));
3404 Float_t beta(fitStatus->GetBinContent(4));
3405 l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3406 } else if (fitStatus) { // only available in SVD fit
3407 Int_t reg((int)fitStatus->GetBinContent(1));
3408 l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3412 if(canvasRatio && canvasV2) {
3414 TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
3415 Double_t _tempx(0), _tempy(0);
3418 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3419 ratioYield->GetPoint(b, _tempx, _tempy);
3420 ratioProfile->Fill(_tempx, _tempy);
3422 ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
3423 ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3424 ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
3425 ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3426 ratioYield->SetFillColor(kRed);
3427 ratioYield->Draw("ap");
3430 TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
3433 for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3434 ratioV2->GetPoint(b, _tempx, _tempy);
3435 v2Profile->Fill(_tempx, _tempy);
3438 v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
3439 v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3440 ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
3441 ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3442 ratioV2->SetFillColor(kRed);
3443 ratioV2->Draw("ap");
3447 TFile output(out.Data(), "RECREATE");
3448 canvasProfiles->cd(1);
3449 Style(ratioProfile);
3450 ratioProfile->DrawCopy();
3451 canvasProfiles->cd(2);
3453 v2Profile->DrawCopy();
3454 SavePadToPDF(canvasProfiles);
3455 canvasProfiles->Write();
3456 SavePadToPDF(canvasIn);
3459 SavePadToPDF(canvasOut);
3462 SavePadToPDF(canvasRatioMeasuredRefoldedIn);
3463 canvasRatioMeasuredRefoldedIn->Write();
3464 if(canvasRatioMeasuredRefoldedOut) {
3465 SavePadToPDF(canvasRatioMeasuredRefoldedOut);
3466 canvasRatioMeasuredRefoldedOut->Write();
3468 SavePadToPDF(canvasSpectraIn);
3469 canvasSpectraIn->Write();
3470 if(canvasSpectraOut) {
3471 SavePadToPDF(canvasSpectraOut);
3472 canvasSpectraOut->Write();
3473 SavePadToPDF(canvasRatio);
3474 canvasRatio->Write();
3475 SavePadToPDF(canvasV2);
3478 SavePadToPDF(canvasMasterIn);
3479 canvasMasterIn->Write();
3480 if(canvasMasterOut) {
3481 SavePadToPDF(canvasMasterOut);
3482 canvasMasterOut->Write();
3484 SavePadToPDF(canvasMISC);
3485 canvasMISC->Write();
3489 //_____________________________________________________________________________
3490 void AliJetFlowTools::BootstrapSpectra(TString def, TString in, TString out) const
3492 // function to interpret results of bootstrapping routine
3493 // TString def should hold the true emperical distribution
3494 if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
3495 TFile readMe(in.Data(), "READ"); // open file read-only
3496 if(readMe.IsZombie()) {
3497 printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
3500 printf("\n\n\n\t\t BootstrapSpectra \n > Recovered the following file structure : \n <");
3502 TList* listOfKeys((TList*)readMe.GetListOfKeys());
3504 printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
3507 // get an estimate of the number of outputs and find the default set
3508 TDirectoryFile* defDir(0x0);
3509 for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
3510 if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
3511 if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3515 // extract the default output this is the 'emperical' distribution
3516 // w.r.t. which the other distributions will be evaluated
3517 TH1D* defaultUnfoldedJetSpectrumIn(0x0);
3518 TH1D* defaultUnfoldedJetSpectrumOut(0x0);
3519 TH1D* defaultInputSpectrumIn(0x0);
3520 TH1D* defaultInputSpectrumOut(0x0);
3521 TGraphErrors* v2Emperical(0x0);
3523 TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
3524 TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
3526 defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
3527 defaultInputSpectrumIn = (TH1D*)defDirIn->Get(Form("InputSpectrum_in_%s", def.Data()));
3530 defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
3531 defaultInputSpectrumOut = (TH1D*)defDirOut->Get(Form("InputSpectrum_out_%s", def.Data()));
3534 if((!defaultUnfoldedJetSpectrumIn && defaultUnfoldedJetSpectrumOut)) {
3535 printf(" BootstrapSpectra: couldn't extract default spectra, aborting! \n");
3538 else v2Emperical = GetV2(defaultUnfoldedJetSpectrumIn, defaultUnfoldedJetSpectrumOut, fEventPlaneRes);
3540 // now that we know for sure that the input is in place, reserve the bookkeeping histograms
3541 TH1F* delta0[fBinsTrue->GetSize()-1];
3542 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3543 delta0[i] = new TH1F(Form("delta%i_%.2f-%.2f_gev", i, fBinsTrue->At(i), fBinsTrue->At(i+1)), Form("#Delta_{0, %i} p_{T} %.2f-%.2f GeV/c", i, fBinsTrue->At(i), fBinsTrue->At(i+1)), 30, -1., 1.);//.15, .15);
3546 // and a canvas for illustration purposes only
3547 TCanvas* spectraIn(new TCanvas("spectraIn", "spectraIn"));
3548 TCanvas* spectraOut(new TCanvas("spectraOut", "spectraOut"));
3549 // common reference (in this case the generated v2)
3550 TF1* commonReference = new TF1("v2_strong_bkg_comb", "(x<=3)*.07+(x>3&&x<5)*(.07-(x-3)*.035)+(x>30&&x<40)*(x-30)*.005+(x>40)*.05", 0, 200);
3552 // loop through the directories, only plot the graphs if the deconvolution converged
3553 TDirectoryFile* tempDir(0x0);
3554 TDirectoryFile* tempIn(0x0);
3555 TDirectoryFile* tempOut(0x0);
3556 TH1D* unfoldedSpectrumIn(0x0);
3557 TH1D* unfoldedSpectrumOut(0x0);
3558 TH1D* measuredSpectrumIn(0x0);
3559 TH1D* measuredSpectrumOut(0x0);
3560 for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
3561 // read the directory by its name
3562 tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3563 if(!tempDir) continue;
3564 TString dirName(tempDir->GetName());
3565 // read only bootstrapped outputs
3566 if(!dirName.Contains("bootstrap")) continue;
3567 // try to read the in- and out of plane subdirs
3568 tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
3569 tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
3571 // extract the unfolded spectra only if both in- and out-of-plane converted (i.e. pearson coefficients were saved)
3573 if(!(TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data()))) continue;
3574 unfoldedSpectrumIn = (TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data()));
3575 measuredSpectrumIn = (TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data()));
3577 (j==1) ? measuredSpectrumIn->DrawCopy() : measuredSpectrumIn->DrawCopy("same");
3580 if(!(TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data()))) continue;
3581 unfoldedSpectrumOut = (TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data()));
3582 measuredSpectrumOut = (TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data()));
3584 (j==1) ? measuredSpectrumOut->DrawCopy() : measuredSpectrumOut->DrawCopy("same");
3586 // get v2 with statistical uncertainties from the extracted spectra
3587 TGraphErrors* v2Bootstrapped(GetV2(unfoldedSpectrumIn, unfoldedSpectrumOut, fEventPlaneRes));
3588 // and loop over all bins to fill the bookkeeping histograms
3589 Double_t yBoot(0), yEmp(0), xDud(0);
3590 // optional: common reference (in this case the sampled v2 value)
3591 for(Int_t k(0); k < fBinsTrue->GetSize()-1; k++) {
3592 // read values point by point (passed by reference)
3593 v2Bootstrapped->GetPoint(k, xDud, yBoot);
3594 v2Emperical->GetPoint(k, xDud, yEmp);
3595 if(commonReference) {
3596 if(!commonReference->Eval(xDud)<1e-10) {
3597 yEmp/=commonReference->Eval(xDud);
3598 yBoot/=commonReference->Eval(xDud);
3599 } else { // if reference equals zero, take emperical distribution as reference
3604 cout << " yEmp " << yEmp << " yBoot " << yBoot << endl;
3605 // save relative difference per pt bin
3606 if(TMath::Abs(yBoot)>1e-10) delta0[k]->Fill(yEmp - yBoot);
3609 // extracting final results now, as first estimate just a gaus fit to the distributions
3610 // (should be changed perhaps to proper rms eventually)
3611 // attach relevant data to current buffer in the same loop
3612 TFile output(out.Data(), "RECREATE");
3613 defaultInputSpectrumIn->SetLineColor(kRed);
3615 defaultInputSpectrumIn->DrawCopy("same");
3616 defaultInputSpectrumOut->SetLineColor(kRed);
3618 defaultInputSpectrumOut->DrawCopy("same");
3620 spectraOut->Write();
3622 TH1F* delta0spread(new TH1F("delta0spread", "#sigma(#Delta_{0})", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3623 TH1F* unfoldingError(new TH1F("unfoldingError", "error from unfolding algorithm", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3624 TF1* gaus(new TF1("gaus", "gaus"/*[0]*exp(-0.5*((x-[1])/[2])**2)"*/, -1., 1.));
3625 Double_t xDud(0), yDud(0);
3626 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3627 delta0[i]->Fit(gaus);
3628 delta0[i]->GetYaxis()->SetTitle("counts");
3629 delta0[i]->GetXaxis()->SetTitle("(v_{2, jet}^{measured} - v_{2, jet}^{generated}) / input v_{2}");
3630 delta0spread->SetBinContent(i+1, gaus->GetParameter(1)); // mean of gaus
3631 delta0spread->SetBinError(i+1, gaus->GetParameter(2)); // sigma of gaus
3632 cout << " mean " << gaus->GetParameter(1) << endl;
3633 cout << " sigm " << gaus->GetParameter(2) << endl;
3635 v2Emperical->GetPoint(i, xDud, yDud);
3636 unfoldingError->SetBinContent(i+1, 1e-10/*gaus->GetParameter(1)*/);
3637 if(commonReference && !commonReference->Eval(xDud)<1e-10) unfoldingError->SetBinError(i+1, v2Emperical->GetErrorY(i)/(commonReference->Eval(xDud)));
3638 else if(yDud>10e-10) unfoldingError->SetBinError(i+1, v2Emperical->GetErrorY(i)/yDud);
3639 else unfoldingError->SetBinError(i+1, 0.);
3641 delta0spread->GetXaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
3642 delta0spread->GetYaxis()->SetTitle("(mean v_{2, jet}^{measured} - v_{2, jet}^{generated}) / input v_{2}");
3643 delta0spread->Write();
3644 unfoldingError->GetXaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
3645 unfoldingError->GetYaxis()->SetTitle("(mean v_{2, jet}^{measured} - v_{2, jet}^{generated}) / input v_{2}");
3646 unfoldingError->Write();
3647 // write the buffer and close the file
3651 //_____________________________________________________________________________
3652 Bool_t AliJetFlowTools::SetRawInput (
3653 TH2D* detectorResponse, // detector response matrix
3654 TH1D* jetPtIn, // in plane jet spectrum
3655 TH1D* jetPtOut, // out of plane jet spectrum
3656 TH1D* dptIn, // in plane delta pt distribution
3657 TH1D* dptOut, // out of plane delta pt distribution
3659 // set input histograms manually
3660 fDetectorResponse = detectorResponse;
3661 fSpectrumIn = jetPtIn;
3662 fSpectrumOut = jetPtOut;
3664 fDptOutDist = dptOut;
3665 fRawInputProvided = kTRUE;
3666 // check if all data is provided
3667 if(!fDetectorResponse) {
3668 printf(" fDetectorResponse not found \n ");
3671 // check if the pt bin for true and rec have been set
3672 if(!fBinsTrue || !fBinsRec) {
3673 printf(" No true or rec bins set, please set binning ! \n");
3676 if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
3677 // procedures, these profiles will be nonsensical, user is responsible
3678 fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3679 fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3680 fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3682 // normalize spectra to event count if requested
3683 if(fNormalizeSpectra) {
3684 fEventCount = eventCount;
3685 if(fEventCount > 0) {
3686 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
3687 fSpectrumOut->Sumw2();
3688 fSpectrumIn->Scale(1./((double)fEventCount));
3689 fSpectrumOut->Scale(1./((double)fEventCount));
3692 if(!fNormalizeSpectra && fEventCount > 0) {
3693 fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
3694 fSpectrumOut->Sumw2();
3695 fSpectrumIn->Scale(1./((double)fEventCount));
3696 fSpectrumOut->Scale(1./((double)fEventCount));
3698 fDptIn = ConstructDPtResponseFromTH1D(fDptInDist, fAvoidRoundingError);
3699 fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
3700 fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
3701 fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
3702 fDptOut = ConstructDPtResponseFromTH1D(fDptOutDist, fAvoidRoundingError);
3703 fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
3704 fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
3705 fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
3709 //_____________________________________________________________________________
3710 TGraphErrors* AliJetFlowTools::GetRatio(TH1 *h1, TH1* h2, TString name, Bool_t appendFit, Int_t xmax)
3712 // return ratio of h1 / h2
3713 // histograms can have different binning. errors are propagated as uncorrelated
3715 printf(" GetRatio called with NULL argument(s) \n ");
3719 TGraphErrors *gr = new TGraphErrors();
3720 gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3721 Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
3722 TH1* dud((TH1*)h1->Clone("dud"));
3726 if(!dud->Divide(h1, h2)) {
3727 printf(" > ROOT failed to divide two histogams, dividing manually \n < ");
3728 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3729 binCent = h1->GetXaxis()->GetBinCenter(i);
3730 if(xmax > 0. && binCent > xmax) continue;
3731 j = h2->FindBin(binCent);
3732 binWidth = h1->GetXaxis()->GetBinWidth(i);
3733 if(h2->GetBinContent(j) > 0.) {
3734 ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
3735 Double_t A = h1->GetBinError(i)/h1->GetBinContent(i);
3736 Double_t B = h2->GetBinError(i)/h2->GetBinContent(i);
3737 error2 = ratio*ratio*A*A+ratio*ratio*B*B;
3738 if(error2 > 0 ) error2 = TMath::Sqrt(error2);
3739 gr->SetPoint(i-1, binCent, ratio);
3740 gr->SetPointError(i-1, 0.5*binWidth, error2);
3744 printf( " > using ROOT to divide two histograms < \n");
3745 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3746 binCent = dud->GetXaxis()->GetBinCenter(i);
3747 if(xmax > 0. && binCent > xmax) continue;
3748 binWidth = dud->GetXaxis()->GetBinWidth(i);
3749 gr->SetPoint(i-1,binCent,dud->GetBinContent(i));
3750 gr->SetPointError(i-1, 0.5*binWidth,dud->GetBinError(i));
3755 TF1* fit(new TF1("lin", "pol0", 10, 100));
3758 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
3762 //_____________________________________________________________________________
3763 TGraphErrors* AliJetFlowTools::GetV2(TH1 *h1, TH1* h2, Double_t r, TString name)
3765 // get v2 from difference of in plane, out of plane yield
3766 // h1 must hold the in-plane yield, h2 holds the out of plane yield
3767 // r is the event plane resolution for the chosen centrality
3769 printf(" GetV2 called with NULL argument(s) \n ");
3772 TGraphErrors *gr = new TGraphErrors();
3773 gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3774 Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
3775 Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
3777 for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
3778 binCent = h1->GetXaxis()->GetBinCenter(i);
3779 binWidth = h1->GetXaxis()->GetBinWidth(i);
3780 if(h2->GetBinContent(i) > 0.) {
3781 in = h1->GetBinContent(i);
3782 ein = h1->GetBinError(i);
3783 out = h2->GetBinContent(i);
3784 eout = h2->GetBinError(i);
3785 ratio = pre*((in-out)/(in+out));
3786 error2 = (4.*out*out/(TMath::Power(in+out, 4)))*ein*ein+(4.*in*in/(TMath::Power(in+out, 4)))*eout*eout;
3787 error2 = error2*pre*pre;
3788 if(error2 > 0) error2 = TMath::Sqrt(error2);
3789 gr->SetPoint(i-1,binCent,ratio);
3790 gr->SetPointError(i-1,0.5*binWidth,error2);
3793 if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
3796 //_____________________________________________________________________________
3797 TGraphAsymmErrors* AliJetFlowTools::GetV2WithSystematicErrors(
3798 TH1* h1, TH1* h2, Double_t r, TString name,
3799 TH1* relativeErrorInUp,
3800 TH1* relativeErrorInLow,
3801 TH1* relativeErrorOutUp,
3802 TH1* relativeErrorOutLow,
3805 // get v2 with asymmetric systematic error
3806 // note that this is ONLY the systematic error, no statistical error!
3807 // rho is the pearson correlation coefficient
3808 TGraphErrors* tempV2(GetV2(h1, h2, r, name));
3809 Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
3810 Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
3811 Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
3812 Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
3813 Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
3814 Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
3815 Double_t in(0.), out(0.), einUp(0.), einLow(0.), eoutUp(0.), eoutLow(0.), error2Up(0.), error2Low(0.);
3816 // loop through the bins and do error propagation
3817 for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3818 // extract the absolute errors
3819 in = h1->GetBinContent(i+1);
3820 einUp = TMath::Abs(in*relativeErrorInUp->GetBinContent(i+1));
3821 einLow = TMath::Abs(in*relativeErrorInLow->GetBinContent(1+i));
3822 out = h2->GetBinContent(i+1);
3823 eoutUp = TMath::Abs(out*relativeErrorOutUp->GetBinContent(1+i));
3824 eoutLow = TMath::Abs(out*relativeErrorOutLow->GetBinContent(1+i));
3825 // get the error squared
3827 error2Up = TMath::Power(((r*4.)/(TMath::Pi())),-2.)*((4.*out*out/(TMath::Power(in+out, 4)))*einUp*einUp+(4.*in*in/(TMath::Power(in+out, 4)))*eoutLow*eoutLow);
3828 error2Low =TMath::Power(((r*4.)/(TMath::Pi())),-2.)*((4.*out*out/(TMath::Power(in+out, 4)))*einLow*einLow+(4.*in*in/(TMath::Power(in+out, 4)))*eoutUp*eoutUp);
3830 error2Up = TMath::Power(((r*4.)/(TMath::Pi())),-2.)*((4.*out*out/(TMath::Power(in+out, 4)))*einUp*einUp+(4.*in*in/(TMath::Power(in+out, 4)))*eoutUp*eoutUp-((8.*out*in)/(TMath::Power(in+out, 4)))*rho*einUp*eoutUp);
3831 error2Low =TMath::Power(((r*4.)/(TMath::Pi())),-2.)*((4.*out*out/(TMath::Power(in+out, 4)))*einLow*einLow+(4.*in*in/(TMath::Power(in+out, 4)))*eoutLow*eoutLow-((8.*out*in)/(TMath::Power(in+out, 4)))*rho*einLow*eoutLow);
3833 if(error2Up > 0) error2Up = TMath::Sqrt(error2Up);
3834 if(error2Low > 0) error2Low = TMath::Sqrt(error2Low);
3838 // get the bin width (which is the 'error' on x)
3839 Double_t binWidth(h1->GetBinWidth(i+1));
3840 axl[i] = binWidth/2.;
3841 axh[i] = binWidth/2.;
3842 // now get the coordinate for the poin
3843 tempV2->GetPoint(i, ax[i], ay[i]);
3845 // save the nominal ratio
3846 TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
3847 nominalError->SetName("v_{2} with shape uncertainty");
3848 // do some memory management
3857 return nominalError;
3859 //_____________________________________________________________________________
3860 void AliJetFlowTools::WriteObject(TObject* object, TString suffix, Bool_t kill) {
3861 // write object, if a unique identifier is given the object is cloned
3862 // and the clone is saved. setting kill to true will delete the original obect from the heap
3864 printf(" > WriteObject:: called with NULL arguments \n ");
3866 } else if(!strcmp("", suffix.Data())) object->Write();
3868 TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
3871 if(kill) delete object;
3873 //_____________________________________________________________________________
3874 TH2D* AliJetFlowTools::ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError) {
3875 // construt a delta pt response matrix from supplied dpt distribution
3876 // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
3877 // do a weighted rebinning to a (coarser) dpt distribution
3878 // be careful with the binning of the dpt response: it should be equal to that
3879 // of the response matrix, otherwise dpt and response matrices cannot be multiplied
3881 // the response matrix will be square and have the same binning
3882 // (min, max and granularity) of the input histogram
3883 Int_t bins(dpt->GetXaxis()->GetNbins()); // number of bins, will also be no of rows, columns
3884 Double_t _bins[bins+1]; // prepare array with bin borders
3885 for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
3886 _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1); // get upper edge
3887 TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
3888 for(Int_t j(0); j < bins+1; j++) { // loop on pt true slices j
3889 Bool_t skip = kFALSE;
3890 for(Int_t k(0); k < bins+1; k++) { // loop on pt gen slices k
3891 (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
3892 if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
3897 //_____________________________________________________________________________
3898 TH2D* AliJetFlowTools::GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix) {
3899 if(!binsTrue || !binsRec) {
3900 printf(" > GetUnityResponse:: function called with NULL arguments < \n");
3903 TString name(Form("unityResponse_%s", suffix.Data()));
3904 TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
3905 for(Int_t i(0); i < binsTrue->GetSize(); i++) {
3906 for(Int_t j(0); j < binsRec->GetSize(); j++) {
3907 if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
3912 //_____________________________________________________________________________
3913 void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
3914 // save configuration parameters to histogram
3915 TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
3916 summary->SetBinContent(1, fBetaIn);
3917 summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
3918 summary->SetBinContent(2, fBetaOut);
3919 summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
3920 summary->SetBinContent(3, fCentralityArray->At(0));
3921 summary->GetXaxis()->SetBinLabel(3, "fCentralityArray[0]");
3922 summary->SetBinContent(4, (int)convergedIn);
3923 summary->GetXaxis()->SetBinLabel(4, "convergedIn");
3924 summary->SetBinContent(5, (int)convergedOut);
3925 summary->GetXaxis()->SetBinLabel(5, "convergedOut");
3926 summary->SetBinContent(6, (int)fAvoidRoundingError);
3927 summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
3928 summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
3929 summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
3930 summary->SetBinContent(8, (int)fPrior);
3931 summary->GetXaxis()->SetBinLabel(8, "fPrior");
3932 summary->SetBinContent(9, fSVDRegIn);
3933 summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
3934 summary->SetBinContent(10, fSVDRegOut);
3935 summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
3936 summary->SetBinContent(11, (int)fSVDToy);
3937 summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
3938 summary->SetBinContent(12, fJetRadius);
3939 summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
3940 summary->SetBinContent(13, (int)fNormalizeSpectra);
3941 summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
3942 summary->SetBinContent(14, (int)fSmoothenPrior);
3943 summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
3944 summary->SetBinContent(15, (int)fTestMode);
3945 summary->GetXaxis()->SetBinLabel(15, "fTestMode");
3946 summary->SetBinContent(16, (int)fUseDetectorResponse);
3947 summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
3948 summary->SetBinContent(17, fBayesianIterIn);
3949 summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
3950 summary->SetBinContent(18, fBayesianIterOut);
3951 summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
3952 summary->SetBinContent(19, fBayesianSmoothIn);
3953 summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
3954 summary->SetBinContent(20, fBayesianSmoothOut);
3955 summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
3957 //_____________________________________________________________________________
3958 void AliJetFlowTools::ResetAliUnfolding() {
3959 // ugly function: reset all unfolding parameters
3960 TVirtualFitter* fitter(TVirtualFitter::GetFitter());
3962 printf(" > Found fitter, will delete it < \n");
3966 printf(" > Found gMinuit, will re-create it < \n");
3968 gMinuit = new TMinuit;
3970 AliUnfolding::fgCorrelationMatrix = 0;
3971 AliUnfolding::fgCorrelationMatrixSquared = 0;
3972 AliUnfolding::fgCorrelationCovarianceMatrix = 0;
3973 AliUnfolding::fgCurrentESDVector = 0;
3974 AliUnfolding::fgEntropyAPriori = 0;
3975 AliUnfolding::fgEfficiency = 0;
3976 AliUnfolding::fgUnfoldedAxis = 0;
3977 AliUnfolding::fgMeasuredAxis = 0;
3978 AliUnfolding::fgFitFunction = 0;
3979 AliUnfolding::fgMaxInput = -1;
3980 AliUnfolding::fgMaxParams = -1;
3981 AliUnfolding::fgOverflowBinLimit = -1;
3982 AliUnfolding::fgRegularizationWeight = 10000;
3983 AliUnfolding::fgSkipBinsBegin = 0;
3984 AliUnfolding::fgMinuitStepSize = 0.1;
3985 AliUnfolding::fgMinuitPrecision = 1e-6;
3986 AliUnfolding::fgMinuitMaxIterations = 1000000;
3987 AliUnfolding::fgMinuitStrategy = 1.;
3988 AliUnfolding::fgMinimumInitialValue = kFALSE;
3989 AliUnfolding::fgMinimumInitialValueFix = -1;
3990 AliUnfolding::fgNormalizeInput = kFALSE;
3991 AliUnfolding::fgNotFoundEvents = 0;
3992 AliUnfolding::fgSkipBin0InChi2 = kFALSE;
3993 AliUnfolding::fgBayesianSmoothing = 1;
3994 AliUnfolding::fgBayesianIterations = 10;
3995 AliUnfolding::fgDebug = kFALSE;
3996 AliUnfolding::fgCallCount = 0;
3997 AliUnfolding::fgPowern = 5;
3998 AliUnfolding::fChi2FromFit = 0.;
3999 AliUnfolding::fPenaltyVal = 0.;
4000 AliUnfolding::fAvgResidual = 0.;
4001 AliUnfolding::fgPrintChi2Details = 0;
4002 AliUnfolding::fgCanvas = 0;
4003 AliUnfolding::fghUnfolded = 0;
4004 AliUnfolding::fghCorrelation = 0;
4005 AliUnfolding::fghEfficiency = 0;
4006 AliUnfolding::fghMeasured = 0;
4007 AliUnfolding::SetMinuitStepSize(1.);
4008 AliUnfolding::SetMinuitPrecision(1e-6);
4009 AliUnfolding::SetMinuitMaxIterations(100000);
4010 AliUnfolding::SetMinuitStrategy(2.);
4011 AliUnfolding::SetDebug(1);
4013 //_____________________________________________________________________________
4014 TH1D* AliJetFlowTools::ProtectHeap(TH1D* protect, Bool_t kill, TString suffix) const {
4015 // protect heap by adding unique qualifier to name
4016 if(!protect) return 0x0;
4017 TH1D* p = (TH1D*)protect->Clone();
4018 TString tempString(fActiveString);
4020 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
4021 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4022 if(kill) delete protect;
4025 //_____________________________________________________________________________
4026 TH2D* AliJetFlowTools::ProtectHeap(TH2D* protect, Bool_t kill, TString suffix) const {
4027 // protect heap by adding unique qualifier to name
4028 if(!protect) return 0x0;
4029 TH2D* p = (TH2D*)protect->Clone();
4030 TString tempString(fActiveString);
4032 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
4033 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4034 if(kill) delete protect;
4037 //_____________________________________________________________________________
4038 TGraphErrors* AliJetFlowTools::ProtectHeap(TGraphErrors* protect, Bool_t kill, TString suffix) const {
4039 // protect heap by adding unique qualifier to name
4040 if(!protect) return 0x0;
4041 TGraphErrors* p = (TGraphErrors*)protect->Clone();
4042 TString tempString(fActiveString);
4044 p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
4045 p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4046 if(kill) delete protect;
4049 //_____________________________________________________________________________
4050 void AliJetFlowTools::MakeAU() {
4051 // === azimuthal unfolding ===
4053 // unfolds the spectrum in delta phi bins, extracts the yield per bin, and does a fit
4054 // in transverse momentum and azimuthal correlation space to extract v2 params
4055 // settings are equal to the ones used for 'Make()'
4057 // basic steps that are followed:
4058 // 1) rebin the raw output of the jet task to the desired binnings
4059 // 2) calls the unfolding routine
4060 // 3) writes output to file
4061 // can be repeated multiple times with different configurations
4063 Int_t low[] = {1, 6, 11, 16, 21, 26, 31, 36};
4064 Int_t up[] = {5, 10, 15, 20, 25, 30, 35, 40};
4065 TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
4066 const Int_t ptBins(fBinsTrue->GetSize()-1);
4067 const Int_t dPhiBins(8);
4068 TH1D* dPtdPhi[fBinsTrue->GetSize()];
4069 for(Int_t i(0); i < ptBins; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), dPhiBins, 0, TMath::Pi());
4071 // for the output initialize a canvas
4072 TCanvas* v2Fits(new TCanvas("v2 fits", "v2 fits"));
4073 v2Fits->Divide(4, TMath::Floor((1+ptBins)/(float)4)+(((1+ptBins)%4)>0));
4075 for(Int_t i(0); i < dPhiBins; i++) {
4076 // 1) manipulation of input histograms
4077 // check if the input variables are present
4078 if(!PrepareForUnfolding(low[i], up[i])) return;
4079 // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
4080 // parts of the spectrum can end up in over or underflow bins
4081 TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, Form("resized_%s", stringArray[i].Data()), kFALSE);
4083 // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
4084 // the template will be used as a prior for the chi2 unfolding
4085 TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, stringArray[i], kFALSE);
4087 // get the full response matrix from the dpt and the detector response
4088 fDetectorResponse = NormalizeTH2D(fDetectorResponse);
4089 // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
4090 // so that unfolding should return the initial spectrum
4091 if(fUseDptResponse && fUseDetectorResponse) fFullResponseIn = MatrixMultiplication(fDptIn, fDetectorResponse);
4092 else if (fUseDptResponse && !fUseDetectorResponse) fFullResponseIn = fDptIn;
4093 else if (!fUseDptResponse && fUseDetectorResponse) fFullResponseIn = fDetectorResponse;
4094 else if (!fUseDptResponse && !fUseDetectorResponse && !fUnfoldingAlgorithm == AliJetFlowTools::kNone) return;
4095 // normalize each slide of the response to one
4096 NormalizeTH2D(fFullResponseIn);
4097 // resize to desired binning scheme
4098 TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, stringArray[i]);
4099 // get the kinematic efficiency
4100 TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
4101 kinematicEfficiencyIn->SetNameTitle(Form("kin_eff_%s", stringArray[i].Data()), Form("kin_eff_%s", stringArray[i].Data()));
4102 // suppress the errors
4103 for(Int_t j(0); j < kinematicEfficiencyIn->GetXaxis()->GetNbins(); j++) kinematicEfficiencyIn->SetBinError(1+j, 0.);
4104 TH1D* jetFindingEfficiency(0x0);
4105 if(fJetFindingEff) {
4106 jetFindingEfficiency = ProtectHeap(fJetFindingEff);
4107 jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
4108 jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
4110 // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
4111 TH1D* unfoldedJetSpectrumIn(0x0);
4112 fActiveDir->cd(); // select active dir
4113 TDirectoryFile* dirIn = new TDirectoryFile(Form("%s___%s", stringArray[i].Data(), fActiveString.Data()), Form("%s___%s", stringArray[i].Data(), fActiveString.Data()));
4114 dirIn->cd(); // select inplane subdir
4115 // select the unfolding method
4116 unfoldedJetSpectrumIn = UnfoldWrapper(
4117 measuredJetSpectrumIn,
4119 kinematicEfficiencyIn,
4120 measuredJetSpectrumTrueBinsIn,
4121 TString("dPtdPhiUnfolding"),
4122 jetFindingEfficiency);
4123 // arbitrarily save one of the full outputs (same for all dphi bins, avoid duplicates)
4125 resizedResponseIn->SetNameTitle(Form("ResponseMatrix_%s", stringArray[i].Data()), Form("response matrix %s", stringArray[i].Data()));
4126 resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
4127 resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
4128 resizedResponseIn = ProtectHeap(resizedResponseIn);
4129 resizedResponseIn->Write();
4130 kinematicEfficiencyIn->SetNameTitle(Form("KinematicEfficiency_%s", stringArray[i].Data()), Form("Kinematic efficiency, %s", stringArray[i].Data()));
4131 kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
4132 kinematicEfficiencyIn->Write();
4133 fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
4134 fDetectorResponse = ProtectHeap(fDetectorResponse, kFALSE);
4135 fDetectorResponse->Write();
4136 // optional histograms
4138 fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", Form("[INPUT] Jet spectrum, %s", stringArray[i].Data()));
4139 fSpectrumIn->Write();
4140 fDptInDist->SetNameTitle("[ORIG]DeltaPt", Form("#delta p_{T} distribution, %s", stringArray[i].Data()));
4141 fDptInDist->Write();
4142 fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix", Form("#delta p_{T} matrix, %s", stringArray[i].Data()));
4144 fFullResponseIn->SetNameTitle("ResponseMatrix", Form("Response matrix, %s", stringArray[i].Data()));
4145 fFullResponseIn->Write();
4149 fDeltaPtDeltaPhi->Write();
4150 fJetPtDeltaPhi->Write();
4152 TH1D* dud(ProtectHeap(unfoldedJetSpectrumIn, kTRUE, stringArray[i]));;
4153 Double_t integralError(0);
4154 // at this point in the code, the spectrum has been unfolded in a certain region of dPhi space
4155 // next step is splitting it in pt space as well to estimate the yield differentially in pt
4156 for(Int_t j(0); j < ptBins; j++) {
4157 // get the integrated
4158 Double_t integral(dud->IntegralAndError(j+1, j+2, integralError));
4159 dPtdPhi[j]->SetBinContent(i+1, integral);
4160 dPtdPhi[j]->SetBinError(i+1, integralError);
4163 // save the current state of the unfolding object
4164 SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, kFALSE);
4166 TF1* fourier = new TF1("fourier", "[0]*(1.+0.5*[1]*(TMath::Cos(2.*x)))", 0, TMath::Pi());
4167 TH1D* v2(new TH1D("v2FromFit", "v2FromFit", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
4168 for(Int_t i(0); i < ptBins; i++) {
4170 dPtdPhi[i]->Fit(fourier, "VI");
4171 Style(gPad, "PEARSON");
4172 Style(dPtdPhi[i], kBlue, kDeltaPhi);
4173 dPtdPhi[i]->DrawCopy();
4174 AliJetFlowTools::AddText(
4175 TString(Form("%.2f #LT p_{T} #LT %.2f", fBinsTrue->At(i), fBinsTrue->At(i+1))),
4182 v2->SetBinContent(1+i, fourier->GetParameter(1));
4183 v2->SetBinError(1+i, fourier->GetParError(1));
4185 v2Fits->cd(1+ptBins);
4186 Style(gPad, "PEARSON");
4187 Style(v2, kBlack, kV2, kTRUE);
4191 //_____________________________________________________________________________
4192 void AliJetFlowTools::ReplaceBins(TArrayI* array, TGraphErrors* graph) {
4194 Double_t x(0), y(0);
4195 graph->GetPoint(0, x, y);
4196 graph->SetPoint(array->At(0)-1, fBinsTrue->At(array->At(0)), y);
4197 graph->SetPointError(array->At(0)-1, 10, graph->GetErrorY(0));
4198 graph->SetPoint(array->At(1)-1, -5, -5);
4200 //_____________________________________________________________________________
4201 void AliJetFlowTools::ReplaceBins(TArrayI* array, TGraphAsymmErrors* graph) {
4203 Double_t x(0), y(0);
4204 graph->GetPoint(0, x, y);
4205 graph->SetPoint(array->At(0)-1, fBinsTrue->At(array->At(0)), y);
4206 Double_t yl = graph->GetErrorYlow(0);
4207 Double_t yh = graph->GetErrorYhigh(0);
4208 graph->SetPointError(array->At(0)-1, 10, 10, yl, yh);
4209 graph->SetPoint(array->At(1)-1, -5, -5);
4211 //_____________________________________________________________________________
4212 void AliJetFlowTools::GetSignificance(
4213 TGraphErrors* n, // points with stat error
4214 TGraphAsymmErrors* shape, // points with shape error
4215 TGraphAsymmErrors* corr, // points with stat error
4216 Int_t low, // lower bin (tgraph starts at 0)
4217 Int_t up // upper bin
4220 // calculate confidence level based on statistical uncertainty only
4221 Double_t statE(0), shapeE(0), corrE(0), totalE(0), y(0), x(0), chi2(0);
4224 for(Int_t i(low); i < up+1; i++) {
4225 n->GetPoint(i, x, y);
4226 printf(" > v2 \t %.4f \n", y);
4228 for(Int_t i(low); i < up+1; i++) {
4229 statE = n->GetErrorYlow(i);
4230 printf(" > stat \t %.4f \n", statE);
4232 for(Int_t i(low); i < up+1; i++) {
4233 shapeE = n->GetErrorYlow(i);
4234 printf(" > shape \t %.4f \n", shapeE);
4236 for(Int_t i(low); i < up+1; i++) {
4237 corrE = n->GetErrorYlow(i);
4238 printf(" > corr \t %.4f \n", corrE);
4242 // get the p value based solely on statistical uncertainties
4243 for(Int_t i(low); i < up+1; i++) {
4244 // set some flags to 0
4247 // get the nominal point
4248 n->GetPoint(i, x, y);
4249 statE = n->GetErrorYlow(i);
4250 // combine the errors
4251 chi2 += TMath::Power(y/statE, 2);
4253 cout << "p-value: p(" << chi2 << ", 6) " << TMath::Prob(chi2, 6) << endl;
4254 cout << " so the probability of finding data at least as imcompatible with 0 as the actually" << endl;
4255 cout << " observed data, using only statistical uncertainties, is " << TMath::Prob(chi2, 2) << endl << endl << endl ;
4261 // to plot the average error as function of number of events
4263 for(Int_t i(low); i < up+1; i++) {
4265 // set some flags to 0
4268 // get the nominal point
4269 n->GetPoint(i, x, y);
4270 statE += n->GetErrorY(i);
4271 shapeE += shape->GetErrorY(i);
4272 corrE += corr->GetErrorY(i);
4273 // combine the errors
4274 // totalE += TMath::Sqrt(corrE*corrE+shapeE*shapeE);// + TMath::Sqrt(corrE*corrE);
4276 cout << " AVERAGE_SHAPE " << shapeE/ctr << endl;
4277 cout << " AVERAGE_CORR " << corrE/ctr << endl;
4278 cout << " AVERAGE_STAT " << statE/ctr << endl;
4280 //_____________________________________________________________________________
4281 void AliJetFlowTools::MinimizeChi22d()
4283 // Choose method upon creation between:
4284 // kMigrad, kSimplex, kCombined,
4286 ROOT::Minuit2::Minuit2Minimizer min ( ROOT::Minuit2::kMigrad );
4287 min.SetMaxFunctionCalls(1000000);
4288 min.SetMaxIterations(100000);
4289 min.SetTolerance(0.001);
4291 ROOT::Math::Functor f(&PhenixChi22d,2);
4292 double step[] = {0.0000001, 0.0000001};
4293 double variable[] = {-1., -1.};
4296 // Set the free variables to be minimized!
4297 min.SetVariable(0,"epsilon_c",variable[0], step[0]);
4298 min.SetVariable(1,"epsilon_b",variable[1], step[1]);
4302 const double *xs = min.X();
4303 cout << endl << endl << "Minimum: f(" << xs[0] << ", " << xs[1] <<"):" << PhenixChi22d(xs) << endl;
4304 cout << "p-value: p(" << PhenixChi22d(xs) << ", 6) " << TMath::Prob(PhenixChi22d(xs), 4) << endl;
4305 cout << " so the probability of finding data at least as imcompatible with 0 as the actually" << endl;
4306 cout << " observed data is " << TMath::Prob(PhenixChi22d(xs), 4) << endl << endl << endl ;
4308 //_____________________________________________________________________________
4309 Double_t AliJetFlowTools::PhenixChi22d(const Double_t *xx )
4311 // define arrays with results and errors
4313 // these points are for 0-5 centrality, 30 - 100 gev (in which data is reported)
4333 Double_t shape[] = {
4350 // these points are for 30 - 50 centrality, 20-90 gev (in which data is reported)
4367 Double_t shape[] = {
4384 // return the function value at certain epsilon
4385 const Double_t epsc = xx[0];
4386 const Double_t epsb = xx[1];
4388 Int_t counts = (Int_t)(sizeof(v2)/sizeof(v2[0]));
4390 // implemtation of eq 3 of arXiv:0801.1665v2
4391 // this will be minimized w.r.t. 'x', which is epsilon_B in the paper
4392 for(Int_t i(1); i < counts-1; i++) {
4393 // quadratic sum of statistical and uncorrelated systematic error
4394 Double_t e = stat[i];
4396 // sum of v2 plus epsilon times correlated error minus hypothesis (0)
4397 // also the numerator of equation 3 of phenix paper
4398 Double_t numerator = TMath::Power(v2[i]+epsc*corr[i]+epsb*shape[i], 2);
4400 // denominator of equation 3 of phenix paper
4401 Double_t denominator = e*e;//TMath::Power((e*(v2[i]+epsc*corr[i]+epsb*shape[i]))/v2[i], 2);
4404 chi2 += numerator/denominator;
4406 // add the square of epsilon to the total chi2 as penalty
4407 chi2 += epsc*epsc + epsb*epsb;
4411 //_____________________________________________________________________________
4412 Double_t AliJetFlowTools::ConstructFunction2d(Double_t *x, Double_t *par)
4414 return AliJetFlowTools::PhenixChi22d(x);
4416 //_____________________________________________________________________________
4417 TF2* AliJetFlowTools::ReturnFunction2d()
4419 TF2 *f1 = new TF2("2dhist",AliJetFlowTools::ConstructFunction2d, -10, 10, -10, 10, 0);
4420 printf(" > locating minima < \n");
4421 Double_t chi2(f1->GetMinimum());
4422 f1->GetXaxis()->SetTitle("#epsilon{b}");
4423 f1->GetXaxis()->SetTitle("#epsilon_{c}");
4424 f1->GetZaxis()->SetTitle("#chi^{2}");
4426 printf(" > minimal chi2 %.8f \n", chi2);
4427 cout << " so the probability of finding data at least as imcompatible with 0 as the actually" << endl;
4428 cout << " observed data is " << TMath::Prob(chi2, 6) << endl;
4432 //_____________________________________________________________________________
4433 void AliJetFlowTools::MinimizeChi2()
4435 // Choose method upon creation between:
4436 // kMigrad, kSimplex, kCombined,
4438 ROOT::Minuit2::Minuit2Minimizer min ( ROOT::Minuit2::kMigrad );
4439 min.SetMaxFunctionCalls(1000000);
4440 min.SetMaxIterations(100000);
4441 min.SetTolerance(0.001);
4443 ROOT::Math::Functor f(&PhenixChi2,1);
4444 double step[] = {0.0000001};
4445 double variable[] = {-1.};
4448 // Set the free variables to be minimized!
4449 min.SetVariable(0,"epsilon_c",variable[0], step[0]);
4453 const double *xs = min.X();
4454 cout << endl << endl << "Minimum: f(" << xs[0] << "):" << PhenixChi2(xs) << endl;
4455 cout << "p-value: p(" << PhenixChi2(xs) << ", 6) " << TMath::Prob(PhenixChi2(xs), 6) << endl;
4456 cout << " so the probability of finding data at least as imcompatible with 0 as the actually" << endl;
4457 cout << " observed data is " << TMath::Prob(PhenixChi2(xs), 6) << endl << endl << endl;
4459 //_____________________________________________________________________________
4460 Double_t AliJetFlowTools::PhenixChi2(const Double_t *xx )
4462 // define arrays with results and errors
4481 Double_t shape[] = {
4498 // these points are for 30 - 50 centrality, 20-90 gev (in which data is reported)
4515 Double_t shape[] = {
4532 // return the function value at certain epsilon
4533 const Double_t epsc = xx[0];
4535 Int_t counts = (Int_t)(sizeof(v2)/sizeof(v2[0]));
4537 // implemtation of eq 3 of arXiv:0801.1665v2
4538 // this will be minimized w.r.t. 'x', which is epsilon_B in the paper
4539 for(Int_t i(0); i < counts; i++) {
4540 // quadratic sum of statistical and uncorrelated systematic error
4541 Double_t e = TMath::Sqrt(stat[i]*stat[i]+shape[i]*shape[i]);
4543 // sum of v2 plus epsilon times correlated error minus hypothesis (0)
4544 // also the numerator of equation 3 of phenix paper
4545 Double_t numerator = TMath::Power(v2[i]+epsc*corr[i], 2);
4547 // denominator of equation 3 of phenix paper
4548 Double_t denominator = e*e;//TMath::Power((e*(v2[i]+epsc*corr[i]))/v2[i], 2);
4551 chi2 += numerator/denominator;
4553 // add the square of epsilon to the total chi2 as penalty
4558 //_____________________________________________________________________________
4559 Double_t AliJetFlowTools::ConstructFunction(Double_t *x, Double_t *par)
4561 return AliJetFlowTools::PhenixChi2(x);
4563 //_____________________________________________________________________________
4564 TF1* AliJetFlowTools::ReturnFunction()
4566 TF1 *f1 = new TF1("1dmyfunc",AliJetFlowTools::ConstructFunction, -10, 10, 0);
4567 printf(" > locating minima < \n");
4568 Double_t chi2(f1->GetMinimum());
4569 f1->GetXaxis()->SetTitle("#epsilon_{c}");
4570 f1->GetYaxis()->SetTitle("#chi^{2}");
4573 //_____________________________________________________________________________
4574 void AliJetFlowTools::MinimizeChi2nd()
4576 // Choose method upon creation between:
4577 // kMigrad, kSimplex, kCombined,
4579 ROOT::Minuit2::Minuit2Minimizer min ( ROOT::Minuit2::kMigrad );
4580 min.SetMaxFunctionCalls(1000000);
4581 min.SetMaxIterations(100000);
4582 min.SetTolerance(0.001);
4584 ROOT::Math::Functor f(&PhenixChi2nd,2);
4585 double step[] = {0.0000001, 0.0000001};
4586 double variable[] = {-1., -1.};
4589 // Set the free variables to be minimized!
4590 min.SetVariable(0,"epsilon_c",variable[0], step[0]);
4591 min.SetVariable(1,"epsilon_b",variable[1], step[1]);
4595 const double *xs = min.X();
4596 cout << endl << endl << "Minimum: Chi2nd(" << xs[0] << ", " << xs[1] <<"):" << PhenixChi2nd(xs) << endl;
4597 cout << "p-value: p(" << PhenixChi2nd(xs) << ", 6) " << TMath::Prob(PhenixChi2nd(xs), 6) << endl;
4598 cout << " so the probability of finding data at least as imcompatible with 0 as the actually" << endl;
4599 cout << " observed data is " << TMath::Prob(PhenixChi2nd(xs), 6) << endl << endl << endl ;
4601 //_____________________________________________________________________________
4602 Double_t AliJetFlowTools::PhenixChi2nd(const Double_t *xx )
4604 // define arrays with results and errors here, see example at PhenixChi2()
4605 // very ugly, but two set of data, for 0-5 and 30-50 pct centrality
4606 // this function has to be static, so this is the easiest way to implement it in the class ...
4625 Double_t shape[] = {
4642 // these points are for 30 - 50 centrality, 20-90 gev (in which data is reported)
4659 Double_t shape[] = {
4676 // return the function value at certain epsilon
4677 const Double_t epsc = xx[0];
4678 const Double_t epsb = xx[1];
4680 Int_t counts = (Int_t)(sizeof(v2)/sizeof(v2[0]));
4682 // implemtation of eq 3 of arXiv:0801.1665v2
4683 // this will be minimized w.r.t. 'x', which is epsilon_B in the paper
4684 for(Int_t i(0); i < counts; i++) {
4685 // quadratic sum of statistical and uncorrelated systematic error
4686 Double_t e = stat[i];
4688 // sum of v2 plus epsilon times correlated error minus hypothesis (0)
4689 // also the numerator of equation 3 of phenix paper
4690 Double_t numerator = TMath::Power(v2[i]+epsc*corr[i]+epsb, 2);
4692 // denominator of equation 3 of phenix paper
4693 Double_t denominator = e*e;//TMath::Power((e*(v2[i]+epsc*corr[i]+epsb*shape[i]))/v2[i], 2);
4696 chi2 += numerator/denominator;
4698 // add the square of epsilon to the total chi2 as penalty
4700 Double_t sumEpsb(0);
4701 for(Int_t j(0); j < counts; j++) sumEpsb += (epsb*epsb)/(shape[j]*shape[j]);
4702 chi2 += epsc*epsc + sumEpsb/((float)counts);
4706 //_____________________________________________________________________________
4707 Double_t AliJetFlowTools::ConstructFunctionnd(Double_t *x, Double_t *par)
4709 return AliJetFlowTools::PhenixChi2nd(x);
4711 //_____________________________________________________________________________
4712 TF2* AliJetFlowTools::ReturnFunctionnd()
4714 TF2 *f1 = new TF2("ndhist",AliJetFlowTools::ConstructFunctionnd, -100, 100, -100, 100, 0);
4715 printf(" > locating minima < \n");
4716 Double_t chi2(f1->GetMinimum());
4717 Double_t x(0), y(0);
4718 f1->GetMinimumXY(x, y);
4719 f1->GetXaxis()->SetTitle("#epsilon{b}");
4720 f1->GetXaxis()->SetTitle("#epsilon_{c}");
4721 f1->GetZaxis()->SetTitle("#chi^{2}");
4723 printf(" > minimal chi2 f(%.8f, %.8f) = %.8f (i'm a wrong value for some reason?) \n", x, y, chi2);
4724 cout << " so the probability of finding data at least as imcompatible with 0 as the actually" << endl;
4725 cout << " observed data is " << TMath::Prob(chi2, 6) << endl;
4727 printf(" > minimal chi2 f(%.8f, %.8f) = %.8f (i should be ok ... ) \n", x, y, f1->Eval(x, y));
4728 cout << " so the probability of finding data at least as imcompatible with 0 as the actually" << endl;
4729 cout << " observed data is " << TMath::Prob(f1->Eval(x, y), 6) << endl;
4734 //_____________________________________________________________________________