AliPhysics  7baac56 (7baac56)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliJetFlowTools.cxx
Go to the documentation of this file.
1 /************************************************************************
2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // Author: Redmer Alexander Bertens, Utrecht University, Utrecht, Netherlands
17 // (rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl)
18 //
19 // Tools class for Jet Flow Analysis, replaces 'extractJetFlow.C' macro
20 //
21 // The task uses input from two analysis tasks:
22 // - $ALICE_PHYSICS/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetV2.cxx
23 // used to retrieve jet spectra and delta pt distributions
24 // - $ALICE_PHYSICS/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 ).
30 //
31 // The weak spot of this class is the function PrepareForUnfolding, which will read
32 // output from two output files and expects histograms with certain names and binning.
33 // Unfolding methods itself are general and should be able to handle any input, therefore one
34 // can forgo the PrepareForUnfolding method, and supply necessary input information via the
35 // SetRawInput() method
36 //
37 // to see an example of how to use this class, see $ALICE_PHYSICS/PWGCF/FLOW/macros/jetFlowTools.C
38 
39 // root includes
40 #include "TH1.h"
41 #include "TF2.h"
42 #include "TH2D.h"
43 #include "TGraph.h"
44 #include "TGraphErrors.h"
45 #include "TGraphAsymmErrors.h"
46 #include "TLine.h"
47 #include "TCanvas.h"
48 #include "TLegend.h"
49 #include "TArrayD.h"
50 #include "TList.h"
51 #include "TMinuit.h"
52 #include "TVirtualFitter.h"
53 #include "TLegend.h"
54 #include "TCanvas.h"
55 #include "TStyle.h"
56 #include "TLine.h"
57 #include "TVirtualFitter.h"
58 #include "TFitResultPtr.h"
59 #include "Minuit2/Minuit2Minimizer.h"
60 #include "Math/Functor.h"
61 // aliroot includes
62 #include "AliUnfolding.h"
63 #include "AliAnaChargedJetResponseMaker.h"
64 // class includes
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"
72 
73 using namespace std;
74 //_____________________________________________________________________________
76  fListPrefix (-1),
77  fResponseMaker (new AliAnaChargedJetResponseMaker()),
78  fRMS (kTRUE),
79  fSymmRMS (0),
80  fConstantUE (kTRUE),
81  fRho0 (kFALSE),
82  fBootstrap (kFALSE),
83  fPower (new TF1("fPower","[0]*TMath::Power(x,-([1]))",0.,300.)),
84  fSaveFull (kTRUE),
85  fActiveString (""),
86  fActiveDir (0x0),
87  fInputList (0x0),
88  fRefreshInput (kTRUE),
89  fOutputFileName ("UnfoldedSpectra.root"),
90  fOutputFile (0x0),
91  fCentralityArray (0x0),
92  fMergeBinsArray (0x0),
93  fCentralityWeights (0x0),
94  fMergeWithList (0x0),
95  fMergeWithCen (-1),
96  fMergeWithWeight (1.),
97  fDetectorResponse (0x0),
98  fJetFindingEff (0x0),
99  fBetaIn (.1),
100  fBetaOut (.1),
101  fBayesianIterIn (4),
102  fBayesianIterOut (4),
103  fBayesianSmoothIn (0.),
104  fBayesianSmoothOut (0.),
105  fAvoidRoundingError (kFALSE),
106  fUnfoldingAlgorithm (kChi2),
107  fPrior (kPriorMeasured),
108  fPriorUser (0x0),
109  fBinsTrue (0x0),
110  fBinsRec (0x0),
111  fBinsTruePrior (0x0),
112  fBinsRecPrior (0x0),
113  fSVDRegIn (5),
114  fSVDRegOut (5),
115  fSVDToy (kTRUE),
116  fJetRadius (0.3),
117  fEventCount (-1),
118  fNormalizeSpectra (kFALSE),
119  fSmoothenPrior (kFALSE),
120  fFitMin (60.),
121  fFitMax (300.),
122  fFitStart (75.),
123  fSmoothenCounts (kTRUE),
124  fTestMode (kFALSE),
125  fRawInputProvided (kFALSE),
126  fEventPlaneRes (.63),
127  fUseDetectorResponse(kTRUE),
128  fUseDptResponse (kTRUE),
129  fTrainPower (kTRUE),
130  fDphiUnfolding (kTRUE),
131  fDphiDptUnfolding (kFALSE),
132  fExLJDpt (kTRUE),
133  fTitleFontSize (-999.),
134  fRMSSpectrumIn (0x0),
135  fRMSSpectrumOut (0x0),
136  fRMSRatio (0x0),
137  fRMSV2 (0x0),
138  fDeltaPtDeltaPhi (0x0),
139  fJetPtDeltaPhi (0x0),
140  fSpectrumIn (0x0),
141  fSpectrumOut (0x0),
142  fDptInDist (0x0),
143  fDptOutDist (0x0),
144  fDptIn (0x0),
145  fDptOut (0x0),
146  fFullResponseIn (0x0),
147  fFullResponseOut (0x0),
148  fPivot (55.),
149  fSubdueError (kTRUE),
150  fUnfoldedSpectrumIn (0x0),
151  fUnfoldedSpectrumOut(0x0),
152  fHarmonic(2) { // class constructor
153 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
154  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
155 #endif
156  fResponseMaker->SetRMMergeWeightFunction(new TF1("weightFunction", "x*TMath::Power(1.+(1./(8.*0.9))*x, -8.)", 0, 200));
157  for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
158 }
159 //_____________________________________________________________________________
160 void AliJetFlowTools::Make(TH1* customIn, TH1* customOut) {
161  // core function of the class
162 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
163  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
164 #endif
165  if(fDphiDptUnfolding) {
166  // to extract the yield as function of Dphi, Dpt - experimental
167  MakeAU();
168  return;
169  }
170  // 1) rebin the raw output of the jet task to the desired binnings
171  // 2) calls the unfolding routine
172  // 3) writes output to file
173  // can be repeated multiple times with different configurations
174 
175  // 1) manipulation of input histograms
176  // check if the input variables are present
177  if(fRefreshInput) {
178  if(!PrepareForUnfolding(customIn, customOut)) {
179  printf(" AliJetFlowTools::Make() Fatal error \n - couldn't prepare for unfolding ! \n");
180  return;
181  }
182  }
183  // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
184  // parts of the spectrum can end up in over or underflow bins
185 
186  // if bootstrap mode is kTRUE, resample the underlying distributions
187  // FIXME think about resampling the rebinned results or raw results, could lead to difference
188  // in smoothness of tail of spectrum (which is probably not used in any case, but still ... )
189 /*
190  if(fBootstrap) {
191  // resample but leave original spectra intact for the next unfolding round
192  fSpectrumIn = reinterpret_cast<TH1D*>(Bootstrap(fSpectrumIn, kFALSE));
193  fSpectrumOut = reinterpret_cast<TH1D*>(Bootstrap(fSpectrumOut, kFALSE));
194  }
195 */
196  TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, TString("resized_in_"), kFALSE);
197  TH1D* measuredJetSpectrumOut = RebinTH1D(fSpectrumOut, fBinsRec, TString("resized_out_"), kFALSE);
198 
199  if(fBootstrap) {
200  measuredJetSpectrumIn = reinterpret_cast<TH1D*>(Bootstrap(measuredJetSpectrumIn, kFALSE));
201  measuredJetSpectrumOut = reinterpret_cast<TH1D*>(Bootstrap(measuredJetSpectrumOut, kFALSE));
202  }
203  // for now do it BEFORE as after gives an issue in Rebin function (counts are wrong)
204 
205 
206 
207  // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
208  // the template will be used as a prior for the chi2 unfolding
209  TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, TString("in"), kFALSE);
210  TH1D* measuredJetSpectrumTrueBinsOut = RebinTH1D(fSpectrumOut, fBinsTrue, TString("out"), kFALSE);
211  // get the full response matrix from the dpt and the detector response
213  // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
214  // so that unfolding should return the initial spectrum
215  if(!fTestMode) {
219  } else if (fUseDptResponse && !fUseDetectorResponse) {
222  } else if (!fUseDptResponse && fUseDetectorResponse) {
226  printf(" > No response, exiting ! < \n" );
227  return;
228  }
229  } else {
232  }
233  // normalize each slide of the response to one
236  // resize to desired binning scheme
237  TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, TString("in"));
238  TH2D* resizedResponseOut = RebinTH2D(fFullResponseOut, fBinsTrue, fBinsRec, TString("out"));
239  // get the kinematic efficiency
240  TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
241  kinematicEfficiencyIn->SetNameTitle("kin_eff_IN","kin_eff_IN");
242  TH1D* kinematicEfficiencyOut = resizedResponseOut->ProjectionX();
243  kinematicEfficiencyOut->SetNameTitle("kin_eff_OUT", "kin_eff_OUT");
244  // suppress the errors
245  for(Int_t i(0); i < kinematicEfficiencyOut->GetXaxis()->GetNbins(); i++) {
246  kinematicEfficiencyIn->SetBinError(1+i, 0.);
247  kinematicEfficiencyOut->SetBinError(1+i, 0.);
248  }
249  TH1D* jetFindingEfficiency(0x0);
250  if(fJetFindingEff) {
251  jetFindingEfficiency = ProtectHeap(fJetFindingEff);
252  jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
253  jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
254  }
255  // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
256  TH1D* unfoldedJetSpectrumIn(0x0);
257  TH1D* unfoldedJetSpectrumOut(0x0);
258  fActiveDir->cd(); // select active dir
259  TDirectoryFile* dirIn = new TDirectoryFile(Form("InPlane___%s", fActiveString.Data()), Form("InPlane___%s", fActiveString.Data()));
260  dirIn->cd(); // select inplane subdir
261  // do the inplane unfolding
262  unfoldedJetSpectrumIn = UnfoldWrapper(
263  measuredJetSpectrumIn,
264  resizedResponseIn,
265  kinematicEfficiencyIn,
266  measuredJetSpectrumTrueBinsIn,
267  TString("in"),
268  jetFindingEfficiency);
269  resizedResponseIn->SetNameTitle("ResponseMatrixIn", "response matrix in plane");
270  resizedResponseIn->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
271  resizedResponseIn->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
272  resizedResponseIn = ProtectHeap(resizedResponseIn);
273  resizedResponseIn->Write();
274  kinematicEfficiencyIn->SetNameTitle("KinematicEfficiencyIn","Kinematic efficiency, in plane");
275  kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
276  kinematicEfficiencyIn->Write();
277  fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
279  fDetectorResponse->Write();
280  // optional histograms
281  if(fSaveFull) {
282  fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", "[INPUT] Jet spectrum, in plane");
283  fSpectrumIn->Write();
284  fDptInDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, in plane");
285  fDptInDist->Write();
286  fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, in plane");
287  fDptIn->Write();
288  fFullResponseIn->SetNameTitle("ResponseMatrix", "Response matrix, in plane");
289  fFullResponseIn->Write();
290  }
291  fActiveDir->cd();
292  if(fDphiUnfolding) {
293  TDirectoryFile* dirOut = new TDirectoryFile(Form("OutOfPlane___%s", fActiveString.Data()), Form("OutOfPlane___%s", fActiveString.Data()));
294  dirOut->cd();
295  // do the out of plane unfolding
296  unfoldedJetSpectrumOut = UnfoldWrapper(
297  measuredJetSpectrumOut,
298  resizedResponseOut,
299  kinematicEfficiencyOut,
300  measuredJetSpectrumTrueBinsOut,
301  TString("out"),
302  jetFindingEfficiency);
303  resizedResponseOut->SetNameTitle("ResponseMatrixOut", "response matrix in plane");
304  resizedResponseOut->SetXTitle("p_{T, jet}^{true} (GeV/#it{c})");
305  resizedResponseOut->SetYTitle("p_{T, jet}^{rec} (GeV/#it{c})");
306  resizedResponseOut = ProtectHeap(resizedResponseOut);
307  resizedResponseOut->Write();
308  kinematicEfficiencyOut->SetNameTitle("KinematicEfficiencyOut","Kinematic efficiency, Out plane");
309  kinematicEfficiencyOut = ProtectHeap(kinematicEfficiencyOut);
310  kinematicEfficiencyOut->Write();
311  fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
313  fDetectorResponse->Write();
314  if(jetFindingEfficiency) jetFindingEfficiency->Write();
315  // optional histograms
316  if(fSaveFull) {
317  fSpectrumOut->SetNameTitle("[ORIG]JetSpectrum", "[INPUT]Jet spectrum, Out plane");
318  fSpectrumOut->Write();
319  fDptOutDist->SetNameTitle("[ORIG]DeltaPt", "#delta p_{T} distribution, Out plane");
320  fDptOutDist->Write();
321  fDptOut->SetNameTitle("[ORIG]DeltaPtMatrix","#delta p_{T} matrix, Out plane");
322  fDptOut->Write();
323  fFullResponseOut->SetNameTitle("[ORIG]ResponseMatrix", "Response matrix, Out plane");
324  fFullResponseOut->Write();
325  }
326 
327  // write general output histograms to file
328  fActiveDir->cd();
329  if(unfoldedJetSpectrumIn && unfoldedJetSpectrumOut && unfoldedJetSpectrumIn && unfoldedJetSpectrumOut) {
330  TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out")));
331  if(ratio) {
332  ratio->SetNameTitle("RatioInOutPlane", "Ratio in plane, out of plane jet spectrum");
333  ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
334  ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
335  ratio = ProtectHeap(ratio);
336  ratio->Write();
337  // write histo values to RMS files if both routines converged
338  // input values are weighted by their uncertainty
339  for(Int_t i(0); i < ratio->GetXaxis()->GetNbins(); i++) {
340  if(unfoldedJetSpectrumIn->GetBinError(i+1) > 0) fRMSSpectrumIn->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumIn->GetBinError(i+1), 2.));
341  if(unfoldedJetSpectrumOut->GetBinError(i+1) > 0) fRMSSpectrumOut->Fill(fRMSSpectrumOut->GetBinCenter(i+1), unfoldedJetSpectrumOut->GetBinContent(i+1), 1./TMath::Power(unfoldedJetSpectrumOut->GetBinError(i+1), 2.));
342  if(unfoldedJetSpectrumOut->GetBinContent(i+1) > 0) fRMSRatio->Fill(fRMSSpectrumIn->GetBinCenter(i+1), unfoldedJetSpectrumIn->GetBinContent(i+1) / unfoldedJetSpectrumOut->GetBinContent(i+1));
343  }
344  }
345  TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
346  if(v2) {
347  v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
348  v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
349  v2->GetYaxis()->SetTitle("v_{2}");
350  v2 = ProtectHeap(v2);
351  v2->Write();
352  }
353  } else if (unfoldedJetSpectrumOut && unfoldedJetSpectrumIn) {
354  TGraphErrors* ratio(GetRatio((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_in"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_out"), TString(""), kTRUE, fBinsRec->At(fBinsRec->GetSize()-1)));
355  if(ratio) {
356  ratio->SetNameTitle("[NC]RatioInOutPlane", "[NC]Ratio in plane, out of plane jet spectrum");
357  ratio->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
358  ratio->GetYaxis()->SetTitle("yield IN / yield OUT");
359  ratio = ProtectHeap(ratio);
360  ratio->Write();
361  }
362  TGraphErrors* v2(GetV2((TH1D*)unfoldedJetSpectrumIn->Clone("unfoldedLocal_inv2"), (TH1D*)unfoldedJetSpectrumOut->Clone("unfoldedLocal_outv2"), fEventPlaneRes));
363  if(v2) {
364  v2->SetNameTitle("v2", "v_{2} from different in, out of plane yield");
365  v2->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
366  v2->GetYaxis()->SetTitle("v_{2}");
367  v2 = ProtectHeap(v2);
368  v2->Write();
369  }
370  }
371  } // end of if(fDphiUnfolding)
372  fDeltaPtDeltaPhi->Write();
373  unfoldedJetSpectrumIn->Sumw2();
374  fUnfoldedSpectrumIn = ProtectHeap(unfoldedJetSpectrumIn, kFALSE);
375  unfoldedJetSpectrumIn->Write();
376  unfoldedJetSpectrumOut->Sumw2();
377  fUnfoldedSpectrumOut = ProtectHeap(unfoldedJetSpectrumOut, kFALSE);
378  unfoldedJetSpectrumOut->Write();
379  fJetPtDeltaPhi->Write();
380  // save the current state of the unfolding object
381  SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, unfoldedJetSpectrumOut ? kTRUE : kFALSE);
382  TH1D* unfoldedJetSpectrumInForSub((TH1D*)unfoldedJetSpectrumIn->Clone("forSubIn"));
383  TH1D* unfoldedJetSpectrumOutForSub((TH1D*)unfoldedJetSpectrumOut->Clone("forSubOut"));
384  unfoldedJetSpectrumInForSub->Add(unfoldedJetSpectrumOutForSub, -1.);
385  unfoldedJetSpectrumInForSub= ProtectHeap(unfoldedJetSpectrumInForSub);
386  unfoldedJetSpectrumInForSub->Write();
387 
388 }
389 //_____________________________________________________________________________
391  const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
392  const TH2D* resizedResponse, // response matrix
393  const TH1D* kinematicEfficiency, // kinematic efficiency
394  const TH1D* measuredJetSpectrumTrueBins, // unfolding template: same binning is pt gen of response
395  const TString suffix, // suffix (in or out of plane)
396  const TH1D* jetFindingEfficiency) // jet finding efficiency
397 {
398 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
399  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
400 #endif
401  // wrapper function to call specific unfolding routine
402  TH1D* (AliJetFlowTools::*myFunction)(const TH1D*, const TH2D*, const TH1D*, const TH1D*, const TString, const TH1D*);
403  // initialize functon pointer
408  else if(fUnfoldingAlgorithm == kFold) myFunction = &AliJetFlowTools::FoldSpectrum;
409  else if(fUnfoldingAlgorithm == kNone) {
410  TH1D* clone((TH1D*)measuredJetSpectrum->Clone("clone"));
411  clone->SetNameTitle(Form("MeasuredJetSpectrum%s", suffix.Data()), Form("measuredJetSpectrum %s", suffix.Data()));
412  return clone;//RebinTH1D(clone, fBinsTrue, clone->GetName(), kFALSE);
413  }
414  else return 0x0;
415  // do the actual unfolding with the selected function
416  return (this->*myFunction)( measuredJetSpectrum, resizedResponse, kinematicEfficiency, measuredJetSpectrumTrueBins, suffix, jetFindingEfficiency);
417 }
418 //_____________________________________________________________________________
420  const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
421  const TH2D* resizedResponse, // response matrix
422  const TH1D* kinematicEfficiency, // kinematic efficiency
423  const TH1D* measuredJetSpectrumTrueBins, // unfolding template: same binning is pt gen of response
424  const TString suffix, // suffix (in or out of plane)
425  const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
426 {
427 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
428  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
429 #endif
430  // unfold the spectrum using chi2 minimization
431 
432  // step 0) setup the static members of AliUnfolding
433  ResetAliUnfolding(); // reset from previous iteration
434  // also deletes and re-creates the global TVirtualFitter
435  AliUnfolding::SetUnfoldingMethod(AliUnfolding::kChi2Minimization);
436  if(!strcmp("in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
437  else if(!strcmp("out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
438  if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaIn);
439  else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetChi2Regularization(AliUnfolding::kLogLog, fBetaOut);
440  AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
441 
442  // step 1) clone all input histograms. the histograms are cloned to make sure that the original histograms
443  // stay intact. a local copy of a histogram (which only exists in the scope of this function) is
444  // denoted by the suffix 'Local'
445 
446  // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
447  TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
448  // unfolded local will be filled with the result of the unfolding
449  TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
450 
451  // full response matrix and kinematic efficiency
452  TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
453  TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
454 
455  // the initial guess for the unfolded pt spectrum, equal to the folded spectrum, but in 'true' bins
456  TH1D *priorLocal = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("priorLocal_%s", suffix.Data()));
457  // optionally, the prior can be smoothened by extrapolating the spectrum using a power law fit
458  if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
459 
460  // step 2) start the unfolding
461  Int_t status(-1), i(0);
462  while(status < 0 && i < 100) {
463  // i > 0 means that the first iteration didn't converge. in that case, the result of the first
464  // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
465  if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
466  status = AliUnfolding::Unfold(
467  resizedResponseLocal, // response matrix
468  kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
469  measuredJetSpectrumLocal, // measured spectrum
470  priorLocal, // initial conditions (set NULL to use measured spectrum)
471  unfoldedLocal); // results
472  // status holds the minuit fit status (where 0 means convergence)
473  i++;
474  }
475  // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
476  TH2D* hPearson(0x0);
477  if(status == 0 && gMinuit->fISW[1] == 3) {
478  // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
479  TVirtualFitter *fitter(TVirtualFitter::GetFitter());
480  if(gMinuit) gMinuit->Command("SET COV");
481  TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
482  TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
483  pearson->Print();
484  hPearson = new TH2D(*pearson);
485  hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
486  hPearson = ProtectHeap(hPearson);
487  hPearson->Write();
488  if(fMergeBinsArray) unfoldedLocal = MergeSpectrumBins(fMergeBinsArray, unfoldedLocal, hPearson);
489  } else status = -1;
490 
491  // step 3) refold the unfolded spectrum and save the ratio measured / refolded
492  TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
493  foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
494  unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
495  TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
496  if(ratio) {
497  ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
498  ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
499  ratio = ProtectHeap(ratio);
500  ratio->Write();
501  }
502 
503  // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
504  // objects are cloned using 'ProtectHeap()'
505  measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
506  measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
507  measuredJetSpectrumLocal->Write();
508 
509  resizedResponseLocal = ProtectHeap(resizedResponseLocal);
510  resizedResponseLocal->Write();
511 
512  unfoldedLocal = ProtectHeap(unfoldedLocal);
513  if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
514  unfoldedLocal->Write();
515 
516  foldedLocal = ProtectHeap(foldedLocal);
517  foldedLocal->Write();
518 
519  priorLocal = ProtectHeap(priorLocal);
520  priorLocal->Write();
521 
522  // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
523  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));
524  fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
525  fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
526  fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
527  fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
528  fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
529  fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
530  fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
531  fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
532  fitStatus->Write();
533 
534  return unfoldedLocal;
535 }
536 //_____________________________________________________________________________
538  const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
539  const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
540  const TH1D* kinematicEfficiency, // kinematic efficiency
541  const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
542  const TString suffix, // suffix (in, out)
543  const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
544 {
545 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
546  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
547 #endif
548  TH1D* priorLocal( GetPrior(
549  measuredJetSpectrum, // jet pt in pt rec bins
550  resizedResponse, // full response matrix, normalized in slides of pt true
551  kinematicEfficiency, // kinematic efficiency
552  measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
553  suffix, // suffix (in, out)
554  jetFindingEfficiency)); // jet finding efficiency (optional)
555  if(!priorLocal) {
556  printf(" > couldn't find prior ! < \n");
557  return 0x0;
558  } else printf(" 1) retrieved prior \n");
559 
560  // go back to the 'root' directory of this instance of the SVD unfolding routine
561  (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
562 
563  // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
564  // measured jets in pt rec binning
565  TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
566  // local copie of the response matrix
567  TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
568  // local copy of response matrix, all true slides normalized to 1
569  // this response matrix will eventually be used in the re-folding routine
570  TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
571  resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
572  // kinematic efficiency
573  TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
574  // place holder histos
575  TH1D *unfoldedLocalSVD(0x0);
576  TH1D *foldedLocalSVD(0x0);
577  cout << " 2) setup necessary input " << endl;
578  // 3) configure routine
579  RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
580  cout << " step 3) configured routine " << endl;
581 
582  // 4) get transpose matrices
583  // a) get the transpose of the full response matrix
584  TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
585  responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
586  // normalize it with the prior. this will ensure that high statistics bins will constrain the
587  // end result most strenuously than bins with limited number of counts
588  responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
589  cout << " 4) retrieved first transpose matrix " << endl;
590 
591  // 5) get response for SVD unfolding
592  RooUnfoldResponse responseSVD(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedSVD_%s", suffix.Data()), Form("respCombinedSVD_%s", suffix.Data()));
593  cout << " 5) retrieved roo unfold response object " << endl;
594 
595  // 6) actualy unfolding loop
596  RooUnfoldSvd unfoldSVD(&responseSVD, measuredJetSpectrumLocal, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
597  unfoldedLocalSVD = (TH1D*)unfoldSVD.Hreco(errorTreatment);
598  // correct the spectrum for the kinematic efficiency
599  unfoldedLocalSVD->Divide(kinematicEfficiencyLocal);
600 
601  // get the pearson coefficients from the covariance matrix
602  TMatrixD covarianceMatrix = unfoldSVD.Ereco(errorTreatment);
603  TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
604  TH2D* hPearson(0x0);
605  if(pearson) {
606  hPearson = new TH2D(*pearson);
607  pearson->Print();
608  hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
609  hPearson = ProtectHeap(hPearson);
610  hPearson->Write();
611  if(fMergeBinsArray) unfoldedLocalSVD = MergeSpectrumBins(fMergeBinsArray, unfoldedLocalSVD, hPearson);
612  } else return 0x0; // return if unfolding didn't converge
613 
614  // plot singular values and d_i vector
615  TSVDUnfold* svdUnfold(unfoldSVD.Impl());
616  TH1* hSVal(svdUnfold->GetSV());
617  TH1D* hdi(svdUnfold->GetD());
618  hSVal->SetNameTitle("SingularValuesOfAC", "Singular values of AC^{-1}");
619  hSVal->SetXTitle("singular values");
620  hSVal->Write();
621  hdi->SetNameTitle("dVector", "d vector after orthogonal transformation");
622  hdi->SetXTitle("|d_{i}^{kreg}|");
623  hdi->Write();
624  cout << " plotted singular values and d_i vector " << endl;
625 
626  // 7) refold the unfolded spectrum
627  foldedLocalSVD = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalSVD, resizedResponseLocalNorm, kinematicEfficiencyLocal);
628  TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalSVD, "ratio measured / re-folded", kTRUE));
629  ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
630  ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
631  ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
632  ratio->Write();
633  cout << " 7) refolded the unfolded spectrum " << endl;
634 
635  // write the measured, unfolded and re-folded spectra to the output directory
636  measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
637  measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
638  measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
639  measuredJetSpectrumLocal->Write(); // input spectrum
640  unfoldedLocalSVD->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
641  unfoldedLocalSVD = ProtectHeap(unfoldedLocalSVD);
642  if(jetFindingEfficiency) unfoldedLocalSVD->Divide(jetFindingEfficiency);
643  unfoldedLocalSVD->Write(); // unfolded spectrum
644  foldedLocalSVD->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
645  foldedLocalSVD = ProtectHeap(foldedLocalSVD);
646  foldedLocalSVD->Write(); // re-folded spectrum
647 
648  // save more general bookkeeeping histograms to the output directory
649  responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
650  responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
651  responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
652  responseMatrixLocalTransposePrior->Write();
653  priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
654  priorLocal->SetXTitle("p_{t} [GeV/c]");
655  priorLocal = ProtectHeap(priorLocal);
656  priorLocal->Write();
657  resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
658  resizedResponseLocalNorm->Write();
659 
660  // save some info
661  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));
662  fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fSVDRegIn : fSVDRegOut);
663  fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fSVDRegIn" : "fSVDRegOut");
664  fitStatus->Write();
665 
666  return unfoldedLocalSVD;
667 }
668 //_____________________________________________________________________________
670  const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
671  const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
672  const TH1D* kinematicEfficiency, // kinematic efficiency
673  const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
674  const TString suffix, // suffix (in, out)
675  const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
676 {
677 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
678  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
679 #endif
680  // unfold the spectrum using the bayesian unfolding impelmented in AliUnfolding
681  // FIXME careful, not tested yet ! (06122013) FIXME
682 
683  // step 0) setup the static members of AliUnfolding
684  ResetAliUnfolding(); // reset from previous iteration
685  // also deletes and re-creates the global TVirtualFitter
686  AliUnfolding::SetUnfoldingMethod(AliUnfolding::kBayesian);
687  if(!strcmp("in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
688  else if(!strcmp("out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
689  else if(!strcmp("prior_in", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothIn, fBayesianIterIn);
690  else if(!strcmp("prior_out", suffix.Data())) AliUnfolding::SetBayesianParameters(fBayesianSmoothOut, fBayesianIterOut);
691  AliUnfolding::SetNbins(fBinsRec->GetSize()-1, fBinsTrue->GetSize()-1);
692 
693  // 1) get a prior for unfolding and clone all the input histograms
694  TH1D* priorLocal( GetPrior(
695  measuredJetSpectrum, // jet pt in pt rec bins
696  resizedResponse, // full response matrix, normalized in slides of pt true
697  kinematicEfficiency, // kinematic efficiency
698  measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
699  suffix, // suffix (in, out)
700  jetFindingEfficiency)); // jet finding efficiency (optional)
701  if(!priorLocal) {
702  printf(" > couldn't find prior ! < \n");
703  return 0x0;
704  } else printf(" 1) retrieved prior \n");
705  // switch back to root dir of this unfolding procedure
706  (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
707 
708  // measuredJetSpectrumLocal holds the spectrum that needs to be unfolded
709  TH1D *measuredJetSpectrumLocal = (TH1D*)measuredJetSpectrum->Clone(Form("measuredJetSpectrumLocal_%s", suffix.Data()));
710  // unfolded local will be filled with the result of the unfolding
711  TH1D *unfoldedLocal(new TH1D(Form("unfoldedLocal_%s", suffix.Data()), Form("unfoldedLocal_%s", suffix.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
712 
713  // full response matrix and kinematic efficiency
714  TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
715  TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
716 
717  // step 2) start the unfolding
718  Int_t status(-1), i(0);
719  while(status < 0 && i < 100) {
720  // i > 0 means that the first iteration didn't converge. in that case, the result of the first
721  // iteration (stored in unfoldedLocal) is cloned and used as a starting point for the
722  if (i > 0) priorLocal = (TH1D*)unfoldedLocal->Clone(Form("priorLocal_%s_%i", suffix.Data(), i));
723  status = AliUnfolding::Unfold(
724  resizedResponseLocal, // response matrix
725  kinematicEfficiencyLocal, // efficiency applied on the unfolded spectrum (can be NULL)
726  measuredJetSpectrumLocal, // measured spectrum
727  priorLocal, // initial conditions (set NULL to use measured spectrum)
728  unfoldedLocal); // results
729  // status holds the minuit fit status (where 0 means convergence)
730  i++;
731  }
732  // get the status of TMinuit::mnhess(), fISW[1] == 3 means the hessian matrix was calculated succesfully
733  TH2D* hPearson(0x0);
734  if(status == 0 && gMinuit->fISW[1] == 3) {
735  // if the unfolding converged and the hessian matrix is reliable, plot the pearson coefficients
736  TVirtualFitter *fitter(TVirtualFitter::GetFitter());
737  if(gMinuit) gMinuit->Command("SET COV");
738  TMatrixD covarianceMatrix(fBinsTrue->GetSize()-1, fBinsTrue->GetSize()-1, fitter->GetCovarianceMatrix());
739  TMatrixD *pearson((TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix));
740  pearson->Print();
741  hPearson= new TH2D(*pearson);
742  hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients, %s plane", suffix.Data()));
743  hPearson = ProtectHeap(hPearson);
744  hPearson->Write();
745  } else status = -1;
746  if(fMergeBinsArray) unfoldedLocal = MergeSpectrumBins(fMergeBinsArray, unfoldedLocal, hPearson);
747 
748  // step 3) refold the unfolded spectrum and save the ratio measured / refolded
749  TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
750  foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
751  unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
752  TGraphErrors* ratio(GetRatio(foldedLocal, measuredJetSpectrumLocal, TString(""), kTRUE, fBinsTrue->At(fBinsTrue->GetSize()-1)));
753  if(ratio) {
754  ratio->SetNameTitle("RatioRefoldedMeasured", Form("Ratio measured, re-folded %s ", suffix.Data()));
755  ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
756  ratio = ProtectHeap(ratio);
757  ratio->Write();
758  }
759 
760  // step 4) write histograms to file. to ensure that these have unique identifiers on the heap,
761  // objects are cloned using 'ProtectHeap()'
762  measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
763  measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
764  measuredJetSpectrumLocal->Write();
765 
766  resizedResponseLocal = ProtectHeap(resizedResponseLocal);
767  resizedResponseLocal->Write();
768 
769  unfoldedLocal = ProtectHeap(unfoldedLocal);
770  if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
771  unfoldedLocal->Write();
772 
773  foldedLocal = ProtectHeap(foldedLocal);
774  foldedLocal->Write();
775 
776  priorLocal = ProtectHeap(priorLocal);
777  priorLocal->Write();
778 
779  // step 5) save the fit status (penalty value, degrees of freedom, chi^2 value)
780  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));
781  fitStatus->SetBinContent(1, AliUnfolding::fChi2FromFit);
782  fitStatus->GetXaxis()->SetBinLabel(1, "fChi2FromFit");
783  fitStatus->SetBinContent(2, AliUnfolding::fPenaltyVal);
784  fitStatus->GetXaxis()->SetBinLabel(2, "fPenaltyVal");
785  fitStatus->SetBinContent(3, fBinsRec->GetSize()-fBinsTrue->GetSize());
786  fitStatus->GetXaxis()->SetBinLabel(3, "DOF");
787  fitStatus->SetBinContent(4, (!strcmp(suffix.Data(), "in")) ? fBetaIn : fBetaOut);
788  fitStatus->GetXaxis()->SetBinLabel(4, (!strcmp(suffix.Data(), "in")) ? "fBetaIn" : "fBetaOut");
789  fitStatus->Write();
790 
791  return unfoldedLocal;
792 }
793 //_____________________________________________________________________________
795  const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
796  const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
797  const TH1D* kinematicEfficiency, // kinematic efficiency
798  const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
799  const TString suffix, // suffix (in, out)
800  const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
801 {
802 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
803  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
804 #endif
805  // use bayesian unfolding from the RooUnfold package to unfold jet spectra
806 
807  // 1) get a prior for unfolding.
808  TH1D* priorLocal( GetPrior(
809  measuredJetSpectrum, // jet pt in pt rec bins
810  resizedResponse, // full response matrix, normalized in slides of pt true
811  kinematicEfficiency, // kinematic efficiency
812  measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
813  suffix, // suffix (in, out)
814  jetFindingEfficiency)); // jet finding efficiency (optional)
815  if(!priorLocal) {
816  printf(" > couldn't find prior ! < \n");
817  return 0x0;
818  } else printf(" 1) retrieved prior \n");
819  (!strcmp(suffix.Data(), "in")) ? fActiveDir->cd(Form("InPlane___%s", fActiveString.Data())) : fActiveDir->cd(Form("OutOfPlane___%s", fActiveString.Data()));
820 
821  // 2) setup all the necessary input for the unfolding routine. all input histograms are copied locally
822  // measured jets in pt rec binning
823  TH1D *measuredJetSpectrumLocal((TH1D*)measuredJetSpectrum->Clone(Form("jets_%s", suffix.Data())));
824  // local copie of the response matrix
825  TH2D *resizedResponseLocal((TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data())));
826  // local copy of response matrix, all true slides normalized to 1
827  // this response matrix will eventually be used in the re-folding routine
828  TH2D *resizedResponseLocalNorm((TH2D*)resizedResponse->Clone(Form("resizedResponseLocalNorm_%s", suffix.Data())));
829  resizedResponseLocalNorm = NormalizeTH2D(resizedResponseLocalNorm);
830  // kinematic efficiency
831  TH1D *kinematicEfficiencyLocal((TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiency_%s", suffix.Data())));
832  // place holder histos
833  TH1D *unfoldedLocalBayes(0x0);
834  TH1D *foldedLocalBayes(0x0);
835  cout << " 2) setup necessary input " << endl;
836  // 4) get transpose matrices
837  // a) get the transpose of the full response matrix
838  TH2* responseMatrixLocalTransposePrior(fResponseMaker->GetTransposeResponsMatrix(resizedResponseLocal));
839  responseMatrixLocalTransposePrior->SetNameTitle(Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()),Form("prior_%s_%s", responseMatrixLocalTransposePrior->GetName(), suffix.Data()));
840  // normalize it with the prior. this will ensure that high statistics bins will constrain the
841  // end result most strenuously than bins with limited number of counts
842  responseMatrixLocalTransposePrior = fResponseMaker->NormalizeResponsMatrixYaxisWithPrior(responseMatrixLocalTransposePrior, priorLocal);
843  // 3) get response for Bayesian unfolding
844  RooUnfoldResponse responseBayes(0, 0, responseMatrixLocalTransposePrior, Form("respCombinedBayes_%s", suffix.Data()), Form("respCombinedBayes_%s", suffix.Data()));
845 
846  // 4) actualy unfolding loop
847  RooUnfoldBayes unfoldBayes(&responseBayes, measuredJetSpectrumLocal, (!strcmp("in", suffix.Data())) ? fBayesianIterIn : fBayesianIterOut);
848  RooUnfold::ErrorTreatment errorTreatment = (fSVDToy) ? RooUnfold::kCovToy : RooUnfold::kCovariance;
849  unfoldedLocalBayes = (TH1D*)unfoldBayes.Hreco(errorTreatment);
850  // correct the spectrum for the kinematic efficiency
851  unfoldedLocalBayes->Divide(kinematicEfficiencyLocal);
852  // get the pearson coefficients from the covariance matrix
853  TMatrixD covarianceMatrix = unfoldBayes.Ereco(errorTreatment);
854  TMatrixD *pearson = (TMatrixD*)CalculatePearsonCoefficients(&covarianceMatrix);
855  TH2D* hPearson(0x0);
856  if(pearson) {
857  hPearson = new TH2D(*pearson);
858  pearson->Print();
859  hPearson->SetNameTitle(Form("PearsonCoefficients_%s", suffix.Data()), Form("Pearson coefficients_%s", suffix.Data()));
860  hPearson = ProtectHeap(hPearson);
861  hPearson->Write();
862  if(fMergeBinsArray) unfoldedLocalBayes = MergeSpectrumBins(fMergeBinsArray, unfoldedLocalBayes, hPearson);
863  } else return 0x0; // return if unfolding didn't converge
864 
865  // 5) refold the unfolded spectrum
866  foldedLocalBayes = fResponseMaker->MultiplyResponseGenerated(unfoldedLocalBayes, resizedResponseLocalNorm, kinematicEfficiencyLocal);
867  TGraphErrors* ratio(GetRatio(measuredJetSpectrumLocal, foldedLocalBayes, "ratio measured / re-folded", kTRUE));
868  ratio->SetNameTitle(Form("RatioRefoldedMeasured_%s", fActiveString.Data()), Form("Ratio measured / re-folded %s", fActiveString.Data()));
869  ratio->GetXaxis()->SetTitle("p_{t}^{rec, rec} [GeV/ c]");
870  ratio->GetYaxis()->SetTitle("ratio measured / re-folded");
871  ratio->Write();
872  cout << " 7) refolded the unfolded spectrum " << endl;
873 
874  // write the measured, unfolded and re-folded spectra to the output directory
875  measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("input spectrum (measured) %s", suffix.Data()));
876  measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
877  measuredJetSpectrumLocal->SetXTitle("p_{t}^{rec} [GeV/c]");
878  measuredJetSpectrumLocal->Write(); // input spectrum
879  unfoldedLocalBayes->SetNameTitle(Form("UnfoldedSpectrum_%s",suffix.Data()), Form("unfolded spectrum %s", suffix.Data()));
880  unfoldedLocalBayes = ProtectHeap(unfoldedLocalBayes);
881  if(jetFindingEfficiency) unfoldedLocalBayes->Divide(jetFindingEfficiency);
882  unfoldedLocalBayes->Write(); // unfolded spectrum
883  foldedLocalBayes->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("refoldedSpectrum_%s", suffix.Data()));
884  foldedLocalBayes = ProtectHeap(foldedLocalBayes);
885  foldedLocalBayes->Write(); // re-folded spectrum
886 
887  // save more general bookkeeeping histograms to the output directory
888  responseMatrixLocalTransposePrior->SetNameTitle("TransposeResponseMatrix", "Transpose of response matrix, normalize with prior");
889  responseMatrixLocalTransposePrior->SetXTitle("p_{T, jet}^{true} [GeV/c]");
890  responseMatrixLocalTransposePrior->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
891  responseMatrixLocalTransposePrior->Write();
892  priorLocal->SetNameTitle("PriorOriginal", "Prior, original");
893  priorLocal->SetXTitle("p_{t} [GeV/c]");
894  priorLocal = ProtectHeap(priorLocal);
895  priorLocal->Write();
896  resizedResponseLocalNorm = ProtectHeap(resizedResponseLocalNorm);
897  resizedResponseLocalNorm->Write();
898 
899  // save some info
900  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));
901  fitStatus->SetBinContent(1, (!strcmp(suffix.Data(), "in")) ? fBayesianIterIn : fBayesianIterOut);
902  fitStatus->GetXaxis()->SetBinLabel(1, (!strcmp(suffix.Data(), "in")) ? "fBayesianIterIn" : "fBayesianIterOut");
903  fitStatus->Write();
904 
905  return unfoldedLocalBayes;
906 }
907 //_____________________________________________________________________________
909  const TH1D* measuredJetSpectrum, // truncated raw jets (same binning as pt rec of response)
910  const TH2D* resizedResponse, // response matrix
911  const TH1D* kinematicEfficiency, // kinematic efficiency
912  const TH1D* measuredJetSpectrumTrueBins, // not used
913  const TString suffix, // suffix (in or out of plane)
914  const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
915 {
916 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
917  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
918 #endif
919  // simple function to fold the given spectrum with the in-plane and out-of-plane response
920 
921  if(!measuredJetSpectrumTrueBins) SquelchWarning();
922  // 0) for consistency with the other methods, keep the same nomenclature which admittedly is a bit confusing
923  // what is 'unfolded' here, is just a clone of the input spectrum, binned to the 'unfolded' binning
924  TH1D* unfoldedLocal((TH1D*)measuredJetSpectrum->Clone(Form("unfoldedLocal_%s", suffix.Data())));
925 
926  // 1) full response matrix and kinematic efficiency
927  TH2D* resizedResponseLocal = (TH2D*)resizedResponse->Clone(Form("resizedResponseLocal_%s", suffix.Data()));
928  TH1D* kinematicEfficiencyLocal = (TH1D*)kinematicEfficiency->Clone(Form("kinematicEfficiencyLocal_%s", suffix.Data()));
929 
930  // step 2) fold the 'unfolded' spectrum and save the ratio measured / refolded
931  TH1D *foldedLocal(fResponseMaker->MultiplyResponseGenerated(unfoldedLocal, resizedResponseLocal,kinematicEfficiencyLocal));
932  foldedLocal->SetNameTitle(Form("RefoldedSpectrum_%s", suffix.Data()), Form("Refolded jet spectrum, %s plane", suffix.Data()));
933  unfoldedLocal->SetNameTitle(Form("UnfoldedSpectrum_%s", suffix.Data()), Form("Unfolded jet spectrum, %s plane", suffix.Data()));
934 
935  // step 3) write histograms to file. to ensure that these have unique identifiers on the heap,
936  // objects are cloned using 'ProtectHeap()'
937  TH1D* measuredJetSpectrumLocal((TH1D*)(measuredJetSpectrum->Clone("tempObject")));
938  measuredJetSpectrumLocal->SetNameTitle(Form("InputSpectrum_%s", suffix.Data()), Form("InputSpectrum_%s", suffix.Data()));
939  measuredJetSpectrumLocal = ProtectHeap(measuredJetSpectrumLocal);
940  measuredJetSpectrumLocal->Write();
941 
942  resizedResponseLocal = ProtectHeap(resizedResponseLocal);
943  resizedResponseLocal->Write();
944 
945  unfoldedLocal = ProtectHeap(unfoldedLocal);
946  if(jetFindingEfficiency) unfoldedLocal->Divide(jetFindingEfficiency);
947  unfoldedLocal->Write();
948 
949  foldedLocal = ProtectHeap(foldedLocal);
950  foldedLocal->Write();
951 
952  // return the folded result
953  return foldedLocal;
954 }
955 //_____________________________________________________________________________
957 {
958 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
959  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
960 #endif
961  // prepare for unfolding. check if all the necessary input data is available, if not, retrieve it
962  if(fRawInputProvided) return kTRUE;
963  if(!fInputList) {
964  printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
965  return kFALSE;
966  }
967  if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
968  // check if the pt bin for true and rec have been set
969  if(!fBinsTrue || !fBinsRec) {
970  printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
971  return kFALSE;
972  }
973  if(!fRMSSpectrumIn && fDphiUnfolding) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
974  // procedures, these profiles will be nonsensical, user is responsible
975  fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
976  fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
977  fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
978  }
979  if(!fTrainPower) {
980  // clear minuit state to avoid constraining the fit with the results of the previous iteration
981  for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
982  }
983  // extract the spectra
984  TString spectrumName(Form("fHistJetPsi3Pt_%i", fCentralityArray->At(0)));
985  if(fRho0) spectrumName = Form("fHistJetPsi3PtRho0_%i", fCentralityArray->At(0));
986  if(!fInputList->FindObject(spectrumName.Data())) {
987  printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
988  return kFALSE;
989  }
990 
991  // get the first scaled spectrum
992  fJetPtDeltaPhi = (TH2D*)fInputList->FindObject(spectrumName.Data());
993  // clone the spectrum on the heap. this is necessary since scale or add change the
994  // contents of the original histogram
996  fJetPtDeltaPhi->Scale(fCentralityWeights->At(0));
997  printf("Extracted %s wight weight %.2f \n", spectrumName.Data(), fCentralityWeights->At(0));
998  // merge subsequent bins (if any)
999  for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1000  spectrumName = Form("fHistJetPsi3Pt_%i", fCentralityArray->At(i));
1001  printf( " Merging with %s with weight %.4f \n", spectrumName.Data(), fCentralityWeights->At(i));
1002  fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())), fCentralityWeights->At(i));
1003  }
1004  // if a second list is passed, merge this with the existing one
1005  if(fMergeWithList) {
1006  spectrumName = Form("fHistJetPsi3Pt_%i", fMergeWithCen);
1007  printf( " Adding additional output %s \n", spectrumName.Data());
1008  fJetPtDeltaPhi->Add(((TH2D*)fMergeWithList->FindObject(spectrumName.Data())), fMergeWithWeight);
1009  }
1010 
1011  // in plane spectrum
1012  if(!fDphiUnfolding) {
1013  fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), 1, 40, "e");
1014  fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 1, 40, "e");
1015  } else {
1016  fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_ina_%s", spectrumName.Data()), 1, 10, "e");
1017  fSpectrumIn->Add(fJetPtDeltaPhi->ProjectionY(Form("_py_inb_%s", spectrumName.Data()), 31, 40, "e"));
1019  // out of plane spectrum
1020  fSpectrumOut = fJetPtDeltaPhi->ProjectionY(Form("_py_out_%s", spectrumName.Data()), 11, 30, "e");
1022  }
1023  // if a custom input is passed, overwrite existing histograms
1024  if(customIn) fSpectrumIn = dynamic_cast<TH1D*>(customIn);
1025  if(customOut) fSpectrumOut = dynamic_cast<TH1D*>(customOut);
1026 
1027  // normalize spectra to event count if requested
1028  if(fNormalizeSpectra) {
1029  TH1* rho((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(0))));
1030  rho->Scale(fCentralityWeights->At(0));
1031  for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1032  rho->Add((TH1*)fInputList->FindObject(Form("fHistRho_%i", fCentralityArray->At(i))), fCentralityWeights->At(i));
1033  }
1034  if(!rho) return 0x0;
1035  Bool_t normalizeToFullSpectrum = (fEventCount < 0) ? kTRUE : kFALSE;
1036  if (normalizeToFullSpectrum) fEventCount = rho->GetEntries();
1037  if(fEventCount > 0) {
1038  fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
1039  fSpectrumOut->Sumw2();
1040  Double_t pt(0);
1041  Double_t error(0); // lots of issues with the errors here ...
1042  for(Int_t i(0); i < fSpectrumIn->GetXaxis()->GetNbins(); i++) {
1043  pt = fSpectrumIn->GetBinContent(1+i)/fEventCount; // normalized count
1044  error = 1./((double)(fEventCount*fEventCount))*fSpectrumIn->GetBinError(1+i)*fSpectrumIn->GetBinError(1+i);
1045  fSpectrumIn->SetBinContent(1+i, pt);
1046  if(pt <= 0 ) fSpectrumIn->SetBinError(1+i, 0.);
1047  if(error > 0) fSpectrumIn->SetBinError(1+i, error);
1048  else fSpectrumIn->SetBinError(1+i, TMath::Sqrt(pt));
1049  }
1050  for(Int_t i(0); i < fSpectrumOut->GetXaxis()->GetNbins(); i++) {
1051  pt = fSpectrumOut->GetBinContent(1+i)/fEventCount; // normalized count
1052  error = 1./((double)(fEventCount*fEventCount))*fSpectrumOut->GetBinError(1+i)*fSpectrumOut->GetBinError(1+i);
1053  fSpectrumOut->SetBinContent(1+i, pt);
1054  if( pt <= 0) fSpectrumOut->SetBinError(1+i, 0.);
1055  if(error > 0) fSpectrumOut->SetBinError(1+i, error);
1056  else fSpectrumOut->SetBinError(1+i, TMath::Sqrt(pt));
1057  }
1058  }
1059  if(normalizeToFullSpectrum) fEventCount = -1;
1060  }
1061  // extract the delta pt matrices
1062  TString deltaptName("");
1063  if(!fRho0) {
1064  deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi3ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi3_%i", fCentralityArray->At(0));
1065  } else {
1066  deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi3ExLJRho0_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi3Rho0_%i", fCentralityArray->At(0));
1067  }
1068  fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1069  if(!fDeltaPtDeltaPhi) {
1070  printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
1071  printf(" > may be ok, depending no what you want to do < \n");
1072  fRefreshInput = kTRUE;
1073  return kTRUE;
1074  }
1075 
1076  // clone the distribution on the heap and if requested merge with other centralities
1078  fDeltaPtDeltaPhi->Scale(fCentralityWeights->At(0));
1079  printf("Extracted %s with weight %.2f \n", deltaptName.Data(), fCentralityWeights->At(0));
1080  for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1081  deltaptName = (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi3ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi3_%i", fCentralityArray->At(i));
1082  printf(" Merging with %s with weight %.4f \n", deltaptName.Data(), fCentralityWeights->At(i));
1083  fDeltaPtDeltaPhi->Add((TH2D*)fInputList->FindObject(deltaptName.Data()), fCentralityWeights->At(i));
1084  }
1085  // if a second list is passed, merge this with the existing one
1086  if(fMergeWithList) {
1087  deltaptName = (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi3ExLJ_%i", fMergeWithCen) : Form("fHistDeltaPtDeltaPhi3_%i", fMergeWithCen);
1088  printf(" Adding additional data %s \n", deltaptName.Data());
1089  fDeltaPtDeltaPhi->Add((TH2D*)fMergeWithList->FindObject(deltaptName.Data()), fMergeWithWeight);
1090  }
1091 
1092  // in plane delta pt distribution
1093  if(!fDphiUnfolding) {
1094  fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), 1, 40, "e");
1095  fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 1, 40, "e");
1096  } else {
1097  fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_ina_%s", deltaptName.Data()), 1, 10, "e");
1098  fDptInDist->Add(fDeltaPtDeltaPhi->ProjectionY(Form("_py_inb_%s", deltaptName.Data()), 31, 40, "e"));
1099  // out of plane delta pt distribution
1100  fDptOutDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_out_%s", deltaptName.Data()), 11, 30, "e");
1103  // TODO get dpt response matrix from ConstructDPtResponseFromTH1D
1104  }
1105 
1106  // create a rec - true smeared response matrix
1107  TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1108  for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
1109  Bool_t skip = kFALSE;
1110  for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
1111  (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1112  if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1113  }
1114  }
1115  TMatrixD* rfOut = new TMatrixD(-50, 249, -50, 249);
1116  for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
1117  Bool_t skip = kFALSE;
1118  for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
1119  (*rfOut)(k, j) = (skip) ? 0. : fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j));
1120  if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptOutDist->GetBinContent(fDptOutDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1121  }
1122  }
1123  fDptIn = new TH2D(*rfIn);
1124  fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
1125  fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1126  fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1128  fDptOut = new TH2D(*rfOut);
1129  fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
1130  fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1131  fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1133 
1134  fRefreshInput = kTRUE; // force cloning of the input
1135  return kTRUE;
1136 }
1137 //_____________________________________________________________________________
1139  // prepare for unfoldingUA - more robust method to extract input spectra from file
1140  // will replace PrepareForUnfolding eventually (09012014)
1141 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1142  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1143 #endif
1144  if(!fInputList) {
1145  printf(" AliJetFlowTools::PrepareForUnfolding() fInputList not found \n - Set a list using AliJetFlowTools::SetInputList() \n");
1146  return kFALSE;
1147  }
1148  if(!fDetectorResponse) printf(" WARNING, no detector response supplied ! May be ok (depending on what you want to do) \n ");
1149  // check if the pt bin for true and rec have been set
1150  if(!fBinsTrue || !fBinsRec) {
1151  printf(" AliJetFlowTools::PrepareForUnfolding() no true or rec bins set, aborting ! \n");
1152  return kFALSE;
1153  }
1154  if(!fTrainPower) {
1155  // clear minuit state to avoid constraining the fit with the results of the previous iteration
1156  for(Int_t i(0); i < fPower->GetNpar(); i++) fPower->SetParameter(i, 0.);
1157  }
1158  // extract the spectra
1159  TString spectrumName(Form("fHistJetPsi3Pt_%i", fCentralityArray->At(0)));
1160  if(fRho0) spectrumName = Form("fHistJetPsi3PtRho0_%i", fCentralityArray->At(0));
1161  fJetPtDeltaPhi = ((TH2D*)fInputList->FindObject(spectrumName.Data()));
1162  if(!fJetPtDeltaPhi) {
1163  printf(" Couldn't find spectrum %s ! \n", spectrumName.Data());
1164  return kFALSE;
1165  }
1166  if(fCentralityArray) {
1167  for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1168  spectrumName = Form("fHistJetPsi3Pt_%i", fCentralityArray->At(i));
1169  fJetPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(spectrumName.Data())));
1170  }
1171  }
1172  fJetPtDeltaPhi = ProtectHeap(fJetPtDeltaPhi, kFALSE);
1173  // in plane spectrum
1174  fSpectrumIn = fJetPtDeltaPhi->ProjectionY(Form("_py_in_%s", spectrumName.Data()), low, up, "e");
1175  // extract the delta pt matrices
1176  TString deltaptName("");
1177  if(!fRho0) {
1178  deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi3ExLJ_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi3_%i", fCentralityArray->At(0));
1179  } else {
1180  deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi3ExLJRho0_%i", fCentralityArray->At(0)) : Form("fHistDeltaPtDeltaPhi3Rho0_%i", fCentralityArray->At(0));
1181  }
1182  fDeltaPtDeltaPhi = ((TH2D*)fInputList->FindObject(deltaptName.Data()));
1183  if(!fDeltaPtDeltaPhi) {
1184  printf(" Couldn't find delta pt matrix %s ! \n", deltaptName.Data());
1185  printf(" > may be ok, depending no what you want to do < \n");
1186  fRefreshInput = kTRUE;
1187  return kTRUE;
1188  }
1189  if(fCentralityArray) {
1190  for(Int_t i(1); i < fCentralityArray->GetSize(); i++) {
1191  deltaptName += (fExLJDpt) ? Form("fHistDeltaPtDeltaPhi3ExLJ_%i", fCentralityArray->At(i)) : Form("fHistDeltaPtDeltaPhi3_%i", fCentralityArray->At(i));
1192  fDeltaPtDeltaPhi->Add(((TH2D*)fInputList->FindObject(deltaptName.Data())));
1193  }
1194  }
1195 
1197  // in plane delta pt distribution
1198  fDptInDist = fDeltaPtDeltaPhi->ProjectionY(Form("_py_in_%s", deltaptName.Data()), low, up, "e");
1199  // create a rec - true smeared response matrix
1200  TMatrixD* rfIn = new TMatrixD(-50, 249, -50, 249);
1201  for(Int_t j(-50); j < 250; j++) { // loop on pt true slices j
1202  Bool_t skip = kFALSE;
1203  for(Int_t k(-50); k < 250; k++) { // loop on pt gen slices k
1204  (*rfIn)(k, j) = (skip) ? 0. : fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j));
1205  if(fAvoidRoundingError && k > j && TMath::AreEqualAbs(fDptInDist->GetBinContent(fDptInDist->GetXaxis()->FindBin(k-j)), 0, 1e-8)) skip = kTRUE;
1206  }
1207  }
1208  fDptIn = new TH2D(*rfIn);
1209  fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
1210  fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
1211  fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
1213 
1214  return kTRUE;
1215 }
1216 //_____________________________________________________________________________
1218  const TH1D* measuredJetSpectrum, // jet pt in pt rec bins
1219  const TH2D* resizedResponse, // full response matrix, normalized in slides of pt true
1220  const TH1D* kinematicEfficiency, // kinematic efficiency
1221  const TH1D* measuredJetSpectrumTrueBins, // jet pt in pt true bins, also the prior when measured is chosen as prior
1222  const TString suffix, // suffix (in, out)
1223  const TH1D* jetFindingEfficiency) // jet finding efficiency (optional)
1224 {
1225 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1226  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1227 #endif
1228  // 1) get a prior for unfolding.
1229  // this can be either an unfolded spectrum from e.g. chi2 unfolding or the measured spectrum
1230  TH1D* unfolded(0x0);
1231  TDirectoryFile* dirOut = new TDirectoryFile(Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()), Form("Prior_%s___%s", suffix.Data(), fActiveString.Data()));
1232  dirOut->cd();
1233  switch (fPrior) { // select the prior for unfolding
1234  case kPriorTF1 : {
1235  if(!fPriorUser) {
1236  printf("> GetPrior:: FATAL ERROR! TF1 prior requested but prior has not been set ! < \n");
1237  return 0x0;
1238  }
1239  return fPriorUser;
1240  } break;
1241  case kPriorPythia : {
1242  if(!fPriorUser) {
1243  printf("> GetPrior:: FATAL ERROR! pythia prior requested but prior has not been set ! < \n");
1244  return 0x0;
1245  }
1246  // rebin the given prior to the true spectrum (creates a new histo)
1247  return RebinTH1D(fPriorUser, fBinsTrue, Form("kPriorPythia_%s", suffix.Data()), kFALSE);
1248  } break;
1249  case kPriorChi2 : {
1250  TArrayI* placeHolder(0x0);
1251  if(fMergeBinsArray) {
1252  placeHolder = fMergeBinsArray;
1253  fMergeBinsArray = 0x0;
1254  }
1255  if(fBinsTruePrior && fBinsRecPrior) { // if set, use different binning for the prior
1256  TArrayD* tempArrayTrue(fBinsTrue); // temporarily cache the original binning
1257  fBinsTrue = fBinsTruePrior; // switch binning schemes (will be used in UnfoldSpectrumChi2())
1258  TArrayD* tempArrayRec(fBinsRec);
1260  // for the prior, do not re-bin the output
1261  TH1D* measuredJetSpectrumChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsRec, TString("resized_chi2"), kFALSE);
1262  TH1D* measuredJetSpectrumTrueBinsChi2 = RebinTH1D((!strcmp("in", suffix.Data())) ? fSpectrumIn : fSpectrumOut, fBinsTruePrior, TString("out"), kFALSE);
1263  TH2D* resizedResponseChi2(RebinTH2D((!strcmp("in", suffix.Data())) ? fFullResponseIn : fFullResponseOut,fBinsTruePrior, fBinsRec, TString("chi2")));
1264  TH1D* kinematicEfficiencyChi2(resizedResponseChi2->ProjectionX());
1265  kinematicEfficiencyChi2->SetNameTitle("kin_eff_chi2","kin_eff_chi2");
1266  for(Int_t i(0); i < kinematicEfficiencyChi2->GetXaxis()->GetNbins(); i++) kinematicEfficiencyChi2->SetBinError(1+i, 0.);
1267  unfolded= UnfoldSpectrumChi2(
1268  measuredJetSpectrumChi2,
1269  resizedResponseChi2,
1270  kinematicEfficiencyChi2,
1271  measuredJetSpectrumTrueBinsChi2, // prior for chi2 unfolding (measured)
1272  TString(Form("prior_%s", suffix.Data())));
1273  if(!unfolded) {
1274  printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1275  printf(" probably Chi2 unfolding did not converge < \n");
1276  if(placeHolder) fMergeBinsArray = placeHolder;
1277  return 0x0;
1278  }
1279  fBinsTrue = tempArrayTrue; // reset bins borders
1280  fBinsRec = tempArrayRec;
1281  // if the chi2 prior has a different binning, rebin to the true binning for the unfolding
1282  unfolded = RebinTH1D(unfolded, fBinsTrue, TString(Form("unfoldedChi2Prior_%s", suffix.Data()))); // rebin unfolded
1283  } else {
1284  unfolded = UnfoldSpectrumChi2(
1285  measuredJetSpectrum,
1286  resizedResponse,
1287  kinematicEfficiency,
1288  measuredJetSpectrumTrueBins, // prior for chi2 unfolding (measured)
1289  TString(Form("prior_%s", suffix.Data())));
1290  if(!unfolded) {
1291  printf(" > GetPrior:: panic, couldn't get prior from Chi2 unfolding! \n");
1292  printf(" probably Chi2 unfolding did not converge < \n");
1293  if(placeHolder) fMergeBinsArray = placeHolder;
1294  return 0x0;
1295  }
1296  if(placeHolder) fMergeBinsArray = placeHolder;
1297  }
1298  break;
1299 
1300  }
1301  case kPriorMeasured : {
1302  unfolded = (TH1D*)measuredJetSpectrumTrueBins->Clone(Form("kPriorMeasured_%s", suffix.Data())); // copy template to unfolded to use as prior
1303  }
1304  default : break;
1305  }
1306  // it can be important that the prior is smooth, this can be achieved by
1307  // extrapolating the spectrum with a fitted power law when bins are sparsely filed
1308  if(jetFindingEfficiency) unfolded->Divide(jetFindingEfficiency);
1309  TH1D *priorLocal((TH1D*)unfolded->Clone(Form("priorUnfolded_%s", suffix.Data())));
1310  if(fSmoothenPrior) priorLocal = SmoothenPrior(priorLocal, fPower, fFitMin, fFitMax, fFitStart, kTRUE, fSmoothenCounts);
1311  return priorLocal;
1312 }
1313 //_____________________________________________________________________________
1315  // resize the x-axis of a th1d
1316 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1317  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1318 #endif
1319  if(!histo) {
1320  printf(" > ResizeXaxisTH!D:: fatal error, NULL pointer passed < \n");
1321  return NULL;
1322  }
1323  // see how many bins we need to copy
1324  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);
1325  // low is the bin number of the first new bin
1326  Int_t l = histo->GetXaxis()->FindBin(low);
1327  // set the values
1328  Double_t _x(0), _xx(0);
1329  for(Int_t i(0); i < up-low; i++) {
1330  _x = histo->GetBinContent(l+i);
1331  _xx=histo->GetBinError(l+i);
1332  resized->SetBinContent(i+1, _x);
1333  resized->SetBinError(i+1, _xx);
1334  }
1335  return resized;
1336 }
1337 //_____________________________________________________________________________
1339  // resize the y-axis of a th2d
1340 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1341  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1342 #endif
1343  if(!histo) {
1344  printf(" > ResizeYaxisTH2D:: fatal error, NULL pointer passed < \n");
1345  return NULL;
1346  }
1347  // see how many bins we need to copy
1348  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());
1349  // assume only the y-axis has changed
1350  // low is the bin number of the first new bin
1351  Int_t low = histo->GetYaxis()->FindBin(y->At(0));
1352  // set the values
1353  Double_t _x(0), _xx(0);
1354  for(Int_t i(0); i < x->GetSize(); i++) {
1355  for(Int_t j(0); j < y->GetSize(); j++) {
1356  _x = histo->GetBinContent(i, low+j);
1357  _xx=histo->GetBinError(i, low+1+j);
1358  resized->SetBinContent(i, j, _x);
1359  resized->SetBinError(i, j, _xx);
1360  }
1361  }
1362  return resized;
1363 }
1364 //_____________________________________________________________________________
1366  // general method to normalize all vertical slices of a th2 to unity
1367  // i.e. get a probability matrix
1368 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1369  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1370 #endif
1371  if(!histo) {
1372  printf(" > NormalizeTH2D:: NULL pointer passed, returning NULL < \n");
1373  return NULL;
1374  }
1375  Int_t binsX = histo->GetXaxis()->GetNbins();
1376  Int_t binsY = histo->GetYaxis()->GetNbins();
1377 
1378  // normalize all slices in x
1379  for(Int_t i(0); i < binsX; i++) { // for each vertical slice
1380  Double_t weight = 0;
1381  for(Int_t j(0); j < binsY; j++) { // loop over all the horizontal components
1382  weight+=histo->GetBinContent(i+1, j+1);
1383  } // now we know the total weight
1384  for(Int_t j(0); j < binsY; j++) {
1385  if (weight <= 0 ) continue;
1386  histo->SetBinContent(1+i, j+1, histo->GetBinContent(1+i, j+1)/weight);
1387  if(noError) histo->SetBinError( 1+i, j+1, 0.);
1388  else histo->SetBinError( 1+i, j+1, histo->GetBinError( 1+i, j+1)/weight);
1389  }
1390  }
1391  return histo;
1392 }
1393 //_____________________________________________________________________________
1395 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1396  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1397 #endif
1398  // resample a TH1
1399  // the returned histogram is new, the original is deleted from the heap if kill is true
1400  if(!hist) {
1401  printf(" > Bootstrap:: fatal error,NULL pointer passed! \n");
1402  return 0x0;
1403  }
1404  else printf(" > Bootstrap:: resampling, may take some time \n");
1405  // clone input histo
1406  TH1* bootstrapped((TH1*)(hist->Clone(Form("%s_bootstrapped", hist->GetName()))));
1407  bootstrapped->Reset();
1408 
1409  /* OLD method - slightly underestimates fluctuations
1410  // reset the content
1411  bootstrapped->Reset();
1412  // resample the input histogram
1413  for(Int_t i(0); i < hist->GetEntries(); i++) bootstrapped->Fill(hist->GetRandom()); */
1414 
1415  // new method
1416  Double_t mean(0), sigma(0);
1417  Int_t sampledMean(0), entries(0);
1418 
1419  for(Int_t i(0); i < hist->GetXaxis()->GetNbins(); i++) {
1420  // for each bin, get the value
1421  mean = hist->GetBinContent(i+1);
1422  sigma = hist->GetBinError(i+1);
1423  // draw a new mean
1424  sampledMean = TMath::Nint(gRandom->Gaus(mean, sigma));
1425  printf(" sampled %i from original number %.2f \n", sampledMean, mean);
1426  // set the new bin content
1427  bootstrapped->SetBinContent(i+1, sampledMean);
1428  if(sampledMean > 0) bootstrapped->SetBinError(i+1, TMath::Sqrt(sampledMean));
1429  entries += sampledMean;
1430  }
1431  printf(" Done bootstrapping, setting number of entries to %i \n", entries);
1432  bootstrapped->SetEntries((double)entries);
1433 
1434 
1435  // if requested kill input histo
1436  if(kill) delete hist;
1437  // return resampled histogram
1438  return bootstrapped;
1439 }
1440 //_____________________________________________________________________________
1442 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1443  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1444 #endif
1445  // return a TH1D with the supplied histogram rebinned to the supplied bins
1446  // the returned histogram is new, the original is deleted from the heap if kill is true
1447  if(!histo || !bins) {
1448  printf(" > RebinTH1D:: fatal error, NULL pointer passed < \n");
1449  return NULL;
1450  }
1451  // create the output histo
1452  TString name = histo->GetName();
1453  name+="_template";
1454  name+=suffix;
1455  TH1D* rebinned = 0x0;
1456  // check if the histogram is normalized. if so, fall back to ROOT style rebinning
1457  // if not, rebin manually
1458  if(histo->GetBinContent(1) > 1 ) {
1459  rebinned = new TH1D(name.Data(), name.Data(), bins->GetSize()-1, bins->GetArray());
1460  for(Int_t i(0); i < histo->GetXaxis()->GetNbins(); i++) {
1461  // loop over the bins of the old histo and fill the new one with its data
1462  rebinned->Fill(histo->GetBinCenter(i+1), histo->GetBinContent(i+1));
1463  }
1464  } else rebinned = reinterpret_cast<TH1D*>(histo->Rebin(bins->GetSize()-1, name.Data(), bins->GetArray()));
1465  if(kill) delete histo;
1466  return rebinned;
1467 }
1468 //_____________________________________________________________________________
1470 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1471  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1472 #endif
1473  // weighted rebinning of a th2d, implementation for function call to AliAnaChargedJetResponseMaker
1474  // not static as it is just a wrapper for the response maker object
1475  if(!fResponseMaker || !binsTrue || !binsRec) {
1476  printf(" > RebinTH2D:: function called with NULL arguments < \n");
1477  return 0x0;
1478  }
1479  TString name(Form("%s_%s", rebinMe->GetName(), suffix.Data()));
1480  return (TH2D*)fResponseMaker->MakeResponseMatrixRebin(rebinMe, (TH2*)(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray())), kTRUE);
1481 }
1482 //_____________________________________________________________________________
1484 {
1485 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1486  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1487 #endif
1488  // multiply two matrices
1489  if (a->GetNbinsX() != b->GetNbinsY()) return 0x0;
1490  TH2D* c = (TH2D*)a->Clone("c");
1491  for (Int_t y1 = 1; y1 <= a->GetNbinsY(); y1++) {
1492  for (Int_t x2 = 1; x2 <= b->GetNbinsX(); x2++) {
1493  Double_t val = 0;
1494  for (Int_t x1 = 1; x1 <= a->GetNbinsX(); x1++) {
1495  Int_t y2 = x1;
1496  val += a->GetBinContent(x1, y1) * b->GetBinContent(x2, y2);
1497  }
1498  c->SetBinContent(x2, y1, val);
1499  c->SetBinError(x2, y1, 0.);
1500  }
1501  }
1502  if(strcmp(name.Data(), "")) c->SetNameTitle(name.Data(), name.Data());
1503  return c;
1504 }
1505 //_____________________________________________________________________________
1507 {
1508 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1509  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1510 #endif
1511  // normalize a th1d to a certain scale
1512  histo->Sumw2();
1513  Double_t integral = histo->Integral()*scale;
1514  if (integral > 0 && scale == 1.) histo->Scale(1./integral, "width");
1515  else if (scale != 1.) histo->Scale(1./scale, "width");
1516  else printf(" > Histogram integral < 0, cannot normalize \n");
1517  return histo;
1518 }
1519 //_____________________________________________________________________________
1521 {
1522 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1523  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1524 #endif
1525  // merge a spectrum histogram taking into account the correlation terms
1526  if(!(bins&&spectrum)) {
1527  printf(" > NULL pointer passed as argument in MergeSpectrumBins ! < \n");
1528  return 0x0;
1529  }
1530  Double_t sum(0), error(0), pearson(0);
1531  // take the sum of the bin content
1532  sum += spectrum->GetBinContent(bins->At(0));
1533  sum += spectrum->GetBinContent(bins->At(1));
1534  // quadratically sum the errors
1535  error += TMath::Power(spectrum->GetBinError(bins->At(0)), 2);
1536  error += TMath::Power(spectrum->GetBinError(bins->At(1)), 2);
1537  // add the covariance term
1538  pearson = corr->GetBinContent(bins->At(0), bins->At(1));
1539  if(!corr) {
1540  printf(" > PANIC ! something is wrong with the covariance matrix, assuming full correlation ! < \n ");
1541  pearson = 1;
1542  }
1543  error += 2.*spectrum->GetBinError(bins->At(0))*spectrum->GetBinError(bins->At(1))*pearson;
1544  spectrum->SetBinContent(1, sum);
1545  if(error <= 0) return spectrum;
1546  spectrum->SetBinError(1, TMath::Sqrt(error));
1547  printf(" > sum is %.2f \t error is %.8f < \n", sum, error);
1548  printf(" > REPLACING BIN CONTENT OF FIRST BIN (%i) WITH MERGED BINS (%i) and (%i) < \n", 1, bins->At(0), bins->At(1));
1549  return spectrum;
1550 }
1551 //_____________________________________________________________________________
1552 TMatrixD* AliJetFlowTools::CalculatePearsonCoefficients(TMatrixD* covarianceMatrix)
1553 {
1554 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1555  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1556 #endif
1557  // Calculate pearson coefficients from covariance matrix
1558  TMatrixD *pearsonCoefficients((TMatrixD*)covarianceMatrix->Clone("pearsonCoefficients"));
1559  Int_t nrows(covarianceMatrix->GetNrows()), ncols(covarianceMatrix->GetNcols());
1560  Double_t pearson(0.);
1561  if(nrows==0 && ncols==0) return 0x0;
1562  for(Int_t row = 0; row < nrows; row++) {
1563  for(Int_t col = 0; col<ncols; col++) {
1564  if((*covarianceMatrix)(row,row)!=0. && (*covarianceMatrix)(col,col)!=0.) pearson = (*covarianceMatrix)(row,col)/TMath::Sqrt((*covarianceMatrix)(row,row)*(*covarianceMatrix)(col,col));
1565  (*pearsonCoefficients)(row,col) = pearson;
1566  }
1567  }
1568  return pearsonCoefficients;
1569 }
1570 //_____________________________________________________________________________
1571 TH1D* AliJetFlowTools::SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill, Bool_t counts) {
1572 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1573  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1574 #endif
1575  // smoothen the spectrum using a user defined function
1576  // returns a clone of the original spectrum if fitting failed
1577  // if kill is kTRUE the input spectrum will be deleted from the heap
1578  // if 'count' is selected, bins are filled with integers (necessary if the
1579  // histogram is interpreted in a routine which accepts only counts)
1580  if(!spectrum || !function) return 0x0;
1581  if(start > max) printf(" > cannot extrapolate fit beyond fit range ! < " );
1582  TH1D* temp = (TH1D*)spectrum->Clone(Form("%s_smoothened", spectrum->GetName()));
1583  temp->Sumw2(); // if already called on the original, this will give off a warning but do nothing
1584  TFitResultPtr r = temp->Fit(function, "", "", min, max);
1585  if((int)r == 0) { // MINUIT status
1586  for(Int_t i(0); i < temp->GetNbinsX() + 1; i++) {
1587  if(temp->GetBinCenter(i) > start) { // from this pt value use extrapolation
1588  (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));
1589  if(temp->GetBinContent(i) > 0) temp->SetBinError(i, TMath::Sqrt(temp->GetBinContent(i)));
1590  }
1591  }
1592  }
1593  if(kill) delete spectrum;
1594  return temp;
1595 }
1596 //_____________________________________________________________________________
1598 {
1599 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1600  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1601 #endif
1602  // set global style for your current aliroot session
1603  if(!gStyle) return;
1604  // legacy style is pleasing to the eye, default is the formal ALICE style
1605  if(legacy) {
1606  gStyle->SetCanvasColor(-1);
1607  gStyle->SetPadColor(-1);
1608  gStyle->SetFrameFillColor(-1);
1609  gStyle->SetHistFillColor(-1);
1610  gStyle->SetTitleFillColor(-1);
1611  gStyle->SetFillColor(-1);
1612  gStyle->SetFillStyle(4000);
1613  gStyle->SetStatStyle(0);
1614  gStyle->SetTitleStyle(0);
1615  gStyle->SetCanvasBorderSize(0);
1616  gStyle->SetFrameBorderSize(0);
1617  gStyle->SetLegendBorderSize(0);
1618  gStyle->SetStatBorderSize(0);
1619  gStyle->SetTitleBorderSize(0);
1620  } else {
1621  gStyle->Reset("Plain");
1622  gStyle->SetOptTitle(0);
1623  gStyle->SetOptStat(0);
1624  gStyle->SetPalette(1);
1625  gStyle->SetCanvasColor(10);
1626  gStyle->SetCanvasBorderMode(0);
1627  gStyle->SetFrameLineWidth(1);
1628  gStyle->SetFrameFillColor(kWhite);
1629  gStyle->SetPadColor(10);
1630  gStyle->SetPadTickX(1);
1631  gStyle->SetPadTickY(1);
1632  gStyle->SetPadBottomMargin(0.17);
1633  gStyle->SetPadLeftMargin(0.15);
1634  gStyle->SetHistLineWidth(1);
1635  gStyle->SetHistLineColor(kRed);
1636  gStyle->SetFuncWidth(2);
1637  gStyle->SetFuncColor(kGreen);
1638  gStyle->SetLineWidth(2);
1639  gStyle->SetLabelSize(0.045,"xyz");
1640  gStyle->SetLabelOffset(0.01,"y");
1641  gStyle->SetLabelOffset(0.01,"x");
1642  gStyle->SetLabelColor(kBlack,"xyz");
1643  gStyle->SetTitleSize(0.05,"xyz");
1644  gStyle->SetTitleOffset(1.25,"y");
1645  gStyle->SetTitleOffset(1.2,"x");
1646  gStyle->SetTitleFillColor(kWhite);
1647  gStyle->SetTextSizePixels(26);
1648  gStyle->SetTextFont(42);
1649  gStyle->SetLegendBorderSize(0);
1650  gStyle->SetLegendFillColor(kWhite);
1651  gStyle->SetLegendFont(42);
1652  }
1653 }
1654 //_____________________________________________________________________________
1655 void AliJetFlowTools::Style(TCanvas* c, TString style)
1656 {
1657 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1658  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1659 #endif
1660  if(c) return;
1661  // set a default style for a canvas
1662  if(!strcmp(style.Data(), "PEARSON")) {
1663  printf(" > style PEARSON canvas < \n");
1664  gStyle->SetOptStat(0);
1665  c->SetGridx();
1666  c->SetGridy();
1667  c->SetTicks();
1668  return;
1669  } else if(!strcmp(style.Data(), "SPECTRUM")) {
1670  printf(" > style SPECTRUM canvas < \n");
1671  gStyle->SetOptStat(0);
1672  c->SetLogy();
1673  c->SetGridx();
1674  c->SetGridy();
1675  c->SetTicks();
1676  return;
1677  } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1678 }
1679 //_____________________________________________________________________________
1680 void AliJetFlowTools::Style(TVirtualPad* c, TString style, Bool_t legacy)
1681 {
1682 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1683  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1684 #endif
1685  // set a default style for a canva
1686  if(!c) return;
1687  if(legacy) {
1688  c->SetLeftMargin(.25);
1689  c->SetBottomMargin(.25);
1690  }
1691  else Style();
1692  if(!strcmp(style.Data(), "PEARSON")) {
1693  printf(" > style PEARSON pad < \n");
1694  gStyle->SetOptStat(0);
1695  c->SetGridx();
1696  c->SetGridy();
1697  c->SetTicks();
1698  return;
1699  } else if(!strcmp(style.Data(), "SPECTRUM")) {
1700  printf(" > style SPECTRUM pad < \n");
1701  gStyle->SetOptStat(0);
1702  c->SetLogy();
1703  c->SetGridx();
1704  c->SetGridy();
1705  c->SetTicks();
1706  return;
1707  } else if (!strcmp(style.Data(), "GRID")) {
1708  printf(" > style GRID pad < \n");
1709  gStyle->SetOptStat(0);
1710  c->SetGridx();
1711  c->SetGridy();
1712  c->SetTicks();
1713  } else printf(" > Style called with unknown option %s \n returning < \n", style.Data());
1714 }
1715 //_____________________________________________________________________________
1716 void AliJetFlowTools::Style(TLegend* l)
1717 {
1718 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1719  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1720 #endif
1721  // set a default style for a legend
1722  if(!l) return;
1723  l->SetFillColor(0);
1724  l->SetBorderSize(0);
1725  if(gStyle) l->SetTextSize(gStyle->GetTextSize()*.08);
1726 }
1727 //_____________________________________________________________________________
1728 void AliJetFlowTools::Style(TH1* h, EColor col, histoType type, Bool_t legacy)
1729 {
1730 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1731  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1732 #endif
1733  // style a histo
1734  if(!h) return;
1735  h->GetYaxis()->SetNdivisions(505);
1736  h->SetLineColor(col);
1737  h->SetMarkerColor(col);
1738  h->SetLineWidth(2);
1739  h->SetMarkerSize(1.5);
1740  if(legacy) {
1741  h->SetTitle("");
1742  h->GetYaxis()->SetLabelSize(0.05);
1743  h->GetXaxis()->SetLabelSize(0.05);
1744  h->GetYaxis()->SetTitleOffset(1.5);
1745  h->GetXaxis()->SetTitleOffset(1.5);
1746  h->GetYaxis()->SetTitleSize(.05);
1747  h->GetXaxis()->SetTitleSize(.05);
1748  } else Style();
1749  switch (type) {
1750  case kInPlaneSpectrum : {
1751  h->SetTitle("IN PLANE");
1752  h->GetXaxis()->SetTitle("#it{p}_{T}^{ch, jet} (GeV/#it{c})");
1753  h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1754  } break;
1755  case kOutPlaneSpectrum : {
1756  h->SetTitle("OUT OF PLANE");
1757  h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1758  h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1759  } break;
1760  case kUnfoldedSpectrum : {
1761  h->SetTitle("UNFOLDED");
1762  h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1763  h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1764  } break;
1765  case kFoldedSpectrum : {
1766  h->SetTitle("FOLDED");
1767  h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1768  h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1769  } break;
1770  case kMeasuredSpectrum : {
1771  h->SetTitle("MEASURED");
1772  h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1773  h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{cht} (GeV/#it{c})");
1774  } break;
1775  case kBar : {
1776  h->SetFillColor(col);
1777  h->SetBarWidth(.6);
1778  h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1779  h->SetBarOffset(0.2);
1780  }
1781  case kRatio : {
1782  h->SetMarkerStyle(8);
1783  h->SetMarkerSize(1.5);
1784  } break;
1785  case kDeltaPhi : {
1786  h->GetYaxis()->SetTitle("[counts]");
1787  h->GetXaxis()->SetTitle("#Delta #phi");
1788  }
1789  default : break;
1790  }
1791 }
1792 //_____________________________________________________________________________
1793 void AliJetFlowTools::Style(TGraph* h, EColor col, histoType type, Bool_t legacy)
1794 {
1795 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1796  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1797 #endif
1798  // style a tgraph
1799  if(!h) return;
1800  h->GetYaxis()->SetNdivisions(505);
1801  h->SetLineColor(col);
1802  h->SetMarkerColor(col);
1803  h->SetLineWidth(2);
1804  h->SetMarkerSize(1.5);
1805  h->SetTitle("");
1806  h->SetFillColor(kCyan);
1807  if(legacy) {
1808  h->GetYaxis()->SetLabelSize(0.05);
1809  h->GetXaxis()->SetLabelSize(0.05);
1810  h->GetYaxis()->SetTitleOffset(1.6);
1811  h->GetXaxis()->SetTitleOffset(1.6);
1812  h->GetYaxis()->SetTitleSize(.05);
1813  h->GetXaxis()->SetTitleSize(.05);
1814  } else Style();
1815  h->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
1816  switch (type) {
1817  case kInPlaneSpectrum : {
1818  h->SetTitle("IN PLANE");
1819  h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1820  } break;
1821  case kOutPlaneSpectrum : {
1822  h->SetTitle("OUT OF PLANE");
1823  h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1824  } break;
1825  case kUnfoldedSpectrum : {
1826  h->SetTitle("UNFOLDED");
1827  h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1828  } break;
1829  case kFoldedSpectrum : {
1830  h->SetTitle("FOLDED");
1831  h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1832  } break;
1833  case kMeasuredSpectrum : {
1834  h->SetTitle("MEASURED");
1835  h->GetYaxis()->SetTitle("#frac{d#it{N}}{d#it{p}_{T}}");
1836  } break;
1837  case kRatio : {
1838  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))");
1839  } break;
1840  case kV2 : {
1841  h->GetYaxis()->SetTitle("#it{v}_{2}^{ch, jet} {EP, |#Delta#eta|>0.9 } ");
1842  h->GetYaxis()->SetRangeUser(-.5, 1.);
1843  h->SetMarkerStyle(8);
1844  h->SetMarkerSize(1.5);
1845  } break;
1846  default : break;
1847  }
1848 }
1849 //_____________________________________________________________________________
1851  TH1D*& ratio, // pointer reference, output of this function
1852  TGraphErrors*& v2, // pointer reference, as output of this function
1853  TArrayI* in,
1854  TArrayI* out,
1855  TString inFile,
1856  TString outFile) const
1857 {
1858 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1859  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1860 #endif
1861  // pass clones of the nominal points and nominal v2 values
1862  if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close(); // if for some weird reason the unfolding output is still mutable
1863  TFile* readMe(new TFile(inFile.Data(), "READ")); // open unfolding output read-only
1864  if(readMe->IsZombie()) {
1865  printf(" > Fatal error, couldn't read %s for post processing ! < \n", inFile.Data());
1866  return;
1867  }
1868  printf("\n\n\n\t\t GetNominalValues \n > Recovered the following file structure : \n <");
1869  readMe->ls();
1870  TFile* output(new TFile(outFile.Data(), "RECREATE")); // create new output file
1871  if(output) SquelchWarning();
1872  // get some placeholders, have to be initialized but will be deleted
1873  ratio = new TH1D("nominal", "nominal", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
1874  TH1D* nominalIn(new TH1D("nominal in", "nominal in", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1875  TH1D* nominalOut(new TH1D("nominal out", "nominal out", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1876  Int_t iIn[] = {in->At(0), in->At(0)}; // first value is the nominal point
1877  Int_t iOut[] = {out->At(0), out->At(0)};
1878 
1879  // call the functions and set the relevant pointer references
1880  TH1D* dud(0x0);
1882  new TArrayI(2, iIn),
1883  new TArrayI(2, iOut),
1884  dud, dud, dud, dud, dud, dud,
1885  ratio, // pointer reference, output of this function
1886  nominalIn,
1887  nominalOut,
1888  1,
1889  fBinsTrue->At(0),
1890  fBinsTrue->At(fBinsTrue->GetSize()-1),
1891  readMe,
1892  "nominal_values");
1893  v2 = GetV2(nominalIn, nominalOut, fEventPlaneRes, "nominal v_{2}");
1894 
1895  // close the open files, reclaim ownership of histograms which are necessary outside of the file scope
1896  ratio->SetDirectory(0); // disassociate from current gDirectory
1897  readMe->Close();
1898 }
1899 //_____________________________________________________________________________
1901  TGraphAsymmErrors*& corrRatio, // correlated uncertainty pointer reference
1902  TGraphAsymmErrors*& corrV2, // correlated uncertainty pointer reference
1903  TArrayI* variationsIn, // variantions in plane
1904  TArrayI* variationsOut, // variantions out of plane
1905  Bool_t sym, // treat as symmmetric
1906  TArrayI* variations2ndIn, // second source of variations
1907  TArrayI* variations2ndOut, // second source of variations
1908  Bool_t sym2nd, // treat as symmetric
1909  TString type, // type of uncertaitny
1910  TString type2,
1911  Int_t columns, // divide the output canvasses in this many columns
1912  Float_t rangeLow, // lower pt range
1913  Float_t rangeUp, // upper pt range
1914  Float_t corr, // correlation strength
1915  TString in, // input file name (created by this unfolding class)
1916  TString out // output file name (which will hold results of the systematic test)
1917  ) const
1918 {
1919 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
1920  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
1921 #endif
1922  // do full systematics
1923  if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close(); // if for some weird reason the unfolding output is still mutable
1924  TFile* readMe(new TFile(in.Data(), "READ")); // open unfolding output read-only
1925  if(readMe->IsZombie()) {
1926  printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
1927  return;
1928  }
1929  printf("\n\n\n\t\t GetCorrelatedUncertainty \n > Recovered the following file structure : \n <");
1930  readMe->ls();
1931  TFile* output(new TFile(out.Data(), "RECREATE")); // create new output file
1932 
1933  if(sym2nd) SquelchWarning();
1934 
1935  // create some null placeholder pointers
1936  TH1D* relativeErrorVariationInLow(0x0);
1937  TH1D* relativeErrorVariationInUp(0x0);
1938  TH1D* relativeErrorVariationOutLow(0x0);
1939  TH1D* relativeErrorVariationOutUp(0x0);
1940  TH1D* relativeError2ndVariationInLow(0x0);
1941  TH1D* relativeError2ndVariationInUp(0x0);
1942  TH1D* relativeError2ndVariationOutLow(0x0);
1943  TH1D* relativeError2ndVariationOutUp(0x0);
1944  TH1D* relativeStatisticalErrorIn(0x0);
1945  TH1D* relativeStatisticalErrorOut(0x0);
1946  // histo for the nominal ratio in / out
1947  TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1948  TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1949  TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
1950 
1951  // call the functions
1952  if(variationsIn && variationsOut) {
1954  variationsIn,
1955  variationsOut,
1956  relativeErrorVariationInUp, // pointer reference
1957  relativeErrorVariationInLow, // pointer reference
1958  relativeErrorVariationOutUp, // pointer reference
1959  relativeErrorVariationOutLow, // pointer reference
1960  relativeStatisticalErrorIn, // pointer reference
1961  relativeStatisticalErrorOut, // pointer reference
1962  nominal,
1963  nominalIn,
1964  nominalOut,
1965  columns,
1966  rangeLow,
1967  rangeUp,
1968  readMe,
1969  type,
1970  kFALSE,
1971  kTRUE);
1972  if(relativeErrorVariationInUp) {
1973  // canvas with the error from variation strength
1974  TCanvas* relativeErrorVariation(new TCanvas(Form("relativeError_%s", type.Data()), Form("relativeError_%s", type.Data())));
1975  relativeErrorVariation->Divide(2);
1976  relativeErrorVariation->cd(1);
1977  Style(gPad, "GRID");
1978  relativeErrorVariationInUp->DrawCopy("b");
1979  relativeErrorVariationInLow->DrawCopy("same b");
1980  Style(AddLegend(gPad));
1981  relativeErrorVariation->cd(2);
1982  Style(gPad, "GRID");
1983  relativeErrorVariationOutUp->DrawCopy("b");
1984  relativeErrorVariationOutLow->DrawCopy("same b");
1985  SavePadToPDF(relativeErrorVariation);
1986  Style(AddLegend(gPad));
1987  relativeErrorVariation->Write();
1988 
1989  // now smoothen the diced response error (as it is expected to be flat)
1990  // this is done by fitting a constant to the diced resonse error histo
1991  //
1992  TF1* lin = new TF1("lin", "[0]", rangeLow, rangeUp);
1993  relativeErrorVariationInUp->Fit(lin, "L", "", rangeLow, rangeUp);
1994  if(gMinuit->fISW[1] != 3) printf(" fit is NOT ok ! " );
1995  for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
1996  relativeErrorVariationInUp->SetBinContent(i+1, lin->GetParameter(0));
1997  }
1998  relativeErrorVariationInLow->Fit(lin, "L", "", rangeLow, rangeUp);
1999  printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
2000  for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
2001  relativeErrorVariationInLow->SetBinContent(i+1, 0);//lin->GetParameter(0));
2002  }
2003  relativeErrorVariationOutUp->Fit(lin, "L", "", rangeLow, rangeUp);
2004  printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
2005  for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
2006  relativeErrorVariationOutUp->SetBinContent(i+1, lin->GetParameter(0));
2007  }
2008  relativeErrorVariationOutLow->Fit(lin, "L", "", rangeLow, rangeUp);
2009  printf(" > Fit over diced resonse, new value for all bins is %.4f < \n ", lin->GetParameter(0));
2010  for(Int_t i(0); i < relativeErrorVariationInUp->GetNbinsX(); i++) {
2011  relativeErrorVariationOutLow->SetBinContent(i+1, 0);//lin->GetParameter(0));
2012  }
2013  }
2014  }
2015  // call the functions for a second set of systematic sources
2016  if(variations2ndIn && variations2ndOut) {
2018  variations2ndIn,
2019  variations2ndOut,
2020  relativeError2ndVariationInUp, // pointer reference
2021  relativeError2ndVariationInLow, // pointer reference
2022  relativeError2ndVariationOutUp, // pointer reference
2023  relativeError2ndVariationOutLow, // pointer reference
2024  relativeStatisticalErrorIn, // pointer reference
2025  relativeStatisticalErrorOut, // pointer reference
2026  nominal,
2027  nominalIn,
2028  nominalOut,
2029  columns,
2030  rangeLow,
2031  rangeUp,
2032  readMe,
2033  type2,
2034  kFALSE,
2035  kTRUE);
2036  if(relativeError2ndVariationInUp) {
2037  // canvas with the error from variation strength
2038  TCanvas* relativeError2ndVariation(new TCanvas(Form("relativeError2nd_%s", type2.Data()), Form("relativeError2nd_%s", type2.Data())));
2039  relativeError2ndVariation->Divide(2);
2040  relativeError2ndVariation->cd(1);
2041  Style(gPad, "GRID");
2042  relativeError2ndVariationInUp->DrawCopy("b");
2043  relativeError2ndVariationInLow->DrawCopy("same b");
2044  Style(AddLegend(gPad));
2045  relativeError2ndVariation->cd(2);
2046  Style(gPad, "GRID");
2047  relativeError2ndVariationOutUp->DrawCopy("b");
2048  relativeError2ndVariationOutLow->DrawCopy("same b");
2049  SavePadToPDF(relativeError2ndVariation);
2050  Style(AddLegend(gPad));
2051  relativeError2ndVariation->Write();
2052  }
2053 
2054  }
2055 
2056  // and the placeholder for the final systematic
2057  TH1D* relativeErrorInUp(new TH1D("max correlated uncertainty in plane", "max correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2058  TH1D* relativeErrorOutUp(new TH1D("max correlated uncertainty out of plane", "max correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2059  TH1D* relativeErrorInLow(new TH1D("min correlated uncertainty in plane", "min correlated uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2060  TH1D* relativeErrorOutLow(new TH1D("min correlated uncertainty out of plane", "min correlated uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2061  relativeErrorInUp->SetYTitle("relative uncertainty");
2062  relativeErrorOutUp->SetYTitle("relative uncertainty");
2063  relativeErrorInLow->SetYTitle("relative uncertainty");
2064  relativeErrorOutLow->SetYTitle("relative uncertainty");
2065 
2066  // merge the correlated errors (FIXME) trivial for one set
2067  Double_t aInUp(0.), bInUp(0.), cInUp(0.), dInUp(0.);
2068  Double_t aOutUp(0.), bOutUp(0.), cOutUp(0.), dOutUp(0.);
2069  Double_t aInLow(0.), bInLow(0.), cInLow(0.), dInLow(0.);
2070  Double_t aOutLow(0.), bOutLow(0.), cOutLow(0.), dOutLow(0.);
2071 
2072  for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2073  // for the upper bound
2074  if(relativeErrorVariationInUp) aInUp = relativeErrorVariationInUp->GetBinContent(b+1);
2075  if(relativeErrorVariationOutUp) aOutUp = relativeErrorVariationOutUp->GetBinContent(b+1);
2076  if(relativeError2ndVariationInUp) bInUp = relativeError2ndVariationInUp->GetBinContent(b+1);
2077  if(relativeError2ndVariationOutUp) bInLow = relativeError2ndVariationOutUp->GetBinContent(b+1);
2078  dInUp = aInUp*aInUp + bInUp*bInUp + cInUp*cInUp;
2079  // for a symmetric first variation
2080  if(sym) dInUp += aInLow*aInLow;
2081  if(dInUp > 0) relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(dInUp));
2082  dOutUp = aOutUp*aOutUp + bOutUp*bOutUp + cOutUp*cOutUp;
2083  if(sym) dOutUp += aOutLow*aOutLow;
2084  if(dOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(dOutUp));
2085  // for the lower bound
2086  if(relativeErrorVariationInLow) aInLow = relativeErrorVariationInLow->GetBinContent(b+1);
2087  if(relativeErrorVariationOutLow) aOutLow = relativeErrorVariationOutLow->GetBinContent(b+1);
2088  if(relativeError2ndVariationInLow) bInLow = relativeError2ndVariationInLow->GetBinContent(b+1);
2089  if(relativeError2ndVariationOutLow) bOutLow = relativeError2ndVariationOutLow->GetBinContent(b+1);
2090  dInLow = aInLow*aInLow + bInLow*bInLow + cInLow*cInLow;
2091  if(sym) dInLow += aInUp*aInUp;
2092  if(dInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1*TMath::Sqrt(dInLow));
2093  dOutLow = aOutLow*aOutLow + bOutLow*bOutLow + cOutLow*cOutLow;
2094  if(sym) dOutLow += aOutUp*aOutUp;
2095  if(dOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(dOutLow));
2096  }
2097  // project the estimated errors on the nominal ratio
2098  if(nominal) {
2099  Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
2100  Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
2101  Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
2102  Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
2103  Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
2104  Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
2105  Double_t lowErr(0.), upErr(0.);
2106  for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
2107  // add the in and out of plane errors in quadrature
2108  // propagation is tricky because we should consider two cases:
2109  // [1] the error is uncorrelated, and we propagate upper 1 with lower 2 and vice versa
2110  // [2] the error is correlated - in this case upper 1 should be propagated with upper 2
2111  // as the fluctuations are bound to always to of same sign
2112  if(corr <= 0) {
2113  lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(1+i)*relativeErrorOutUp->GetBinContent(i+1);
2114  upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
2115  } else {
2116  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);
2117  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);
2118  }
2119  ayl[i] = TMath::Sqrt(TMath::Abs(lowErr))*nominal->GetBinContent(i+1);
2120  ayh[i] = TMath::Sqrt(TMath::Abs(upErr))*nominal->GetBinContent(i+1);
2121  // get the bin width (which is the 'error' on x
2122  Double_t binWidth(nominal->GetBinWidth(i+1));
2123  axl[i] = binWidth/2.;
2124  axh[i] = binWidth/2.;
2125  // now get the coordinate for the point
2126  ax[i] = nominal->GetBinCenter(i+1);
2127  ay[i] = nominal->GetBinContent(i+1);
2128  }
2129  // save the nominal ratio
2130  TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
2131  corrRatio = (TGraphAsymmErrors*)nominalError->Clone();
2132  nominalError->SetName("correlated uncertainty");
2133  TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
2134  nominalCanvas->Divide(2);
2135  nominalCanvas->cd(1);
2136  Style(nominal, kBlack);
2137  Style(nominalError, kCyan, kRatio);
2138  nominalError->SetLineColor(kCyan);
2139  nominalError->SetMarkerColor(kCyan);
2140  nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2141  nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
2142  nominalError->DrawClone("a2");
2143  nominal->DrawCopy("same E1");
2144  Style(AddLegend(gPad));
2145  Style(gPad, "GRID");
2146  Style(nominalCanvas);
2147  // save nominal v2 and systematic errors
2149  nominalIn,
2150  nominalOut,
2152  "v_{2} with shape uncertainty",
2153  relativeErrorInUp,
2154  relativeErrorInLow,
2155  relativeErrorOutUp,
2156  relativeErrorOutLow,
2157  corr));
2158  // pass the nominal values to the pointer references
2159  corrV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
2160  TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
2161  nominalCanvas->cd(2);
2162  Style(nominalV2, kBlack);
2163  Style(nominalV2Error, kCyan, kV2);
2164  nominalV2Error->SetLineColor(kCyan);
2165  nominalV2Error->SetMarkerColor(kCyan);
2166  nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2167  nominalV2Error->DrawClone("a2");
2168  nominalV2->DrawClone("same E1");
2169  Style(AddLegend(gPad));
2170  Style(gPad, "GRID");
2171  Style(nominalCanvas);
2172  SavePadToPDF(nominalCanvas);
2173  nominalCanvas->Write();
2174  }
2175 
2176  TCanvas* relativeError(new TCanvas("relativeCorrelatedError"," relativeCorrelatedError"));
2177  relativeError->Divide(2);
2178  relativeError->cd(1);
2179  Style(gPad, "GRID");
2180  relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2181  Style(relativeErrorInUp, kBlue, kBar);
2182  Style(relativeErrorInLow, kGreen, kBar);
2183  relativeErrorInUp->DrawCopy("b");
2184  relativeErrorInLow->DrawCopy("same b");
2185  Style(relativeStatisticalErrorIn, kRed);
2186  relativeStatisticalErrorIn->DrawCopy("same");
2187  Style(AddLegend(gPad));
2188  relativeError->cd(2);
2189  Style(gPad, "GRID");
2190  relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2191  Style(relativeErrorOutUp, kBlue, kBar);
2192  Style(relativeErrorOutLow, kGreen, kBar);
2193  relativeErrorOutUp->DrawCopy("b");
2194  relativeErrorOutLow->DrawCopy("same b");
2195  Style(relativeStatisticalErrorOut, kRed);
2196  relativeStatisticalErrorOut->DrawCopy("same");
2197  Style(AddLegend(gPad));
2198 
2199  // write the buffered file to disk and close the file
2200  SavePadToPDF(relativeError);
2201  relativeError->Write();
2202  output->Write();
2203 }
2204 //_____________________________________________________________________________
2206  TGraphAsymmErrors*& shapeRatio, // pointer reference to final shape uncertainty of ratio
2207  TGraphAsymmErrors*& shapeV2, // pointer reference to final shape uncertainty of v2
2208  TArrayI* regularizationIn, // regularization strength systematics, in plane
2209  TArrayI* regularizationOut, // regularization strenght systematics, out of plane
2210  TArrayI* recBinIn, // rec bin ranges, in plane
2211  TArrayI* recBinOut, // rec bin ranges, out of plane
2212  TArrayI* methodIn, // method in
2213  TArrayI* methodOut, // method out
2214  Int_t columns, // divide the output canvasses in this many columns
2215  Float_t rangeLow, // lower pt range
2216  Float_t rangeUp, // upper pt range
2217  Float_t corr, // correlation strength
2218  TString in, // input file name (created by this unfolding class)
2219  TString out, // output file name (which will hold results of the systematic test)
2220  Bool_t regularizationOnV2 // get uncertainty on yields separately or v2 directly
2221  ) const
2222 {
2223 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
2224  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2225 #endif
2226  // do full systematics
2227  if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close(); // if for some weird reason the unfolding output is still mutable
2228  TFile* readMe(new TFile(in.Data(), "READ")); // open unfolding output read-only
2229  if(readMe->IsZombie()) {
2230  printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
2231  return;
2232  }
2233  printf("\n\n\n\t\t DOSYSTEMATICS \n > Recovered the following file structure : \n <");
2234  readMe->ls();
2235  TFile* output(new TFile(out.Data(), "RECREATE")); // create new output file
2236  // create some null placeholder pointers
2237  TH1D* relativeErrorRegularizationInLow(0x0);
2238  TH1D* relativeErrorRegularizationInUp(0x0);
2239  TH1D* relativeErrorRecBinInLow(0x0);
2240  TH1D* relativeErrorRecBinInUp(0x0);
2241  TH1D* relativeErrorMethodInLow(0x0);
2242  TH1D* relativeErrorMethodInUp(0x0);
2243  TH1D* relativeErrorRegularizationOutLow(0x0);
2244  TH1D* relativeErrorRegularizationOutUp(0x0);
2245  TH1D* relativeErrorRecBinOutLow(0x0);
2246  TH1D* relativeErrorRecBinOutUp(0x0);
2247  TH1D* relativeStatisticalErrorIn(0x0);
2248  TH1D* relativeStatisticalErrorOut(0x0);
2249  TH1D* relativeErrorMethodOutLow(0x0);
2250  TH1D* relativeErrorMethodOutUp(0x0);
2251  // histo for the nominal ratio in / out
2252  TH1D* nominal(new TH1D("ratio in plane, out of plane", "ratio in plane, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2253  TH1D* nominalIn(new TH1D("in plane jet yield", "in plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2254  TH1D* nominalOut(new TH1D("out of plane jet yield", "out of plane jet yield", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2255 
2256  // call the functions
2257  if(regularizationIn && regularizationOut) {
2259  regularizationIn,
2260  regularizationOut,
2261  relativeErrorRegularizationInUp, // pointer reference
2262  relativeErrorRegularizationInLow, // pointer reference
2263  relativeErrorRegularizationOutUp, // pointer reference
2264  relativeErrorRegularizationOutLow, // pointer reference
2265  relativeStatisticalErrorIn, // pointer reference
2266  relativeStatisticalErrorOut, // pointer reference
2267  nominal,
2268  nominalIn,
2269  nominalOut,
2270  columns,
2271  rangeLow,
2272  rangeUp,
2273  readMe,
2274  "regularization",
2275  fRMS,
2276  !regularizationOnV2); // set to false means NO undertainty from ratio, but v2
2277  if(relativeErrorRegularizationInUp && !regularizationOnV2 ) {
2278  // canvas with the error from regularization strength
2279  TCanvas* relativeErrorRegularization(new TCanvas("relativeErrorRegularization", "relativeErrorRegularization"));
2280  relativeErrorRegularization->Divide(2);
2281  relativeErrorRegularization->cd(1);
2282  Style(gPad, "GRID");
2283  relativeErrorRegularizationInUp->DrawCopy("b");
2284  relativeErrorRegularizationInLow->DrawCopy("same b");
2285  Style(AddLegend(gPad));
2286  relativeErrorRegularization->cd(2);
2287  Style(gPad, "GRID");
2288  relativeErrorRegularizationOutUp->DrawCopy("b");
2289  relativeErrorRegularizationOutLow->DrawCopy("same b");
2290  SavePadToPDF(relativeErrorRegularization);
2291  Style(AddLegend(gPad));
2292  relativeErrorRegularization->Write();
2293  } else if(relativeErrorRegularizationInUp && regularizationOnV2 ) {
2294  // canvas with the error from regularization strength
2295  TCanvas* relativeErrorRegularization(new TCanvas("relativeErrorRegularization", "relativeErrorRegularization"));
2296  Style(gPad, "GRID");
2297  relativeErrorRegularization->cd(1);
2298  TH1F* relativeErrorRegularizationInUpErrors = (TH1F*)relativeErrorRegularizationInUp->Clone();
2299  for(Int_t i(1); i < relativeErrorRegularizationInUp->GetNbinsX() + 1; i++) {
2300  relativeErrorRegularizationInUpErrors->SetBinContent(i, relativeErrorRegularizationInUp->GetBinError(i));
2301  relativeErrorRegularizationInUpErrors->SetBinError(i, 0);
2302  }
2303  relativeErrorRegularizationInUpErrors->GetYaxis()->SetTitle("absolute error on v_{2} from unfolding");
2304  relativeErrorRegularizationInUpErrors->GetXaxis()->SetTitle("#it{p}_{T, jet}^{ch} (GeV/#it{c})");
2305  relativeErrorRegularizationInUpErrors->DrawCopy("b");
2306  Style(gPad, "GRID");
2307  relativeErrorRegularization->Write();
2308  }
2309  }
2310  if(recBinIn && recBinOut) {
2312  recBinIn,
2313  recBinOut,
2314  relativeErrorRecBinInUp, // pointer reference
2315  relativeErrorRecBinInLow, // pointer reference
2316  relativeErrorRecBinOutUp, // pointer reference
2317  relativeErrorRecBinOutLow, // pointer reference
2318  relativeStatisticalErrorIn,
2319  relativeStatisticalErrorOut,
2320  nominal,
2321  nominalIn,
2322  nominalOut,
2323  columns,
2324  rangeLow,
2325  rangeUp,
2326  readMe,
2327  "recBin",
2328  kFALSE,//fRMS, // was RMS< i think not ok ...
2329  kTRUE); // undertainty from ratio
2330  if(relativeErrorRecBinOutUp) {
2331  // canvas with the error from regularization strength
2332  TCanvas* relativeErrorRecBin(new TCanvas("relativeErrorRecBin"," relativeErrorRecBin"));
2333  relativeErrorRecBin->Divide(2);
2334  relativeErrorRecBin->cd(1);
2335  Style(gPad, "GRID");
2336  relativeErrorRecBinInUp->DrawCopy("b");
2337  relativeErrorRecBinInLow->DrawCopy("same b");
2338  Style(AddLegend(gPad));
2339  relativeErrorRecBin->cd(2);
2340  Style(gPad, "GRID");
2341  relativeErrorRecBinOutUp->DrawCopy("b");
2342  relativeErrorRecBinOutLow->DrawCopy("same b");
2343  SavePadToPDF(relativeErrorRecBin);
2344  Style(AddLegend(gPad));
2345  relativeErrorRecBin->Write();
2346  }
2347  }
2348  if(methodIn && methodOut) {
2350  methodIn,
2351  methodOut,
2352  relativeErrorMethodInUp, // pointer reference
2353  relativeErrorMethodInLow, // pointer reference
2354  relativeErrorMethodOutUp, // pointer reference
2355  relativeErrorMethodOutLow, // pointer reference
2356  relativeStatisticalErrorIn,
2357  relativeStatisticalErrorOut,
2358  nominal,
2359  nominalIn,
2360  nominalOut,
2361  columns,
2362  rangeLow,
2363  rangeUp,
2364  readMe,
2365  "method",
2366  kFALSE,
2367  kTRUE); // undertaitny from ratio
2368  if(relativeErrorMethodInUp) {
2369  TCanvas* relativeErrorMethod(new TCanvas("relativeErrorMethod", "relativeErrorMethod"));
2370  relativeErrorMethod->Divide(2);
2371  relativeErrorMethod->cd(1);
2372  Style(gPad, "GRID");
2373  relativeErrorMethodInUp->DrawCopy("b");
2374  relativeErrorMethodInLow->DrawCopy("same b");
2375  Style(AddLegend(gPad));
2376  relativeErrorMethod->cd(2);
2377  Style(gPad, "GRID");
2378  relativeErrorMethodOutUp->DrawCopy("b");
2379  relativeErrorMethodOutLow->DrawCopy("same b");
2380  SavePadToPDF(relativeErrorMethod);
2381  Style(AddLegend(gPad));
2382  relativeErrorMethod->Write();
2383  }
2384  }
2385 
2386  // and the placeholder for the final systematic
2387  TH1D* relativeErrorInUp(new TH1D("max shape uncertainty in plane", "max shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2388  TH1D* relativeErrorOutUp(new TH1D("max shape uncertainty out of plane", "max shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2389  TH1D* relativeErrorInLow(new TH1D("min shape uncertainty in plane", "min shape uncertainty in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2390  TH1D* relativeErrorOutLow(new TH1D("min shape uncertainty out of plane", "min shape uncertainty out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2391  relativeErrorInUp->SetYTitle("relative uncertainty");
2392  relativeErrorOutUp->SetYTitle("relative uncertainty");
2393  relativeErrorInLow->SetYTitle("relative uncertainty");
2394  relativeErrorOutLow->SetYTitle("relative uncertainty");
2395 
2396  // sum of squares for the total systematic error
2397  Double_t aInUp(0.), cInUp(0.), dInUp(0.), eInUp(0.);
2398  Double_t aOutUp(0.), cOutUp(0.), dOutUp(0.), eOutUp(0.);
2399  Double_t aInLow(0.), cInLow(0.), dInLow(0.), eInLow(0.);
2400  Double_t aOutLow(0.), cOutLow(0.), dOutLow(0.), eOutLow(0.);
2401 
2402  GetErrorFromFit(relativeErrorRecBinInUp, relativeErrorRecBinInLow, rangeLow, rangeUp, fPivot, fSubdueError, "measured, in-plane");
2403  GetErrorFromFit(relativeErrorRecBinOutUp, relativeErrorRecBinOutLow, rangeLow, rangeUp, fPivot, fSubdueError, "measured, out-of-plane");
2404  GetErrorFromFit(relativeErrorMethodInUp, relativeErrorMethodInLow, rangeLow, rangeUp, fPivot, fSubdueError, "background, in-plane");
2405  GetErrorFromFit(relativeErrorMethodOutUp, relativeErrorMethodOutLow, rangeLow, rangeUp, fPivot, fSubdueError, "background, out-of-plane");
2406 
2407  for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2408  // for the upper bound only add regularization stuff if it is NOT done on v2 directly
2409  if(relativeErrorRegularizationInUp && !regularizationOnV2) aInUp = relativeErrorRegularizationInUp->GetBinContent(b+1);
2410  if(relativeErrorRegularizationOutUp && !regularizationOnV2) aOutUp = relativeErrorRegularizationOutUp->GetBinContent(b+1);
2411  if(relativeErrorRecBinInUp) cInUp = relativeErrorRecBinInUp->GetBinContent(b+1);
2412  if(relativeErrorRecBinOutUp) cOutUp = relativeErrorRecBinOutUp->GetBinContent(b+1);
2413  if(relativeErrorMethodInUp) dInUp = relativeErrorMethodInUp->GetBinContent(b+1);
2414  if(relativeErrorMethodOutUp) dOutUp = relativeErrorMethodOutUp->GetBinContent(b+1);
2415  // for the lower bound
2416  if(relativeErrorRegularizationInLow) aInLow = relativeErrorRegularizationInLow->GetBinContent(b+1);
2417  if(relativeErrorRegularizationOutLow) aOutLow = relativeErrorRegularizationOutLow->GetBinContent(b+1);
2418  if(relativeErrorRecBinInLow) cInLow = relativeErrorRecBinInLow->GetBinContent(b+1);
2419  if(relativeErrorRecBinOutLow) cOutLow = relativeErrorRecBinOutLow->GetBinContent(b+1);
2420  if(relativeErrorMethodInLow) dInLow = relativeErrorMethodInLow->GetBinContent(b+1);
2421  if(relativeErrorMethodOutLow) dOutLow = relativeErrorMethodOutLow->GetBinContent(b+1);
2422 
2423  if(fSymmRMS) { // since there is no a-priori assumption, we take these errors as symmetric
2424  RemoveSign(aInLow);
2425  RemoveSign(cInLow);
2426  RemoveSign(dInLow);
2427  RemoveSign(aOutLow);
2428  RemoveSign(cOutLow);
2429  RemoveSign(dOutLow);
2430  if(aInLow > aInUp) { aInUp = aInLow;
2431  } else {aInLow = aInUp;}
2432  if(aOutLow > aOutUp) { aOutUp = aOutLow;
2433  } else { aOutLow = aOutUp;}
2434  if(cInLow > cInUp) { cInUp = cInLow;
2435  } else {cInLow = cInUp; };
2436  if(cOutLow > cOutUp ) { cOutUp = cOutLow;
2437  } else { cOutLow = cOutUp; }
2438  // for confusio the logic changes here ...
2439  if(dInLow < dInUp) {dInLow = dInUp;
2440  } else { dInUp = dInLow;}
2441  if(dOutLow < dOutUp) { dOutLow = dOutUp;
2442  } else { dOutUp = dOutLow; }
2443  }
2444  // get the quadratic sums
2445  eInUp = aInUp*aInUp + cInUp*cInUp + dInUp*dInUp;
2446  if(eInUp > 0) relativeErrorInUp->SetBinContent(b+1, 1.*TMath::Sqrt(eInUp));
2447  eOutUp = aOutUp*aOutUp + cOutUp*cOutUp + dOutUp*dOutUp;
2448  if(eOutUp > 0) relativeErrorOutUp->SetBinContent(b+1, 1.*TMath::Sqrt(eOutUp));
2449 
2450  eInLow = aInLow*aInLow + cInLow*cInLow + dInLow*dInLow;
2451  if(eInLow > 0) relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(eInLow));
2452  eOutLow = aOutLow*aOutLow + cOutLow*cOutLow + dOutLow*dOutLow;
2453  if(eOutLow > 0) relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(eOutLow));
2454 
2455  printf("%i) \teInUp %.4f \t eInLow %.4f \t eOutUp %4.f \t eOutLow %4.f \n", b, eInUp, eInLow, eOutUp, eOutLow);
2456  }
2457 
2458  // project the estimated errors on the nominal ratio
2459  if(nominal) {
2460  Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
2461  Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
2462  Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
2463  Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
2464  Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
2465  Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
2466  Double_t lowErr(0.), upErr(0.);
2467  for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
2468  // add the in and out of plane errors in quadrature
2469  // take special care here: to propagate the assymetric error, we need to correlate the
2470  // InLow with OutUp (minimum value of ratio) and InUp with OutLow (maximum value of ratio)
2471  lowErr = relativeErrorInLow->GetBinContent(i+1)*relativeErrorInLow->GetBinContent(i+1)+relativeErrorOutUp->GetBinContent(1+i)*relativeErrorOutUp->GetBinContent(i+1);
2472  upErr = relativeErrorInUp->GetBinContent(i+1)*relativeErrorInUp->GetBinContent(i+1)+relativeErrorOutLow->GetBinContent(i+1)*relativeErrorOutLow->GetBinContent(i+1);
2473  // set the errors
2474  ayl[i] = TMath::Sqrt(lowErr)*nominal->GetBinContent(i+1);
2475  ayh[i] = TMath::Sqrt(upErr)*nominal->GetBinContent(i+1);
2476  // get the bin width (which is the 'error' on x
2477  Double_t binWidth(nominal->GetBinWidth(i+1));
2478  axl[i] = binWidth/2.;
2479  axh[i] = binWidth/2.;
2480  // now get the coordinate for the point
2481  ax[i] = nominal->GetBinCenter(i+1);
2482  ay[i] = nominal->GetBinContent(i+1);
2483  }
2484 
2485  // save the nominal ratio
2486  TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
2487  shapeRatio = (TGraphAsymmErrors*)nominalError->Clone();
2488  nominalError->SetName("shape uncertainty");
2489  TCanvas* nominalCanvas(new TCanvas("nominalCanvas", "nominalCanvas"));
2490  nominalCanvas->Divide(2);
2491  nominalCanvas->cd(1);
2492  Style(nominal, kBlack);
2493  Style(nominalError, kCyan, kRatio);
2494  nominalError->SetLineColor(kCyan);
2495  nominalError->SetMarkerColor(kCyan);
2496  nominalError->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2497  nominalError->GetYaxis()->SetRangeUser(.7, 2.2);
2498  nominalError->DrawClone("a2");
2499  nominal->DrawCopy("same E1");
2500  Style(AddLegend(gPad));
2501  Style(gPad, "GRID");
2502  Style(nominalCanvas);
2503  // save nominal v2 and systematic errors
2505  nominalIn,
2506  nominalOut,
2508  "v_{2} with shape uncertainty",
2509  relativeErrorInUp,
2510  relativeErrorInLow,
2511  relativeErrorOutUp,
2512  relativeErrorOutLow,
2513  corr));
2514  shapeV2 = (TGraphAsymmErrors*)nominalV2Error->Clone();
2515  if(regularizationOnV2) shapeV2 = AddHistoErrorsToGraphErrors(shapeV2, relativeErrorRegularizationInUp);
2516  TGraphErrors* nominalV2(GetV2(nominalIn, nominalOut, fEventPlaneRes, "v_{2}"));
2517  // add in quadratufe (not very nice, rethink this because it may add additional weight to
2518  // the rms unfolded piece) the additional error
2519  nominalCanvas->cd(2);
2520  Style(nominalV2, kBlack);
2521  Style(nominalV2Error, kCyan, kV2);
2522  nominalV2Error->SetLineColor(kCyan);
2523  nominalV2Error->SetMarkerColor(kCyan);
2524  nominalV2Error->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
2525  nominalV2Error->DrawClone("a2");
2526  nominalV2->DrawClone("same E1");
2527  Style(AddLegend(gPad));
2528  Style(gPad, "GRID");
2529  Style(nominalCanvas);
2530  SavePadToPDF(nominalCanvas);
2531  nominalCanvas->Write();
2532  }
2533  TCanvas* relativeError(new TCanvas("relativeError"," relativeError"));
2534  relativeError->Divide(2);
2535  relativeError->cd(1);
2536  Style(gPad, "GRID");
2537  relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
2538  Style(relativeErrorInUp, kBlue, kBar);
2539  Style(relativeErrorInLow, kGreen, kBar);
2540  relativeErrorInUp->DrawCopy("b");
2541  relativeErrorInLow->DrawCopy("same b");
2542  Style(relativeStatisticalErrorIn, kRed);
2543  relativeStatisticalErrorIn->DrawCopy("same");
2544  Style(AddLegend(gPad));
2545  relativeError->cd(2);
2546  Style(gPad, "GRID");
2547  Style(relativeErrorOutUp, kBlue, kBar);
2548  Style(relativeErrorOutLow, kGreen, kBar);
2549  relativeErrorOutUp->DrawCopy("b");
2550  if(relativeErrorOutLow) relativeErrorOutLow->DrawCopy("same b");
2551  if(relativeStatisticalErrorOut) Style(relativeStatisticalErrorOut, kRed);
2552  if(relativeStatisticalErrorOut) relativeStatisticalErrorOut->DrawCopy("same");
2553  Style(AddLegend(gPad));
2554 
2555  // write the buffered file to disk
2556  SavePadToPDF(relativeError);
2557  relativeError->Write();
2558  output->Write();
2559 }
2560 //_____________________________________________________________________________
2562  TArrayI* variationsIn, // variantions in plane
2563  TArrayI* variationsOut, // variantions out of plane
2564  TH1D*& relativeErrorInUp, // pointer reference to minimum relative error histogram in plane
2565  TH1D*& relativeErrorInLow, // pointer reference to maximum relative error histogram in plane
2566  TH1D*& relativeErrorOutUp, // pointer reference to minimum relative error histogram out of plane
2567  TH1D*& relativeErrorOutLow, // pointer reference to maximum relative error histogram out of plane
2568  TH1D*& relativeStatisticalErrorIn, // relative systematic error on ratio
2569  TH1D*& relativeStatisticalErrorOut, // relative systematic error on ratio
2570  TH1D*& nominal, // clone of the nominal ratio in / out of plane
2571  TH1D*& nominalIn, // clone of the nominal in plane yield
2572  TH1D*& nominalOut, // clone of the nominal out of plane yield
2573  Int_t columns, // divide the output canvasses in this many columns
2574  Float_t rangeLow, // lower pt range
2575  Float_t rangeUp, // upper pt range
2576  TFile* readMe, // input file name (created by this unfolding class)
2577  TString source, // source of the variation
2578  Bool_t RMS, // return RMS of distribution of variations as error
2579  Bool_t onRatio // use ratio or v2 directly for assessing systematic uncertainty
2580  ) const
2581 {
2582 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
2583  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2584 #endif
2585  // wrapper function to call specific systematics function using function pointers
2586  void (AliJetFlowTools::*myFunction)(TArrayI*, TArrayI*, TH1D*&, TH1D*&, TH1D*&, TH1D*&, TH1D*&, TH1D*&,
2587  TH1D*&, TH1D*&, TH1D*&, Int_t, Float_t, Float_t, TFile*, TString, Bool_t) const;
2588 
2589  printf(" >>> systematic wrapper called for %s <<< \n", source.Data());
2590 
2591  // assign function pointer
2593 
2594  // do the actual unfolding with the selected function
2595  return (this->*myFunction)(
2596  variationsIn,
2597  variationsOut,
2598  relativeErrorInUp,
2599  relativeErrorInLow,
2600  relativeErrorOutUp,
2601  relativeErrorOutLow,
2602  relativeStatisticalErrorIn,
2603  relativeStatisticalErrorOut,
2604  nominal,
2605  nominalIn,
2606  nominalOut,
2607  columns,
2608  rangeLow,
2609  rangeUp,
2610  readMe,
2611  source,
2612  RMS);
2613 }
2614 //_____________________________________________________________________________
2616  TArrayI* variationsIn, // variantions in plane
2617  TArrayI* variationsOut, // variantions out of plane
2618  TH1D*& relativeErrorInUp, // pointer reference to minimum relative error histogram in plane
2619  TH1D*& relativeErrorInLow, // pointer reference to maximum relative error histogram in plane
2620  TH1D*& relativeErrorOutUp, // pointer reference to minimum relative error histogram out of plane
2621  TH1D*& relativeErrorOutLow, // pointer reference to maximum relative error histogram out of plane
2622  TH1D*& relativeStatisticalErrorIn, // relative systematic error on ratio
2623  TH1D*& relativeStatisticalErrorOut, // relative systematic error on ratio
2624  TH1D*& nominal, // clone of the nominal ratio in / out of plane
2625  TH1D*& nominalIn, // clone of the nominal in plane yield
2626  TH1D*& nominalOut, // clone of the nominal out of plane yield
2627  Int_t columns, // divide the output canvasses in this many columns
2628  Float_t rangeLow, // lower pt range
2629  Float_t rangeUp, // upper pt range
2630  TFile* readMe, // input file name (created by this unfolding class)
2631  TString source, // source of the variation
2632  Bool_t RMS // return RMS of distribution of variations as error
2633  ) const
2634 {
2635 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
2636  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
2637 #endif
2638  // intermediate systematic check function. first index of supplied array is nominal value
2639  TList* listOfKeys((TList*)readMe->GetListOfKeys());
2640  if(!listOfKeys) {
2641  printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
2642  return;
2643  }
2644  // check input params
2645  if(variationsIn->GetSize() != variationsOut->GetSize()) {
2646  printf(" > DoSystematics: fatal error, input arrays have different sizes ! < \n ");
2647  return;
2648  }
2649  TDirectoryFile* defRootDirIn(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(0))->GetName())));
2650  TDirectoryFile* defRootDirOut(dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(0))->GetName())));
2651  if(!(defRootDirIn && defRootDirOut)) {
2652  printf(" > DoSystematics: fatal error, couldn't retrieve nominal values ! < \n ");
2653  return;
2654  }
2655  TString defIn(defRootDirIn->GetName());
2656  TString defOut(defRootDirOut->GetName());
2657 
2658  // define lines to make the output prettier
2659  TLine* lineLow(new TLine(rangeLow, 0., rangeLow, 2.));
2660  TLine* lineUp(new TLine(rangeUp, 0., rangeUp, 2.));
2661  lineLow->SetLineColor(11);
2662  lineUp->SetLineColor(11);
2663  lineLow->SetLineWidth(3);
2664  lineUp->SetLineWidth(3);
2665 
2666  // define an output histogram with the maximum relative error from this systematic constribution
2667  // 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
2668  // reached in this function call.
2669  // 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
2670  // which should correspond to a 68% confidence level
2671  relativeErrorInUp = new TH1D(Form("relative error (up) from %s", source.Data()), Form("relative error (up) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2672  relativeErrorInLow = new TH1D(Form("relative error (low) from %s", source.Data()), Form("relative error (low) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2673  relativeErrorOutUp = new TH1D(Form("relative error (up) from %s", source.Data()), Form("relative error (up) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2674  relativeErrorOutLow = new TH1D(Form("relative error (low) from %s", source.Data()), Form("relative error (low) from %s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2675  for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2676  if(!RMS) {
2677  relativeErrorInUp->SetBinContent(b+1, 1.);
2678  relativeErrorInUp->SetBinError(b+1, 0.);
2679  relativeErrorOutUp->SetBinContent(b+1, 1.);
2680  relativeErrorOutUp->SetBinError(b+1, .0);
2681  relativeErrorInLow->SetBinContent(b+1, 1.);
2682  relativeErrorInLow->SetBinError(b+1, 0.);
2683  relativeErrorOutLow->SetBinContent(b+1, 1.);
2684  relativeErrorOutLow->SetBinError(b+1, .0);
2685  } else if(RMS) {
2686  relativeErrorInUp->SetBinContent(b+1, 0.);
2687  relativeErrorInUp->SetBinError(b+1, 0.);
2688  relativeErrorOutUp->SetBinContent(b+1, 0.);
2689  relativeErrorOutUp->SetBinError(b+1, 0.);
2690  relativeErrorInLow->SetBinContent(b+1, 0.);
2691  relativeErrorInLow->SetBinError(b+1, 0.);
2692  relativeErrorOutLow->SetBinContent(b+1, 0.);
2693  relativeErrorOutLow->SetBinError(b+1, 0.);
2694  }
2695  }
2696  Int_t relativeErrorInUpN[100] = {0};
2697  Int_t relativeErrorOutUpN[100] = {0};
2698  Int_t relativeErrorInLowN[100] = {0};
2699  Int_t relativeErrorOutLowN[100] = {0};
2700 
2701  // define an output histogram with the systematic error from this systematic constribution
2702  if(!relativeStatisticalErrorIn || !relativeStatisticalErrorOut) {
2703  relativeStatisticalErrorIn = new TH1D("relative statistical error, in plane", "relative statistical error, in plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2704  relativeStatisticalErrorOut = new TH1D("relative statistical error, out of plane", "relative statistital error, out of plane", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
2705  }
2706 
2707  // prepare necessary canvasses
2708  TCanvas* canvasIn(new TCanvas(Form("SYST_%s_PearsonIn", source.Data()), Form("SYST_%s_PearsonIn", source.Data())));
2709  TCanvas* canvasOut(0x0);
2710  if(fDphiUnfolding) canvasOut = new TCanvas(Form("SYST_%s_PearsonOut", source.Data()), Form("SYST_%s_PearsonOut", source.Data()));
2711  TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas(Form("SYST_%s_RefoldedIn", source.Data()), Form("SYST_%s_RefoldedIn", source.Data())));
2712  TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
2713  if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas(Form("SYST_%s_RefoldedOut", source.Data()), Form("SYST_%s_RefoldedOut", source.Data()));
2714  TCanvas* canvasSpectraIn(new TCanvas(Form("SYST_%s_SpectraIn", source.Data()), Form("SYST_%s_SpectraIn", source.Data())));
2715  TCanvas* canvasSpectraOut(0x0);
2716  if(fDphiUnfolding) canvasSpectraOut = new TCanvas(Form("SYST_%s_SpectraOut", source.Data()), Form("SYST_%s_SpectraOut", source.Data()));
2717  TCanvas* canvasRatio(0x0);
2718  if(fDphiUnfolding) canvasRatio = new TCanvas(Form("SYST_%s_Ratio", source.Data()), Form("SYST_%s_Ratio", source.Data()));
2719  TCanvas* canvasV2(0x0);
2720  if(fDphiUnfolding) canvasV2 = new TCanvas(Form("SYST_%s_V2", source.Data()), Form("SYST_%s_V2", source.Data()));
2721  TCanvas* canvasMISC(new TCanvas(Form("SYST_%s_MISC", source.Data()), Form("SYST_%s_MISC", source.Data())));
2722  TCanvas* canvasMasterIn(new TCanvas(Form("SYST_%s_defaultIn", source.Data()), Form("SYST_%s_defaultIn", source.Data())));
2723  TCanvas* canvasMasterOut(0x0);
2724  if(fDphiUnfolding) canvasMasterOut = new TCanvas(Form("SYST_%s_defaultOut", source.Data()), Form("SYST_%s_defaultOut", source.Data()));
2725  (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
2726 
2727  TCanvas* canvasProfiles(new TCanvas(Form("SYST_%s_canvasProfiles", source.Data()), Form("SYST_%s_canvasProfiles", source.Data())));
2728  canvasProfiles->Divide(2);
2729  TProfile* ratioProfile(new TProfile(Form("SYST_%s_ratioProfile", source.Data()), Form("SYST_%s_ratioProfile", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2730  TProfile* v2Profile(new TProfile(Form("SYST_%s_v2Profile", source.Data()), Form("SYST_%s_v2Profile", source.Data()),fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
2731  // get an estimate of the number of outputs and find the default set
2732 
2733  Int_t rows = 1;
2734  columns = variationsIn->GetSize()-1;
2735  (TMath::Floor(variationsIn->GetSize()/(float)columns)+((variationsIn->GetSize()%columns)>0));
2736  canvasIn->Divide(columns, rows);
2737  if(canvasOut) canvasOut->Divide(columns, rows);
2738  canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
2739  if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
2740  canvasSpectraIn->Divide(columns, rows);
2741  if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
2742  if(canvasRatio) canvasRatio->Divide(columns, rows);
2743  if(canvasV2) canvasV2->Divide(columns, rows);
2744  canvasMasterIn->Divide(columns, rows);
2745  if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
2746 
2747  // prepare a separate set of canvases to hold the nominal points
2748  TCanvas* canvasNominalIn(new TCanvas(Form("Nominal_%s_PearsonIn", source.Data()), Form("Nominal_%s_PearsonIn", source.Data())));
2749  TCanvas* canvasNominalOut(0x0);
2750  if(fDphiUnfolding) canvasNominalOut = new TCanvas(Form("Nominal_%s_PearsonOut", source.Data()), Form("Nominal_%s_PearsonOut", source.Data()));
2751  TCanvas* canvasNominalRatioMeasuredRefoldedIn(new TCanvas(Form("Nominal_%s_RefoldedIn", source.Data()), Form("Nominal_%s_RefoldedIn", source.Data())));
2752  TCanvas* canvasNominalRatioMeasuredRefoldedOut(0x0);
2753  if(fDphiUnfolding) canvasNominalRatioMeasuredRefoldedOut = new TCanvas(Form("Nominal_%s_RefoldedOut", source.Data()), Form("Nominal_%s_RefoldedOut", source.Data()));
2754  TCanvas* canvasNominalSpectraIn(new TCanvas(Form("Nominal_%s_SpectraIn", source.Data()), Form("Nominal_%s_SpectraIn", source.Data())));
2755  TCanvas* canvasNominalSpectraOut(0x0);
2756  if(fDphiUnfolding) canvasNominalSpectraOut = new TCanvas(Form("Nominal_%s_SpectraOut", source.Data()), Form("Nominal_%s_SpectraOut", source.Data()));
2757  TCanvas* canvasNominalRatio(0x0);
2758  if(fDphiUnfolding) canvasNominalRatio = new TCanvas(Form("Nominal_%s_Ratio", source.Data()), Form("Nominal_%s_Ratio", source.Data()));
2759  TCanvas* canvasNominalV2(0x0);
2760  if(fDphiUnfolding) canvasNominalV2 = new TCanvas(Form("Nominal_%s_V2", source.Data()), Form("Nominal_%s_V2", source.Data()));
2761  TCanvas* canvasNominalMISC(new TCanvas(Form("Nominal_%s_MISC", source.Data()), Form("Nominal_%s_MISC", source.Data())));
2762  TCanvas* canvasNominalMasterIn(new TCanvas(Form("Nominal_%s_defaultIn", source.Data()), Form("Nominal_%s_defaultIn", source.Data())));
2763  TCanvas* canvasNominalMasterOut(0x0);
2764  if(fDphiUnfolding) canvasNominalMasterOut = new TCanvas(Form("Nominal_%s_defaultOut", source.Data()), Form("Nominal_%s_defaultOut", source.Data()));
2765  (fDphiUnfolding) ? canvasNominalMISC->Divide(4, 2) : canvasNominalMISC->Divide(4, 1);
2766 
2767  canvasNominalSpectraIn->Divide(2);
2768  if(canvasNominalSpectraOut) canvasNominalSpectraOut->Divide(2);
2769 
2770  canvasNominalMasterIn->Divide(2);
2771  if(canvasNominalMasterOut) canvasNominalMasterOut->Divide(2);
2772 
2773  // extract the default output
2774  TH1D* defaultUnfoldedJetSpectrumIn(0x0);
2775  TH1D* defaultUnfoldedJetSpectrumOut(0x0);
2776  TDirectoryFile* defDirIn = (TDirectoryFile*)defRootDirIn->Get(Form("InPlane___%s", defIn.Data()));
2777  TDirectoryFile* defDirOut = (TDirectoryFile*)defRootDirOut->Get(Form("OutOfPlane___%s", defOut.Data()));
2778  if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", defIn.Data()));
2779  if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", defOut.Data()));
2780  printf(" > succesfully extracted default results < \n");
2781 
2782  // loop through the directories, only plot the graphs if the deconvolution converged
2783  TDirectoryFile* tempDirIn(0x0);
2784  TDirectoryFile* tempDirOut(0x0);
2785  TDirectoryFile* tempIn(0x0);
2786  TDirectoryFile* tempOut(0x0);
2787  TH1D* unfoldedSpectrumInForRatio(0x0);
2788  TH1D* unfoldedSpectrumOutForRatio(0x0);
2789  for(Int_t i(0), j(-1); i < variationsIn->GetSize(); i++) {
2790  tempDirIn = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsIn->At(i))->GetName())));
2791  tempDirOut = (dynamic_cast<TDirectoryFile*>(readMe->Get(listOfKeys->At(variationsOut->At(i))->GetName())));
2792  if(!(tempDirIn && tempDirOut)) {
2793  printf(" > DoSystematics: couldn't get a set of variations < \n");
2794  continue;
2795  }
2796  TString dirNameIn(tempDirIn->GetName());
2797  TString dirNameOut(tempDirOut->GetName());
2798  // try to read the in- and out of plane subdirs
2799  tempIn = (TDirectoryFile*)tempDirIn->Get(Form("InPlane___%s", dirNameIn.Data()));
2800  tempOut = (TDirectoryFile*)tempDirOut->Get(Form("OutOfPlane___%s", dirNameOut.Data()));
2801  j++;
2802  if(tempIn) {
2803  // to see if the unfolding converged try to extract the pearson coefficients
2804  TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirNameIn.Data())));
2805  if(pIn) {
2806  printf(" - %s in plane converged \n", dirNameIn.Data());
2807  canvasIn->cd(j);
2808  if(i==0) canvasNominalIn->cd(j);
2809  Style(gPad, "PEARSON");
2810  pIn->DrawCopy("colz");
2811  TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirNameIn.Data())));
2812  if(rIn) {
2813  Style(rIn);
2814  printf(" > found RatioRefoldedMeasured < \n");
2815  canvasRatioMeasuredRefoldedIn->cd(j);
2816  if(i==0) canvasNominalRatioMeasuredRefoldedIn->cd(j);
2817  Style(gPad, "GRID");
2818  rIn->SetFillColor(kRed);
2819  rIn->Draw("ap");
2820  }
2821  TH1D* dvector((TH1D*)tempIn->Get("dVector"));
2822  TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
2823  TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirNameIn.Data())));
2824  TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirNameIn.Data())));
2825  if(dvector && avalue && rm && eff) {
2826  Style(dvector);
2827  Style(avalue);
2828  Style(rm);
2829  Style(eff);
2830  canvasMISC->cd(1);
2831  if(i==0) canvasNominalMISC->cd(1);
2832  Style(gPad, "SPECTRUM");
2833  dvector->DrawCopy();
2834  canvasMISC->cd(2);
2835  if(i==0) canvasNominalMISC->cd(2);
2836  Style(gPad, "SPECTRUM");
2837  avalue->DrawCopy();
2838  canvasMISC->cd(3);
2839  if(i==0) canvasNominalMISC->cd(3);
2840  Style(gPad, "PEARSON");
2841  rm->DrawCopy("colz");
2842  canvasMISC->cd(4);
2843  if(i==0) canvasNominalMISC->cd(4);
2844  Style(gPad, "GRID");
2845  eff->DrawCopy();
2846  } else if(rm && eff) {
2847  Style(rm);
2848  Style(eff);
2849  canvasMISC->cd(3);
2850  if(i==0) canvasNominalMISC->cd(3);
2851  Style(gPad, "PEARSON");
2852  rm->DrawCopy("colz");
2853  canvasMISC->cd(4);
2854  if(i==0) canvasNominalMISC->cd(4);
2855  Style(gPad, "GRID");
2856  eff->DrawCopy();
2857  }
2858  }
2859  TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirNameIn.Data())));
2860  TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirNameIn.Data())));
2861  unfoldedSpectrumInForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
2862  TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirNameIn.Data())));
2863  if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
2864  if(defaultUnfoldedJetSpectrumIn) {
2865  Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
2866  TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirNameIn.Data())));
2867  temp->Divide(unfoldedSpectrum);
2868  // get the absolute relative error
2869  for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
2870  if(!RMS) { // save the maximum deviation that a variation can cause
2871  // the variation is HIGHER than the nominal point, so the bar goes UP
2872  if( temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorInUp->GetBinContent(b+1)) {
2873  relativeErrorInUp->SetBinContent(b+1, temp->GetBinContent(b+1));
2874  relativeErrorInUp->SetBinError(b+1, temp->GetBinError(b+1));
2875  }
2876  // the variation is LOWER than the nominal point, so the bar goes DOWN
2877  else if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorInLow->GetBinContent(b+1)) {
2878  relativeErrorInLow->SetBinContent(b+1, temp->GetBinContent(b+1));
2879  relativeErrorInLow->SetBinError(b+1, temp->GetBinError(b+1));
2880  }
2881  } else if (RMS && !fSymmRMS) { // save info necessary for evaluating the RMS of a distribution of variations
2882  printf(" oops shouldnt be here \n " );
2883  if(temp->GetBinContent(b+1) < 1) {
2884  relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2885  relativeErrorInUpN[b]++;
2886  }
2887  // the variation is LOWER than the nominal point, so the bar goes DOWN
2888  else if(temp->GetBinContent(b+1) > 1) {
2889  relativeErrorInLow->SetBinContent(b+1, relativeErrorInLow->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
2890  relativeErrorInLowN[b]++;
2891  }
2892  } else if (fSymmRMS) {
2893  // save symmetric sum of square to get a symmetric rms
2894  relativeErrorInUp->SetBinContent(b+1, relativeErrorInUp->GetBinContent(b+1)+TMath::Power(temp->GetBinContent(b+1)-1., 2));
2895  relativeErrorInUpN[b]++;
2896  }
2897  if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorIn->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
2898  }
2899  temp->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2900  temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
2901  temp->GetYaxis()->SetTitle("ratio");
2902  canvasMasterIn->cd(j);
2903  temp->GetYaxis()->SetRangeUser(0., 2);
2904  Style(gPad, "GRID");
2905  temp->DrawCopy();
2906  canvasNominalMasterIn->cd(1);
2907  Style(gPad, "GRID");
2908  if(i > 0 ) {
2909  TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
2910  tempSyst->SetTitle(Form("[%s] / [%s]", defIn.Data(), dirNameIn.Data()));
2911  Style(tempSyst, (EColor)(i+2));
2912  if(i==1) tempSyst->DrawCopy();
2913  else tempSyst->DrawCopy("same");
2914  }
2915  }
2916  TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirNameIn.Data())));
2917  canvasSpectraIn->cd(j);
2918  if(i==0) canvasNominalSpectraIn->cd(1);
2919  Style(gPad);
2920  Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
2921  unfoldedSpectrum->DrawCopy();
2922  Style(inputSpectrum, kGreen, kMeasuredSpectrum);
2923  inputSpectrum->DrawCopy("same");
2924  Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
2925  refoldedSpectrum->DrawCopy("same");
2926  TLegend* l(AddLegend(gPad));
2927  Style(l);
2928  if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
2929  Float_t chi(fitStatus->GetBinContent(1));
2930  Float_t pen(fitStatus->GetBinContent(2));
2931  Int_t dof((int)fitStatus->GetBinContent(3));
2932  Float_t beta(fitStatus->GetBinContent(4));
2933  l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
2934  } else if (fitStatus) { // only available in SVD fit
2935  Int_t reg((int)fitStatus->GetBinContent(1));
2936  l->AddEntry((TObject*)0, Form("REG %i", reg), "");
2937  }
2938  canvasNominalSpectraIn->cd(2);
2939  TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
2940  tempSyst->SetTitle(Form("[%s]", dirNameIn.Data()));
2941  Style(tempSyst, (EColor)(i+2));
2942  Style(gPad, "SPECTRUM");
2943  if(i==0) tempSyst->DrawCopy();
2944  else tempSyst->DrawCopy("same");
2945  }
2946  }
2947  if(tempOut) {
2948  TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirNameOut.Data())));
2949  if(pOut) {
2950  printf(" - %s out of plane converged \n", dirNameOut.Data());
2951  canvasOut->cd(j);
2952  if(i==0) canvasNominalOut->cd(j);
2953  Style(gPad, "PEARSON");
2954  pOut->DrawCopy("colz");
2955  TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirNameOut.Data())));
2956  if(rOut) {
2957  Style(rOut);
2958  printf(" > found RatioRefoldedMeasured < \n");
2959  canvasRatioMeasuredRefoldedOut->cd(j);
2960  if(i==0) canvasNominalRatioMeasuredRefoldedOut->cd(j);
2961  Style(gPad, "GRID");
2962  rOut->SetFillColor(kRed);
2963  rOut->Draw("ap");
2964  }
2965  TH1D* dvector((TH1D*)tempOut->Get("dVector"));
2966  TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
2967  TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirNameOut.Data())));
2968  TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirNameOut.Data())));
2969  if(dvector && avalue && rm && eff) {
2970  Style(dvector);
2971  Style(avalue);
2972  Style(rm);
2973  Style(eff);
2974  canvasMISC->cd(5);
2975  if(i==0) canvasNominalMISC->cd(5);
2976  Style(gPad, "SPECTRUM");
2977  dvector->DrawCopy();
2978  canvasMISC->cd(6);
2979  if(i==0) canvasNominalMISC->cd(6);
2980  Style(gPad, "SPECTRUM");
2981  avalue->DrawCopy();
2982  canvasMISC->cd(7);
2983  if(i==0) canvasNominalMISC->cd(7);
2984  Style(gPad, "PEARSON");
2985  rm->DrawCopy("colz");
2986  canvasMISC->cd(8);
2987  if(i==0) canvasNominalMISC->cd(8);
2988  Style(gPad, "GRID");
2989  eff->DrawCopy();
2990  } else if(rm && eff) {
2991  Style(rm);
2992  Style(eff);
2993  canvasMISC->cd(7);
2994  if(i==0) canvasNominalMISC->cd(7);
2995  Style(gPad, "PEARSON");
2996  rm->DrawCopy("colz");
2997  canvasMISC->cd(8);
2998  if(i==0) canvasNominalMISC->cd(8);
2999  Style(gPad, "GRID");
3000  eff->DrawCopy();
3001  }
3002  }
3003  TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirNameOut.Data())));
3004  TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirNameOut.Data())));
3005  unfoldedSpectrumOutForRatio = ProtectHeap(unfoldedSpectrum, kFALSE, TString("ForRatio"));
3006  TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirNameOut.Data())));
3007  if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3008  if(defaultUnfoldedJetSpectrumOut) {
3009  Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
3010  TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirNameOut.Data())));
3011  temp->Divide(unfoldedSpectrum);
3012  // get the absolute relative error
3013  for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
3014  if(!RMS) {
3015  // check if the error is larger than the current maximum
3016  if(temp->GetBinContent(b+1) < 1 && temp->GetBinContent(b+1) < relativeErrorOutUp->GetBinContent(b+1)) {
3017  relativeErrorOutUp->SetBinContent(b+1, temp->GetBinContent(b+1));
3018  relativeErrorOutUp->SetBinError(b+1, temp->GetBinError(b+1));
3019  }
3020  // check if the error is smaller than the current minimum
3021  else if(temp->GetBinContent(b+1) > 1 && temp->GetBinContent(b+1) > relativeErrorOutLow->GetBinContent(b+1)) {
3022  relativeErrorOutLow->SetBinContent(b+1, temp->GetBinContent(b+1));
3023  relativeErrorOutLow->SetBinError(b+1, temp->GetBinError(b+1));
3024  }
3025  } else if (RMS && !fSymmRMS) {
3026  printf(" OOps \n ");
3027  if(temp->GetBinContent(b+1) < 1) {
3028  relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
3029  relativeErrorOutUpN[b]++;
3030  }
3031  else if(temp->GetBinContent(b+1) > 1) {
3032  relativeErrorOutLow->SetBinContent(b+1, relativeErrorOutLow->GetBinContent(b+1)+TMath::Power(1.-temp->GetBinContent(b+1), 2));
3033  relativeErrorOutLowN[b]++;
3034  }
3035  } else if (fSymmRMS) {
3036  // save symmetric rms value
3037  relativeErrorOutUp->SetBinContent(b+1, relativeErrorOutUp->GetBinContent(b+1)+TMath::Power(temp->GetBinContent(b+1)-1., 2));
3038  relativeErrorOutUpN[b]++;
3039  }
3040  if(temp->GetBinError(b+1) > 0) relativeStatisticalErrorOut->SetBinContent(b+1, temp->GetBinError(b+1)/temp->GetBinContent(b+1));
3041  }
3042  temp->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
3043  temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3044  temp->GetYaxis()->SetTitle("ratio");
3045  canvasMasterOut->cd(j);
3046  temp->GetYaxis()->SetRangeUser(0., 2);
3047  Style(gPad, "GRID");
3048  temp->DrawCopy();
3049  canvasNominalMasterOut->cd(1);
3050  Style(gPad, "GRID");
3051  if(i > 0 ) {
3052  TH1D* tempSyst((TH1D*)temp->Clone(Form("%s_syst", temp->GetName())));
3053  tempSyst->SetTitle(Form("[%s] / [%s]", defOut.Data(), dirNameOut.Data()));
3054  Style(tempSyst, (EColor)(i+2));
3055  if(i==1) tempSyst->DrawCopy();
3056  else tempSyst->DrawCopy("same");
3057  }
3058  }
3059  TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirNameOut.Data())));
3060  canvasSpectraOut->cd(j);
3061  if(i==0) canvasNominalSpectraOut->cd(1);
3062  Style(gPad);
3063  Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3064  unfoldedSpectrum->DrawCopy();
3065  Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3066  inputSpectrum->DrawCopy("same");
3067  Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3068  refoldedSpectrum->DrawCopy("same");
3069  TLegend* l(AddLegend(gPad));
3070  Style(l);
3071  if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3072  Float_t chi(fitStatus->GetBinContent(1));
3073  Float_t pen(fitStatus->GetBinContent(2));
3074  Int_t dof((int)fitStatus->GetBinContent(3));
3075  Float_t beta(fitStatus->GetBinContent(4));
3076  l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3077  } else if (fitStatus) { // only available in SVD fit
3078  Int_t reg((int)fitStatus->GetBinContent(1));
3079  l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3080  }
3081  canvasNominalSpectraOut->cd(2);
3082  TH1D* tempSyst((TH1D*)unfoldedSpectrum->Clone(Form("%s_syst", unfoldedSpectrum->GetName())));
3083  tempSyst->SetTitle(Form("[%s]", dirNameOut.Data()));
3084  Style(tempSyst, (EColor)(i+2));
3085  Style(gPad, "SPECTRUM");
3086  if(i==0) tempSyst->DrawCopy();
3087  else tempSyst->DrawCopy("same");
3088  }
3089  }
3090  if(canvasRatio && canvasV2) {
3091  canvasRatio->cd(j);
3092  if(i==0) {
3093  canvasNominalRatio->cd(j);
3094  if(nominal && nominalIn && nominalOut) {
3095  // if a nominal ratio is requested, delete the placeholder and update the nominal point
3096  delete nominal;
3097  delete nominalIn;
3098  delete nominalOut;
3099  nominalIn = (TH1D*)unfoldedSpectrumInForRatio->Clone("in plane jet yield");
3100  nominalOut = (TH1D*)unfoldedSpectrumOutForRatio->Clone("out of plane jet yield");
3101  nominal = (TH1D*)unfoldedSpectrumInForRatio->Clone("ratio in plane out of plane");
3102  nominal->Divide(nominal, unfoldedSpectrumOutForRatio); // standard root divide for uncorrelated histos
3103  }
3104  }
3105  TGraphErrors* ratioYield(GetRatio(unfoldedSpectrumInForRatio, unfoldedSpectrumOutForRatio, TString(Form("ratio [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
3106  Double_t _tempx(0), _tempy(0);
3107  if(ratioYield) {
3108  Style(ratioYield);
3109  for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3110  ratioYield->GetPoint(b, _tempx, _tempy);
3111  ratioProfile->Fill(_tempx, _tempy);
3112  }
3113  ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
3114  ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3115  ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
3116  ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3117  ratioYield->SetFillColor(kRed);
3118  ratioYield->Draw("ap");
3119  }
3120  canvasV2->cd(j);
3121  if(i==0) canvasNominalV2->cd(j);
3122  TGraphErrors* ratioV2(GetV2(unfoldedSpectrumInForRatio,unfoldedSpectrumOutForRatio, fEventPlaneRes, TString(Form("v_{2} [in=%s, out=%s]", dirNameIn.Data(), dirNameOut.Data()))));
3123  if(ratioV2) {
3124  Style(ratioV2);
3125  for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3126  ratioV2->GetPoint(b, _tempx, _tempy);
3127  v2Profile->Fill(_tempx, _tempy);
3128 
3129  }
3130  v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
3131  v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3132  ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
3133  ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3134  ratioV2->SetFillColor(kRed);
3135  ratioV2->Draw("ap");
3136  }
3137  }
3138  delete unfoldedSpectrumInForRatio;
3139  delete unfoldedSpectrumOutForRatio;
3140  }
3141  // save the canvasses
3142  canvasProfiles->cd(1);
3143  Style(ratioProfile);
3144  ratioProfile->DrawCopy();
3145  canvasProfiles->cd(2);
3146  Style(v2Profile);
3147  v2Profile->DrawCopy();
3148  SavePadToPDF(canvasProfiles);
3149  canvasProfiles->Write();
3150  SavePadToPDF(canvasIn);
3151  canvasIn->Write();
3152  if(canvasOut) {
3153  SavePadToPDF(canvasOut);
3154  canvasOut->Write();
3155  }
3156  SavePadToPDF(canvasRatioMeasuredRefoldedIn);
3157  canvasRatioMeasuredRefoldedIn->Write();
3158  if(canvasRatioMeasuredRefoldedOut) {
3159  SavePadToPDF(canvasRatioMeasuredRefoldedOut);
3160  canvasRatioMeasuredRefoldedOut->Write();
3161  }
3162  SavePadToPDF(canvasSpectraIn);
3163  canvasSpectraIn->Write();
3164  if(canvasSpectraOut) {
3165  SavePadToPDF(canvasSpectraOut);
3166  canvasSpectraOut->Write();
3167  SavePadToPDF(canvasRatio);
3168  canvasRatio->Write();
3169  SavePadToPDF(canvasV2);
3170  canvasV2->Write();
3171  }
3172  SavePadToPDF(canvasMasterIn);
3173  canvasMasterIn->Write();
3174  if(canvasMasterOut) {
3175  SavePadToPDF(canvasMasterOut);
3176  canvasMasterOut->Write();
3177  }
3178  SavePadToPDF(canvasMISC);
3179  canvasMISC->Write();
3180  // save the nomial canvasses
3181  SavePadToPDF(canvasNominalIn);
3182  canvasNominalIn->Write();
3183  if(canvasNominalOut) {
3184  SavePadToPDF(canvasNominalOut);
3185  canvasNominalOut->Write();
3186  }
3187  SavePadToPDF(canvasNominalRatioMeasuredRefoldedIn);
3188  canvasNominalRatioMeasuredRefoldedIn->Write();
3189  if(canvasNominalRatioMeasuredRefoldedOut) {
3190  SavePadToPDF(canvasNominalRatioMeasuredRefoldedOut);
3191  canvasNominalRatioMeasuredRefoldedOut->Write();
3192  }
3193  canvasNominalSpectraIn->cd(2);
3194  Style(AddLegend(gPad));
3195  SavePadToPDF(canvasNominalSpectraIn);
3196  canvasNominalSpectraIn->Write();
3197  if(canvasNominalSpectraOut) {
3198  canvasNominalSpectraOut->cd(2);
3199  Style(AddLegend(gPad));
3200  SavePadToPDF(canvasNominalSpectraOut);
3201  canvasNominalSpectraOut->Write();
3202  SavePadToPDF(canvasNominalRatio);
3203  canvasNominalRatio->Write();
3204  SavePadToPDF(canvasNominalV2);
3205  canvasNominalV2->Write();
3206  }
3207  canvasNominalMasterIn->cd(1);
3208  Style(AddLegend(gPad));
3209  lineUp->DrawClone("same");
3210  lineLow->DrawClone("same");
3211  SavePadToPDF(canvasNominalMasterIn);
3212  canvasNominalMasterIn->Write();
3213  if(canvasNominalMasterOut) {
3214  canvasNominalMasterOut->cd(1);
3215  Style(AddLegend(gPad));
3216  lineUp->DrawClone("same");
3217  lineLow->DrawClone("same");
3218  SavePadToPDF(canvasNominalMasterOut);
3219  canvasNominalMasterOut->Write();
3220  }
3221  SavePadToPDF(canvasNominalMISC);
3222  canvasNominalMISC->Write();
3223 
3224  // save the relative errors
3225  for(Int_t b(0); b < fBinsTrue->GetSize()-1; b++) {
3226  // to arrive at a min and max from here, combine in up and out low
3227  if(!RMS) {
3228  relativeErrorInUp->SetBinContent(b+1, -1.*(relativeErrorInUp->GetBinContent(b+1)-1));
3229 // relativeErrorInUp->SetBinError(b+1, 0.);
3230  relativeErrorOutUp->SetBinContent(b+1, -1.*(relativeErrorOutUp->GetBinContent(b+1)-1));
3231 // relativeErrorOutUp->SetBinError(b+1, .0);
3232  relativeErrorInLow->SetBinContent(b+1, -1.*(relativeErrorInLow->GetBinContent(b+1)-1));
3233 // relativeErrorInLow->SetBinError(b+1, 0.);
3234  relativeErrorOutLow->SetBinContent(b+1, -1.*(relativeErrorOutLow->GetBinContent(b+1)-1));
3235 // relativeErrorOutLow->SetBinError(b+1, .0);
3236  } else if (RMS) {
3237  // these guys are already stored as percentages, so no need to remove the offset of 1
3238  // RMS is defined as sqrt(sum(squared))/N
3239  // min is defined as negative, max is defined as positive
3240  if(!fSymmRMS) {
3241  if(relativeErrorInUpN[b] < 1) relativeErrorInUpN[b] = 1;
3242  if(relativeErrorInLowN[b] < 1) relativeErrorInLowN[b] = 1;
3243  if(relativeErrorOutUpN[b] < 1) relativeErrorOutUpN[b] = 1;
3244  if(relativeErrorOutLowN[b] < 1) relativeErrorOutLowN[b] = 1;
3245  relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(b+1)/relativeErrorInUpN[b]));
3246 // relativeErrorInUp->SetBinError(b+1, 0.);
3247  relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorOutUp->GetBinContent(b+1)/relativeErrorOutUpN[b]));
3248 // relativeErrorOutUp->SetBinError(b+1, .0);
3249  relativeErrorInLow->SetBinContent(b+1, -1.*TMath::Sqrt(relativeErrorInLow->GetBinContent(b+1)/relativeErrorInLowN[b]));
3250 // relativeErrorInLow->SetBinError(b+1, 0.);
3251  relativeErrorOutLow->SetBinContent(b+1, -1.*TMath::Sqrt(relativeErrorOutLow->GetBinContent(b+1)/relativeErrorOutLowN[b]));
3252 // relativeErrorOutLow->SetBinError(b+1, .0);
3253  } else if (fSymmRMS) {
3254  if(relativeErrorInUpN[b] < 1) relativeErrorInUpN[b] = 1;
3255  if(relativeErrorOutUpN[b] < 1) relativeErrorOutUpN[b] = 1;
3256  relativeErrorInUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(b+1)/relativeErrorInUpN[b]));
3257  relativeErrorOutUp->SetBinContent(b+1, TMath::Sqrt(relativeErrorOutUp->GetBinContent(b+1)/relativeErrorOutUpN[b]));
3258  }
3259  }
3260  }
3261  relativeErrorInUp->SetYTitle("relative uncertainty");
3262  relativeErrorOutUp->SetYTitle("relative uncertainty");
3263  relativeErrorInLow->SetYTitle("relative uncertainty");
3264  relativeErrorOutLow->SetYTitle("relative uncertainty");
3265  relativeErrorInUp->GetYaxis()->SetRangeUser(-1.5, 3.);
3266  relativeErrorInLow->GetYaxis()->SetRangeUser(-1.5, 3.);
3267  relativeErrorOutUp->GetYaxis()->SetRangeUser(-1.5, 3.);
3268  relativeErrorOutLow->GetYaxis()->SetRangeUser(-1.5, 3.);
3269 
3270  canvasNominalMasterIn->cd(2);
3271  Style(gPad, "GRID");
3272  Style(relativeErrorInUp, kBlue, kBar);
3273  Style(relativeErrorInLow, kGreen, kBar);
3274  relativeErrorInUp->DrawCopy("b");
3275  relativeErrorInLow->DrawCopy("same b");
3276  Style(AddLegend(gPad));
3277  SavePadToPDF(canvasNominalMasterIn);
3278  canvasNominalMasterIn->Write();
3279  canvasNominalMasterOut->cd(2);
3280  Style(gPad, "GRID");
3281  Style(relativeErrorOutUp, kBlue, kBar);
3282  Style(relativeErrorOutLow, kGreen, kBar);
3283  relativeErrorOutUp->DrawCopy("b");
3284  relativeErrorOutLow->DrawCopy("same b");
3285  Style(AddLegend(gPad));
3286  SavePadToPDF(canvasNominalMasterOut);
3287  canvasNominalMasterOut->Write();
3288 }
3289 //_____________________________________________________________________________
3291  TArrayI* variationsIn, // variantions in plane
3292  TArrayI* variationsOut, // variantions out of plane
3293  TH1D*& relativeErrorInUp, // will store systematic error on v2
3294  TH1D*& relativeErrorInLow, // remains NULL
3295  TH1D*& relativeErrorOutUp, // remains NULL
3296  TH1D*& relativeErrorOutLow, // remains NULL
3297  TH1D*& relativeStatisticalErrorIn, // remains NULL
3298  TH1D*& relativeStatisticalErrorOut, // remains NULL
3299  TH1D*& nominal, // nominal v2
3300  TH1D*& nominalIn, // remains NULL
3301  TH1D*& nominalOut, // remains NULL
3302  Int_t columns, // trivial
3303  Float_t rangeLow, // trivial
3304  Float_t rangeUp, // trivial
3305  TFile* readMe, // trivial
3306  TString source, // trivial
3307  Bool_t RMS // NOT trivial
3308  ) const
3309 {
3310 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
3311  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
3312 #endif
3313  // 0) intermediate systematic check function. first make sure a placeholder is ready
3314  if(relativeErrorInUp) delete relativeErrorInUp;
3315  relativeErrorInUp = new TH1D(Form("absolute_systematic_uncertainty_%s", source.Data()), Form("absolute_systematic_uncertainty_%s", source.Data()), fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
3316 
3317  // 1) then loop over the variations
3318  for(Int_t i(0); i < variationsIn->GetSize(); i++) {
3319  printf(" variation %i of %i \n", i, variationsIn->GetSize());
3320  // the frist variation will be stored as nominal, others are looped over
3321  Int_t iIn[] = {variationsIn->At(i), variationsIn->At(i)};
3322  Int_t iOut[] = {variationsOut->At(i), variationsOut->At(i)};
3323 
3324  // call the functions and set the relevant pointer references
3325  TH1D* dud(0x0);
3326  DoIntermediateSystematics( // i guess this is the point where i retreive the variation
3327  new TArrayI(2, iIn), // except for i is one, that's the nominal measurement
3328  new TArrayI(2, iOut),
3329  dud, dud, dud, dud, dud, dud, dud,
3330  nominalIn,
3331  nominalOut,
3332  1,
3333  fBinsTrue->At(0),
3334  fBinsTrue->At(fBinsTrue->GetSize()-1),
3335  readMe,
3336  Form("error_on_v2_variation_%i", i));
3337 
3338  printf(" nominal in %i binsX \n", nominalIn->GetNbinsX());
3339  printf(" nominal out %i binsX \n", nominalOut->GetNbinsX());
3340 
3341  // 2) get a histo with the v2
3342  TH1D* v2(GetV2Histo(nominalIn, nominalOut, fEventPlaneRes, "nominal v_{2}"));
3343 
3344  printf(" v2 in %i binsX \n", v2->GetNbinsX());
3345  // 3) bounce if nominal measurement ...
3346  if(i==0) {
3347  nominal = v2;
3348  } else { // 4) calculate the uncertainty. for now only RMS style, easier and the first check
3349  for(Int_t k(0); k < nominal->GetNbinsX(); k++) {
3350  // 5) store the sum of squares of difference (so orig + square of new point) differentially
3351  relativeErrorInUp->SetBinContent(k+1, relativeErrorInUp->GetBinContent(k+1) + TMath::Power(v2->GetBinContent(k+1)-nominal->GetBinContent(k+1), 2));
3352  printf(" nominal bin %i is %.4f \n", k+1, nominal->GetBinContent(k+1));
3353  printf(" variati bin %i is %.4f \n", k+1, v2->GetBinContent(k+1));
3354  printf(" sum of squared difference is %.4f \n", relativeErrorInUp->GetBinContent(k+1));
3355  }
3356  delete v2;
3357  }
3358  }
3359  // 6) looped over all the variations, now divide by N and take the square for the RMS
3360  // to make life fun, the function that will add all uncertainties uses bin ERRORS, not content,
3361  // so we'll set the error here and set the relative error as bin content
3362  for(Int_t k(0); k < nominal->GetNbinsX(); k++) {
3363  relativeErrorInUp->SetBinError(k+1, TMath::Sqrt(relativeErrorInUp->GetBinContent(k+1)/((double)variationsIn->GetSize())));
3364  printf(" total absolute error in bin %i is %.4f \n", k+1, relativeErrorInUp->GetBinError(k+1));
3365  relativeErrorInUp->SetBinContent(k+1, relativeErrorInUp->GetBinError(k+1)/nominal->GetBinContent(k+1));
3366  }
3367  if(fConstantUE) {
3368  Int_t ctr(0);
3369  Double_t av(0), avct(0);
3370  for(Int_t k(3); k < 8; k++) {
3371  av+=relativeErrorInUp->GetBinError(k+1);
3372  avct+=relativeErrorInUp->GetBinContent(k+1);
3373  ctr++;
3374  }
3375  av/=((double)ctr);
3376  avct/=((double)ctr);
3377  for(Int_t k(8); k < nominal->GetNbinsX(); k++) {
3378  relativeErrorInUp->SetBinError(k+1, av);
3379  relativeErrorInUp->SetBinContent(k+1, avct);
3380  }
3381  }
3382 }
3383 //_____________________________________________________________________________
3384 void AliJetFlowTools::PostProcess(TString def, Int_t columns, Float_t rangeLow, Float_t rangeUp, TString in, TString out) const
3385 {
3386 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
3387  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
3388 #endif
3389  // go through the output file and perform post processing routines
3390  // can either be performed in one go with the unfolding, or at a later stage
3391  if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
3392  TFile readMe(in.Data(), "READ"); // open file read-only
3393  if(readMe.IsZombie()) {
3394  printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
3395  return;
3396  }
3397  printf("\n\n\n\t\t POSTPROCESSING \n > Recovered the following file structure : \n <");
3398  readMe.ls();
3399  TList* listOfKeys((TList*)readMe.GetListOfKeys());
3400  if(!listOfKeys) {
3401  printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
3402  return;
3403  }
3404  // prepare necessary canvasses
3405  TCanvas* canvasIn(new TCanvas("PearsonIn", "PearsonIn"));
3406  TCanvas* canvasOut(0x0);
3407  if(fDphiUnfolding) canvasOut = new TCanvas("PearsonOut", "PearsonOut");
3408  TCanvas* canvasRatioMeasuredRefoldedIn(new TCanvas("RefoldedIn", "RefoldedIn"));
3409  TCanvas* canvasRatioMeasuredRefoldedOut(0x0);
3410  if(fDphiUnfolding) canvasRatioMeasuredRefoldedOut = new TCanvas("RefoldedOut", "RefoldedOut");
3411  TCanvas* canvasSpectraIn(new TCanvas("SpectraIn", "SpectraIn"));
3412  TCanvas* canvasSpectraOut(0x0);
3413  if(fDphiUnfolding) canvasSpectraOut = new TCanvas("SpectraOut", "SpectraOut");
3414  TCanvas* canvasRatio(0x0);
3415  if(fDphiUnfolding) canvasRatio = new TCanvas("Ratio", "Ratio");
3416  TCanvas* canvasV2(0x0);
3417  if(fDphiUnfolding) canvasV2 = new TCanvas("V2", "V2");
3418  TCanvas* canvasMISC(new TCanvas("MISC", "MISC"));
3419  TCanvas* canvasMasterIn(new TCanvas("defaultIn", "defaultIn"));
3420  TCanvas* canvasMasterOut(0x0);
3421  if(fDphiUnfolding) canvasMasterOut = new TCanvas("defaultOut", "defaultOut");
3422  (fDphiUnfolding) ? canvasMISC->Divide(4, 2) : canvasMISC->Divide(4, 1);
3423  TDirectoryFile* defDir(0x0);
3424  TCanvas* canvasProfiles(new TCanvas("canvasProfiles", "canvasProfiles"));
3425  canvasProfiles->Divide(2);
3426  TProfile* ratioProfile(new TProfile("ratioProfile", "ratioProfile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3427  TProfile* v2Profile(new TProfile("v2Profile", "v2Profile", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3428  // get an estimate of the number of outputs and find the default set
3429  Int_t cacheMe(0);
3430  for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
3431  if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
3432  if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3433  cacheMe++;
3434  }
3435  }
3436  Int_t rows(TMath::Floor(cacheMe/(float)columns)+((cacheMe%columns)>0));
3437  canvasIn->Divide(columns, rows);
3438  if(canvasOut) canvasOut->Divide(columns, rows);
3439  canvasRatioMeasuredRefoldedIn->Divide(columns, rows);
3440  if(canvasRatioMeasuredRefoldedOut) canvasRatioMeasuredRefoldedOut->Divide(columns, rows);
3441  canvasSpectraIn->Divide(columns, rows);
3442  if(canvasSpectraOut) canvasSpectraOut->Divide(columns, rows);
3443  if(canvasRatio) canvasRatio->Divide(columns, rows);
3444  if(canvasV2) canvasV2->Divide(columns, rows);
3445 
3446  canvasMasterIn->Divide(columns, rows);
3447  if(canvasMasterOut) canvasMasterOut->Divide(columns, rows);
3448  // extract the default output
3449  TH1D* defaultUnfoldedJetSpectrumIn(0x0);
3450  TH1D* defaultUnfoldedJetSpectrumOut(0x0);
3451  if(defDir) {
3452  TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
3453  TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
3454  if(defDirIn) defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
3455  if(defDirOut) defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
3456  printf(" > succesfully extracted default results < \n");
3457  }
3458 
3459  // loop through the directories, only plot the graphs if the deconvolution converged
3460  TDirectoryFile* tempDir(0x0);
3461  TDirectoryFile* tempIn(0x0);
3462  TDirectoryFile* tempOut(0x0);
3463  for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
3464  // read the directory by its name
3465  tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3466  if(!tempDir) continue;
3467  TString dirName(tempDir->GetName());
3468  // try to read the in- and out of plane subdirs
3469  tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
3470  tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
3471  j++;
3472  if(!tempIn) { // attempt to get MakeAU output
3473  TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
3474  TCanvas* tempPearson(new TCanvas(Form("pearson_%s", dirName.Data()), Form("pearson_%s", dirName.Data())));
3475  tempPearson->Divide(4,2);
3476  TCanvas* tempRatio(new TCanvas(Form("ratio_%s", dirName.Data()), Form("ratio_%s", dirName.Data())));
3477  tempRatio->Divide(4,2);
3478  TCanvas* tempResult(new TCanvas(Form("result_%s", dirName.Data()), Form("result_%s", dirName.Data())));
3479  tempResult->Divide(4,2);
3480  for(Int_t q(0); q < 8; q++) {
3481  tempIn = (TDirectoryFile*)tempDir->Get(Form("%s___%s", stringArray[q].Data(), dirName.Data()));
3482  if(tempIn) {
3483  // to see if the unfolding converged try to extract the pearson coefficients
3484  TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
3485  if(pIn) {
3486  printf(" - %s in plane converged \n", dirName.Data());
3487  tempPearson->cd(1+q);
3488  Style(gPad, "PEARSON");
3489  pIn->DrawCopy("colz");
3490  TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3491  if(rIn) {
3492  Style(rIn);
3493  printf(" > found RatioRefoldedMeasured < \n");
3494  tempRatio->cd(q+1);
3495  rIn->SetFillColor(kRed);
3496  rIn->Draw("ap");
3497  }
3498  TH1D* dvector((TH1D*)tempIn->Get("dVector"));
3499  TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
3500  TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
3501  TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
3502  if(dvector && avalue && rm && eff) {
3503  Style(dvector);
3504  Style(avalue);
3505  Style(rm);
3506  Style(eff);
3507  canvasMISC->cd(1);
3508  Style(gPad, "SPECTRUM");
3509  dvector->DrawCopy();
3510  canvasMISC->cd(2);
3511  Style(gPad, "SPECTRUM");
3512  avalue->DrawCopy();
3513  canvasMISC->cd(3);
3514  Style(gPad, "PEARSON");
3515  rm->DrawCopy("colz");
3516  canvasMISC->cd(4);
3517  Style(gPad, "GRID");
3518  eff->DrawCopy();
3519  } else if(rm && eff) {
3520  Style(rm);
3521  Style(eff);
3522  canvasMISC->cd(3);
3523  Style(gPad, "PEARSON");
3524  rm->DrawCopy("colz");
3525  canvasMISC->cd(4);
3526  Style(gPad, "GRID");
3527  eff->DrawCopy();
3528  }
3529  }
3530  TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
3531  TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
3532  TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
3533  if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3534  if(defaultUnfoldedJetSpectrumIn) {
3535  Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
3536  TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
3537  temp->Divide(unfoldedSpectrum);
3538  temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3539  temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3540  temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3541  canvasMasterIn->cd(j);
3542  temp->GetYaxis()->SetRangeUser(0., 2);
3543  temp->DrawCopy();
3544  }
3545  TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
3546  tempResult->cd(q+1);
3547  Style(gPad);
3548  Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3549  unfoldedSpectrum->DrawCopy();
3550  Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3551  inputSpectrum->DrawCopy("same");
3552  Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3553  refoldedSpectrum->DrawCopy("same");
3554  TLegend* l(AddLegend(gPad));
3555  Style(l);
3556  if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3557  Float_t chi(fitStatus->GetBinContent(1));
3558  Float_t pen(fitStatus->GetBinContent(2));
3559  Int_t dof((int)fitStatus->GetBinContent(3));
3560  Float_t beta(fitStatus->GetBinContent(4));
3561  l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3562  } else if (fitStatus) { // only available in SVD fit
3563  Int_t reg((int)fitStatus->GetBinContent(1));
3564  l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3565  }
3566  }
3567  }
3568  }
3569  }
3570  if(tempIn) {
3571  // to see if the unfolding converged try to extract the pearson coefficients
3572  TH2D* pIn((TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data())));
3573  if(pIn) {
3574  printf(" - %s in plane converged \n", dirName.Data());
3575  canvasIn->cd(j);
3576  Style(gPad, "PEARSON");
3577  pIn->DrawCopy("colz");
3578  TGraphErrors* rIn((TGraphErrors*)tempIn->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3579  if(rIn) {
3580  Style(rIn);
3581  printf(" > found RatioRefoldedMeasured < \n");
3582  canvasRatioMeasuredRefoldedIn->cd(j);
3583  rIn->SetFillColor(kRed);
3584  rIn->Draw("ap");
3585  }
3586  TH1D* dvector((TH1D*)tempIn->Get("dVector"));
3587  TH1D* avalue((TH1D*)tempIn->Get("SingularValuesOfAC"));
3588  TH2D* rm((TH2D*)tempIn->Get(Form("ResponseMatrixIn_%s", dirName.Data())));
3589  TH1D* eff((TH1D*)tempIn->Get(Form("KinematicEfficiencyIn_%s", dirName.Data())));
3590  if(dvector && avalue && rm && eff) {
3591  Style(dvector);
3592  Style(avalue);
3593  Style(rm);
3594  Style(eff);
3595  canvasMISC->cd(1);
3596  Style(gPad, "SPECTRUM");
3597  dvector->DrawCopy();
3598  canvasMISC->cd(2);
3599  Style(gPad, "SPECTRUM");
3600  avalue->DrawCopy();
3601  canvasMISC->cd(3);
3602  Style(gPad, "PEARSON");
3603  rm->DrawCopy("colz");
3604  canvasMISC->cd(4);
3605  Style(gPad, "GRID");
3606  eff->DrawCopy();
3607  } else if(rm && eff) {
3608  Style(rm);
3609  Style(eff);
3610  canvasMISC->cd(3);
3611  Style(gPad, "PEARSON");
3612  rm->DrawCopy("colz");
3613  canvasMISC->cd(4);
3614  Style(gPad, "GRID");
3615  eff->DrawCopy();
3616  }
3617  }
3618  TH1D* inputSpectrum((TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data())));
3619  TH1D* unfoldedSpectrum((TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data())));
3620  TH1D* refoldedSpectrum((TH1D*)tempIn->Get(Form("RefoldedSpectrum_in_%s", dirName.Data())));
3621  if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3622  if(defaultUnfoldedJetSpectrumIn) {
3623  Style(defaultUnfoldedJetSpectrumIn, kBlue, kUnfoldedSpectrum);
3624  TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumIn->Clone(Form("defaultUnfoldedJetSpectrumIn_%s", dirName.Data())));
3625  temp->Divide(unfoldedSpectrum);
3626  temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3627  temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3628  temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3629  canvasMasterIn->cd(j);
3630  temp->GetYaxis()->SetRangeUser(0., 2);
3631  temp->DrawCopy();
3632  }
3633  TH1F* fitStatus((TH1F*)tempIn->Get(Form("fitStatus_%s_in", dirName.Data())));
3634  canvasSpectraIn->cd(j);
3635  Style(gPad);
3636  Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3637  unfoldedSpectrum->DrawCopy();
3638  Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3639  inputSpectrum->DrawCopy("same");
3640  Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3641  refoldedSpectrum->DrawCopy("same");
3642  TLegend* l(AddLegend(gPad));
3643  Style(l);
3644  if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3645  Float_t chi(fitStatus->GetBinContent(1));
3646  Float_t pen(fitStatus->GetBinContent(2));
3647  Int_t dof((int)fitStatus->GetBinContent(3));
3648  Float_t beta(fitStatus->GetBinContent(4));
3649  l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3650  } else if (fitStatus) { // only available in SVD fit
3651  Int_t reg((int)fitStatus->GetBinContent(1));
3652  l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3653  }
3654  }
3655  }
3656  if(tempOut) {
3657  TH2D* pOut((TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data())));
3658  if(pOut) {
3659  printf(" - %s out of plane converged \n", dirName.Data());
3660  canvasOut->cd(j);
3661  Style(gPad, "PEARSON");
3662  pOut->DrawCopy("colz");
3663  TGraphErrors* rOut((TGraphErrors*)tempOut->Get(Form("RatioRefoldedMeasured_%s", dirName.Data())));
3664  if(rOut) {
3665  Style(rOut);
3666  printf(" > found RatioRefoldedMeasured < \n");
3667  canvasRatioMeasuredRefoldedOut->cd(j);
3668  rOut->SetFillColor(kRed);
3669  rOut->Draw("ap");
3670  }
3671  TH1D* dvector((TH1D*)tempOut->Get("dVector"));
3672  TH1D* avalue((TH1D*)tempOut->Get("SingularValuesOfAC"));
3673  TH2D* rm((TH2D*)tempOut->Get(Form("ResponseMatrixOut_%s", dirName.Data())));
3674  TH1D* eff((TH1D*)tempOut->Get(Form("KinematicEfficiencyOut_%s", dirName.Data())));
3675  if(dvector && avalue && rm && eff) {
3676  Style(dvector);
3677  Style(avalue);
3678  Style(rm);
3679  Style(eff);
3680  canvasMISC->cd(5);
3681  Style(gPad, "SPECTRUM");
3682  dvector->DrawCopy();
3683  canvasMISC->cd(6);
3684  Style(gPad, "SPECTRUM");
3685  avalue->DrawCopy();
3686  canvasMISC->cd(7);
3687  Style(gPad, "PEARSON");
3688  rm->DrawCopy("colz");
3689  canvasMISC->cd(8);
3690  Style(gPad, "GRID");
3691  eff->DrawCopy();
3692  } else if(rm && eff) {
3693  Style(rm);
3694  Style(eff);
3695  canvasMISC->cd(7);
3696  Style(gPad, "PEARSON");
3697  rm->DrawCopy("colz");
3698  canvasMISC->cd(8);
3699  Style(gPad, "GRID");
3700  eff->DrawCopy();
3701  }
3702  }
3703  TH1D* inputSpectrum((TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data())));
3704  TH1D* unfoldedSpectrum((TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data())));
3705  TH1D* refoldedSpectrum((TH1D*)tempOut->Get(Form("RefoldedSpectrum_out_%s", dirName.Data())));
3706  if(inputSpectrum && unfoldedSpectrum && refoldedSpectrum) {
3707  if(defaultUnfoldedJetSpectrumOut) {
3708  Style(defaultUnfoldedJetSpectrumOut, kBlue, kUnfoldedSpectrum);
3709  TH1D* temp((TH1D*)defaultUnfoldedJetSpectrumOut->Clone(Form("defaultUnfoldedJetSpectrumOut_%s", dirName.Data())));
3710  temp->Divide(unfoldedSpectrum);
3711  temp->SetTitle(Form("ratio default unfolded / %s", dirName.Data()));
3712  temp->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
3713  temp->GetYaxis()->SetTitle(Form("%s / %s", def.Data(), dirName.Data()));
3714  canvasMasterOut->cd(j);
3715  temp->GetYaxis()->SetRangeUser(0., 2.);
3716  temp->DrawCopy();
3717  }
3718  TH1F* fitStatus((TH1F*)tempOut->Get(Form("fitStatus_%s_out", dirName.Data())));
3719  canvasSpectraOut->cd(j);
3720  Style(gPad);
3721  Style(unfoldedSpectrum, kRed, kUnfoldedSpectrum);
3722  unfoldedSpectrum->DrawCopy();
3723  Style(inputSpectrum, kGreen, kMeasuredSpectrum);
3724  inputSpectrum->DrawCopy("same");
3725  Style(refoldedSpectrum, kBlue, kFoldedSpectrum);
3726  refoldedSpectrum->DrawCopy("same");
3727  TLegend* l(AddLegend(gPad));
3728  Style(l);
3729  if(fitStatus && fitStatus->GetNbinsX() == 4) { // only available in chi2 fit
3730  Float_t chi(fitStatus->GetBinContent(1));
3731  Float_t pen(fitStatus->GetBinContent(2));
3732  Int_t dof((int)fitStatus->GetBinContent(3));
3733  Float_t beta(fitStatus->GetBinContent(4));
3734  l->AddEntry((TObject*)0, Form("#chi %.2f \tP %.2f \tDOF %i, #beta %.2f", chi, pen, dof, beta), "");
3735  } else if (fitStatus) { // only available in SVD fit
3736  Int_t reg((int)fitStatus->GetBinContent(1));
3737  l->AddEntry((TObject*)0, Form("REG %i", reg), "");
3738  }
3739  }
3740  }
3741  if(canvasRatio && canvasV2) {
3742  canvasRatio->cd(j);
3743  TGraphErrors* ratioYield((TGraphErrors*)tempDir->Get(Form("RatioInOutPlane_%s", dirName.Data())));
3744  Double_t _tempx(0), _tempy(0);
3745  if(ratioYield) {
3746  Style(ratioYield);
3747  for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3748  ratioYield->GetPoint(b, _tempx, _tempy);
3749  ratioProfile->Fill(_tempx, _tempy);
3750  }
3751  ratioProfile->GetYaxis()->SetRangeUser(-0., 2.);
3752  ratioProfile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3753  ratioYield->GetYaxis()->SetRangeUser(-0., 2.);
3754  ratioYield->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3755  ratioYield->SetFillColor(kRed);
3756  ratioYield->Draw("ap");
3757  }
3758  canvasV2->cd(j);
3759  TGraphErrors* ratioV2((TGraphErrors*)tempDir->Get(Form("v2_%s", dirName.Data())));
3760  if(ratioV2) {
3761  Style(ratioV2);
3762  for(Int_t b(0); b < fBinsTrue->GetSize(); b++) {
3763  ratioV2->GetPoint(b, _tempx, _tempy);
3764  v2Profile->Fill(_tempx, _tempy);
3765 
3766  }
3767  v2Profile->GetYaxis()->SetRangeUser(-0., 2.);
3768  v2Profile->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3769  ratioV2->GetYaxis()->SetRangeUser(-.25, .75);
3770  ratioV2->GetXaxis()->SetRangeUser(rangeLow, rangeUp);
3771  ratioV2->SetFillColor(kRed);
3772  ratioV2->Draw("ap");
3773  }
3774  }
3775  }
3776  TFile output(out.Data(), "RECREATE");
3777  canvasProfiles->cd(1);
3778  Style(ratioProfile);
3779  ratioProfile->DrawCopy();
3780  canvasProfiles->cd(2);
3781  Style(v2Profile);
3782  v2Profile->DrawCopy();
3783  SavePadToPDF(canvasProfiles);
3784  canvasProfiles->Write();
3785  SavePadToPDF(canvasIn);
3786  canvasIn->Write();
3787  if(canvasOut) {
3788  SavePadToPDF(canvasOut);
3789  canvasOut->Write();
3790  }
3791  SavePadToPDF(canvasRatioMeasuredRefoldedIn);
3792  canvasRatioMeasuredRefoldedIn->Write();
3793  if(canvasRatioMeasuredRefoldedOut) {
3794  SavePadToPDF(canvasRatioMeasuredRefoldedOut);
3795  canvasRatioMeasuredRefoldedOut->Write();
3796  }
3797  SavePadToPDF(canvasSpectraIn);
3798  canvasSpectraIn->Write();
3799  if(canvasSpectraOut) {
3800  SavePadToPDF(canvasSpectraOut);
3801  canvasSpectraOut->Write();
3802  SavePadToPDF(canvasRatio);
3803  canvasRatio->Write();
3804  SavePadToPDF(canvasV2);
3805  canvasV2->Write();
3806  }
3807  SavePadToPDF(canvasMasterIn);
3808  canvasMasterIn->Write();
3809  if(canvasMasterOut) {
3810  SavePadToPDF(canvasMasterOut);
3811  canvasMasterOut->Write();
3812  }
3813  SavePadToPDF(canvasMISC);
3814  canvasMISC->Write();
3815  output.Write();
3816  output.Close();
3817 }
3818 //_____________________________________________________________________________
3820 {
3821 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
3822  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
3823 #endif
3824  // function to interpret results of bootstrapping routine
3825  // TString def should hold the true emperical distribution
3826  if(fOutputFile && !fOutputFile->IsZombie()) fOutputFile->Close();
3827  TFile readMe(in.Data(), "READ"); // open file read-only
3828  if(readMe.IsZombie()) {
3829  printf(" > Fatal error, couldn't read %s for post processing ! < \n", in.Data());
3830  return;
3831  }
3832  printf("\n\n\n\t\t BootstrapSpectra \n > Recovered the following file structure : \n <");
3833  readMe.ls();
3834  TList* listOfKeys((TList*)readMe.GetListOfKeys());
3835  if(!listOfKeys) {
3836  printf(" > Fatal error, couldn't retrieve list of keys. Input file might have been corrupted ! < \n");
3837  return;
3838  }
3839  // get an estimate of the number of outputs and find the default set
3840  TDirectoryFile* defDir(0x0);
3841  for(Int_t i(0); i < listOfKeys->GetSize(); i++) {
3842  if(dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()))) {
3843  if(!strcmp(listOfKeys->At(i)->GetName(), def.Data())) defDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3844  }
3845  }
3846 
3847  // extract the default output this is the 'emperical' distribution
3848  // w.r.t. which the other distributions will be evaluated
3849  TH1D* defaultUnfoldedJetSpectrumIn(0x0);
3850  TH1D* defaultUnfoldedJetSpectrumOut(0x0);
3851  TH1D* defaultInputSpectrumIn(0x0);
3852  TH1D* defaultInputSpectrumOut(0x0);
3853  TGraphErrors* v2Emperical(0x0);
3854  if(defDir) {
3855  TDirectoryFile* defDirIn = (TDirectoryFile*)defDir->Get(Form("InPlane___%s", def.Data()));
3856  TDirectoryFile* defDirOut = (TDirectoryFile*)defDir->Get(Form("OutOfPlane___%s", def.Data()));
3857  if(defDirIn) {
3858  defaultUnfoldedJetSpectrumIn = (TH1D*)defDirIn->Get(Form("UnfoldedSpectrum_in_%s", def.Data()));
3859  defaultInputSpectrumIn = (TH1D*)defDirIn->Get(Form("InputSpectrum_in_%s", def.Data()));
3860  }
3861  if(defDirOut) {
3862  defaultUnfoldedJetSpectrumOut = (TH1D*)defDirOut->Get(Form("UnfoldedSpectrum_out_%s", def.Data()));
3863  defaultInputSpectrumOut = (TH1D*)defDirOut->Get(Form("InputSpectrum_out_%s", def.Data()));
3864  }
3865  }
3866  if((!defaultUnfoldedJetSpectrumIn && defaultUnfoldedJetSpectrumOut)) {
3867  printf(" BootstrapSpectra: couldn't extract default spectra, aborting! \n");
3868  return;
3869  }
3870  else v2Emperical = GetV2(defaultUnfoldedJetSpectrumIn, defaultUnfoldedJetSpectrumOut, fEventPlaneRes);
3871 
3872  // now that we know for sure that the input is in place, reserve the bookkeeping histograms
3873  TH1F* delta0[fBinsTrue->GetSize()-1];
3874  for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3875  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);
3876  delta0[i]->Sumw2();
3877  }
3878  // and a canvas for illustration purposes only
3879  TCanvas* spectraIn(new TCanvas("spectraIn", "spectraIn"));
3880  TCanvas* spectraOut(new TCanvas("spectraOut", "spectraOut"));
3881  // common reference (in this case the generated v2)
3882  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);
3883 
3884  // loop through the directories, only plot the graphs if the deconvolution converged
3885  TDirectoryFile* tempDir(0x0);
3886  TDirectoryFile* tempIn(0x0);
3887  TDirectoryFile* tempOut(0x0);
3888  TH1D* unfoldedSpectrumIn(0x0);
3889  TH1D* unfoldedSpectrumOut(0x0);
3890  TH1D* measuredSpectrumIn(0x0);
3891  TH1D* measuredSpectrumOut(0x0);
3892  for(Int_t i(0), j(0); i < listOfKeys->GetSize(); i++) {
3893  // read the directory by its name
3894  tempDir = dynamic_cast<TDirectoryFile*>(readMe.Get(listOfKeys->At(i)->GetName()));
3895  if(!tempDir) continue;
3896  TString dirName(tempDir->GetName());
3897  // read only bootstrapped outputs
3898  if(!dirName.Contains("bootstrap")) continue;
3899  // try to read the in- and out of plane subdirs
3900  tempIn = (TDirectoryFile*)tempDir->Get(Form("InPlane___%s", dirName.Data()));
3901  tempOut = (TDirectoryFile*)tempDir->Get(Form("OutOfPlane___%s", dirName.Data()));
3902  j++;
3903  // extract the unfolded spectra only if both in- and out-of-plane converted (i.e. pearson coefficients were saved)
3904  if(tempIn) {
3905  if(!(TH2D*)tempIn->Get(Form("PearsonCoefficients_in_%s", dirName.Data()))) continue;
3906  unfoldedSpectrumIn = (TH1D*)tempIn->Get(Form("UnfoldedSpectrum_in_%s", dirName.Data()));
3907  measuredSpectrumIn = (TH1D*)tempIn->Get(Form("InputSpectrum_in_%s", dirName.Data()));
3908  spectraIn->cd();
3909  (j==1) ? measuredSpectrumIn->DrawCopy() : measuredSpectrumIn->DrawCopy("same");
3910  }
3911  if(tempOut) {
3912  if(!(TH2D*)tempOut->Get(Form("PearsonCoefficients_out_%s", dirName.Data()))) continue;
3913  unfoldedSpectrumOut = (TH1D*)tempOut->Get(Form("UnfoldedSpectrum_out_%s", dirName.Data()));
3914  measuredSpectrumOut = (TH1D*)tempOut->Get(Form("InputSpectrum_out_%s", dirName.Data()));
3915  spectraOut->cd();
3916  (j==1) ? measuredSpectrumOut->DrawCopy() : measuredSpectrumOut->DrawCopy("same");
3917  }
3918  // get v2 with statistical uncertainties from the extracted spectra
3919  TGraphErrors* v2Bootstrapped(GetV2(unfoldedSpectrumIn, unfoldedSpectrumOut, fEventPlaneRes));
3920  // and loop over all bins to fill the bookkeeping histograms
3921  Double_t yBoot(0), yEmp(0), xDud(0);
3922  // optional: common reference (in this case the sampled v2 value)
3923  for(Int_t k(0); k < fBinsTrue->GetSize()-1; k++) {
3924  // read values point by point (passed by reference)
3925  v2Bootstrapped->GetPoint(k, xDud, yBoot);
3926  v2Emperical->GetPoint(k, xDud, yEmp);
3927  if(commonReference) {
3928  if(!commonReference->Eval(xDud)<1e-10) {
3929  yEmp/=commonReference->Eval(xDud);
3930  yBoot/=commonReference->Eval(xDud);
3931  } else { // if reference equals zero, take emperical distribution as reference
3932  yEmp/=yEmp; // 1
3933  yBoot/=yEmp;
3934  }
3935  }
3936  cout << " yEmp " << yEmp << " yBoot " << yBoot << endl;
3937  // save relative difference per pt bin
3938  if(TMath::Abs(yBoot)>1e-10) delta0[k]->Fill(yEmp - yBoot);
3939  }
3940  }
3941  // extracting final results now, as first estimate just a gaus fit to the distributions
3942  // (should be changed perhaps to proper rms eventually)
3943  // attach relevant data to current buffer in the same loop
3944  TFile output(out.Data(), "RECREATE");
3945  defaultInputSpectrumIn->SetLineColor(kRed);
3946  spectraIn->cd();
3947  defaultInputSpectrumIn->DrawCopy("same");
3948  defaultInputSpectrumOut->SetLineColor(kRed);
3949  spectraOut->cd();
3950  defaultInputSpectrumOut->DrawCopy("same");
3951  spectraIn->Write();
3952  spectraOut->Write();
3953 
3954  TH1F* delta0spread(new TH1F("delta0spread", "#sigma(#Delta_{0})", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3955  TH1F* unfoldingError(new TH1F("unfoldingError", "error from unfolding algorithm", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
3956  TF1* gaus(new TF1("gaus", "gaus"/*[0]*exp(-0.5*((x-[1])/[2])**2)"*/, -1., 1.));
3957  Double_t xDud(0), yDud(0);
3958  for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
3959  delta0[i]->Fit(gaus);
3960  delta0[i]->GetYaxis()->SetTitle("counts");
3961  delta0[i]->GetXaxis()->SetTitle("(v_{2, jet}^{measured} - v_{2, jet}^{generated}) / input v_{2}");
3962  delta0spread->SetBinContent(i+1, gaus->GetParameter(1)); // mean of gaus
3963  delta0spread->SetBinError(i+1, gaus->GetParameter(2)); // sigma of gaus
3964  cout << " mean " << gaus->GetParameter(1) << endl;
3965  cout << " sigm " << gaus->GetParameter(2) << endl;
3966  delta0[i]->Write();
3967  v2Emperical->GetPoint(i, xDud, yDud);
3968  unfoldingError->SetBinContent(i+1, 1e-10/*gaus->GetParameter(1)*/);
3969  if(commonReference && !commonReference->Eval(xDud)<1e-10) unfoldingError->SetBinError(i+1, v2Emperical->GetErrorY(i)/(commonReference->Eval(xDud)));
3970  else if(yDud>10e-10) unfoldingError->SetBinError(i+1, v2Emperical->GetErrorY(i)/yDud);
3971  else unfoldingError->SetBinError(i+1, 0.);
3972  }
3973  delta0spread->GetXaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
3974  delta0spread->GetYaxis()->SetTitle("(mean v_{2, jet}^{measured} - v_{2, jet}^{generated}) / input v_{2}");
3975  delta0spread->Write();
3976  unfoldingError->GetXaxis()->SetTitle("p_{T}^{jet} (GeV/c)");
3977  unfoldingError->GetYaxis()->SetTitle("(mean v_{2, jet}^{measured} - v_{2, jet}^{generated}) / input v_{2}");
3978  unfoldingError->Write();
3979  // write the buffer and close the file
3980  output.Write();
3981  output.Close();
3982 }
3983 //_____________________________________________________________________________
3985  TH2D* detectorResponse, // detector response matrix
3986  TH1D* jetPtIn, // in plane jet spectrum
3987  TH1D* jetPtOut, // out of plane jet spectrum
3988  TH1D* dptIn, // in plane delta pt distribution
3989  TH1D* dptOut, // out of plane delta pt distribution
3990  Int_t eventCount) {
3991 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
3992  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
3993 #endif
3994  // set input histograms manually
3995  fDetectorResponse = detectorResponse;
3996  fSpectrumIn = jetPtIn;
3997  fSpectrumOut = jetPtOut;
3998  fDptInDist = dptIn;
3999  fDptOutDist = dptOut;
4000  // check if all data is provided
4001  if(!fDetectorResponse) {
4002  printf(" fDetectorResponse not found \n ");
4003  return kFALSE;
4004  }
4005  // check if the pt bin for true and rec have been set
4006  if(!fBinsTrue || !fBinsRec) {
4007  printf(" No true or rec bins set, please set binning ! \n");
4008  return kFALSE;
4009  }
4010  if(!fRMSSpectrumIn) { // initialie the profiles which will hold the RMS values. if binning changes in between unfolding
4011  // procedures, these profiles will be nonsensical, user is responsible
4012  fRMSSpectrumIn = new TProfile("fRMSSpectrumIn", "fRMSSpectrumIn", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
4013  fRMSSpectrumOut = new TProfile("fRMSSpectrumOut", "fRMSSpectrumOut", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
4014  fRMSRatio = new TProfile("fRMSRatio", "fRMSRatio", fBinsTrue->GetSize()-1, fBinsTrue->GetArray());
4015  }
4016  // normalize spectra to event count if requested
4017  if(fNormalizeSpectra) {
4018  fEventCount = eventCount;
4019  if(fEventCount > 0) {
4020  fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
4021  fSpectrumOut->Sumw2();
4022  fSpectrumIn->Scale(1./((double)fEventCount));
4023  fSpectrumOut->Scale(1./((double)fEventCount));
4024  }
4025  }
4026  if(!fNormalizeSpectra && fEventCount > 0) {
4027  fSpectrumIn->Sumw2(); // necessary for correct error propagation of scale
4028  fSpectrumOut->Sumw2();
4029  fSpectrumIn->Scale(1./((double)fEventCount));
4030  fSpectrumOut->Scale(1./((double)fEventCount));
4031  }
4033  fDptIn->SetNameTitle(Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_INPLANE_%i", fCentralityArray->At(0)));
4034  fDptIn->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
4035  fDptIn->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
4037  fDptOut->SetNameTitle(Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)), Form("dpt_response_OUTOFPLANE_%i", fCentralityArray->At(0)));
4038  fDptOut->GetXaxis()->SetTitle("p_{T, jet}^{gen} [GeV/c]");
4039  fDptOut->GetYaxis()->SetTitle("p_{T, jet}^{rec} [GeV/c]");
4040 
4041 
4042  // set the flag to let other functions to read data from input file
4043  fRawInputProvided = kTRUE;
4044  return fRawInputProvided;
4045 }
4046 //_____________________________________________________________________________
4048 {
4049 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4050  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4051 #endif
4052  // return ratio of h1 / h2
4053  // histograms can have different binning. errors are propagated as uncorrelated
4054  if(!(h1 && h2) ) {
4055  printf(" GetRatio called with NULL argument(s) \n ");
4056  return 0x0;
4057  }
4058  Int_t j(0);
4059  TGraphErrors *gr = new TGraphErrors();
4060  gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
4061  Double_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
4062  TH1* dud((TH1*)h1->Clone("dud"));
4063  dud->Sumw2();
4064  h1->Sumw2();
4065  h2->Sumw2();
4066  if(!dud->Divide(h1, h2)) {
4067  printf(" > ROOT failed to divide two histogams, dividing manually \n < ");
4068  for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
4069  binCent = h1->GetXaxis()->GetBinCenter(i);
4070  if(xmax > 0. && binCent > xmax) continue;
4071  j = h2->FindBin(binCent);
4072  binWidth = h1->GetXaxis()->GetBinWidth(i);
4073  if(h2->GetBinContent(j) > 0.) {
4074  ratio = h1->GetBinContent(i)/h2->GetBinContent(j);
4075  Double_t A = h1->GetBinError(i)/h1->GetBinContent(i);
4076  Double_t B = h2->GetBinError(i)/h2->GetBinContent(i);
4077  error2 = ratio*ratio*A*A+ratio*ratio*B*B;
4078  if(error2 > 0 ) error2 = TMath::Sqrt(error2);
4079  gr->SetPoint(i-1, binCent, ratio);
4080  gr->SetPointError(i-1, 0.5*binWidth, error2);
4081  }
4082  }
4083  } else {
4084  printf( " > using ROOT to divide two histograms < \n");
4085  for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
4086  binCent = dud->GetXaxis()->GetBinCenter(i);
4087  if(xmax > 0. && binCent > xmax) continue;
4088  binWidth = dud->GetXaxis()->GetBinWidth(i);
4089  gr->SetPoint(i-1,binCent,dud->GetBinContent(i));
4090  gr->SetPointError(i-1, 0.5*binWidth,dud->GetBinError(i));
4091  }
4092  }
4093 
4094  if(appendFit) {
4095  TF1* fit(new TF1("lin", "pol0", 10, 100));
4096  gr->Fit(fit);
4097  }
4098  if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
4099  if(dud) delete dud;
4100  return gr;
4101 }
4102 //_____________________________________________________________________________
4104 {
4105 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4106  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4107 #endif
4108  // get v2 from difference of in plane, out of plane yield
4109  // h1 must hold the in-plane yield, h2 holds the out of plane yield
4110  // r is the event plane resolution for the chosen centrality
4111  if(!(h1 && h2) ) {
4112  printf(" GetV2 called with NULL argument(s) \n ");
4113  return 0x0;
4114  }
4115  TGraphErrors *gr = new TGraphErrors();
4116  gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
4117  Float_t binCent(0.), ratio(0.), error2(0.), binWidth(0.);
4118  Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
4119 
4120  for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
4121  binCent = h1->GetXaxis()->GetBinCenter(i);
4122  binWidth = h1->GetXaxis()->GetBinWidth(i);
4123  if(h2->GetBinContent(i) > 0.) {
4124  in = h1->GetBinContent(i);
4125  ein = h1->GetBinError(i);
4126  out = h2->GetBinContent(i);
4127  eout = h2->GetBinError(i);
4128  ratio = pre*((in-out)/(in+out));
4129  error2 = (4.*out*out/(TMath::Power(in+out, 4)))*ein*ein+(4.*in*in/(TMath::Power(in+out, 4)))*eout*eout;
4130  error2 = error2*pre*pre;
4131  if(error2 > 0) error2 = TMath::Sqrt(error2);
4132  gr->SetPoint(i-1,binCent,ratio);
4133  gr->SetPointError(i-1,0.5*binWidth,error2);
4134  }
4135  }
4136  if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
4137  return gr;
4138 }
4139 //_____________________________________________________________________________
4141 {
4142 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4143  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4144 #endif
4145  // get v2 from difference of in plane, out of plane yield
4146  // h1 must hold the in-plane yield, h2 holds the out of plane yield
4147  // r is the event plane resolution for the chosen centrality
4148  if(!(h1 && h2) ) {
4149  printf(" GetV2 called with NULL argument(s) \n ");
4150  return 0x0;
4151  }
4152  TH1D* gr((TH1D*)h1->Clone(name.Data()));
4153  gr->GetXaxis()->SetTitle("p_{T, jet} [GeV/c]");
4154  Float_t ratio(0.), error2(0.);
4155  Double_t pre(TMath::Pi()/(4.*r)), in(0.), out(0.), ein(0.), eout(0.);
4156 
4157  for(Int_t i(1); i <= h1->GetNbinsX(); i++) {
4158  if(h2->GetBinContent(i) > 0.) {
4159  in = h1->GetBinContent(i);
4160  ein = h1->GetBinError(i);
4161  out = h2->GetBinContent(i);
4162  eout = h2->GetBinError(i);
4163  ratio = pre*((in-out)/(in+out));
4164  error2 = (4.*out*out/(TMath::Power(in+out, 4)))*ein*ein+(4.*in*in/(TMath::Power(in+out, 4)))*eout*eout;
4165  error2 = error2*pre*pre;
4166  if(error2 > 0) error2 = TMath::Sqrt(error2);
4167  gr->SetBinContent(i,ratio);
4168  gr->SetBinError(i,error2);
4169  }
4170  }
4171  if(strcmp(name, "")) gr->SetNameTitle(name.Data(), name.Data());
4172  return gr;
4173 }
4174 //_____________________________________________________________________________
4176 {
4177 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4178  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4179 #endif
4180  // convert a tgrapherrors to a histogram
4181  if(!g) {
4182  printf(" > ConvertGraphToHistogram recevied a NULL pointer > \n");
4183  return 0x0;
4184  }
4185  // first get the frame which we'll use to build the histogram
4186  TH1F* hist(g->GetHistogram());
4187  Double_t xref(0), yref(0), yerr(0);
4188  // then copy the points and errors (remember graph starts at bin 0 and hist at 1)
4189  for(Int_t i(0); i < hist->GetNbinsX(); i++) {
4190  g->GetPoint(i, xref, yref);
4191  yerr = g->GetErrorY(i);
4192  hist->SetBinContent(i+1, yref);
4193  hist->SetBinError(i+1, yerr);
4194  }
4195  return hist;
4196 }
4197 //_____________________________________________________________________________
4199 {
4200 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4201  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4202 #endif
4203  // add the errors of a tgra histogram, assuming oncorrelated uncertainties
4204  // the points on the graph are returned !
4205  if(!(g&&h)) {
4206  printf(" > ConvertGraphToHistogram recevied a NULL pointer > \n");
4207  return 0x0;
4208  }
4209 
4210  printf("\n\n\n adding histo error to graph error \n\n\n\n");
4211 
4212 
4213  Double_t yerrL(0), yerrH(0), herr(0);
4214  // quadratic sum of the errors, assyming here symmetric ones
4215  for(Int_t i(0); i < h->GetNbinsX(); i++) {
4216  yerrH = g->GetErrorY(i);//high(i);
4217  yerrL = g->GetErrorY(i);//low(i);
4218  herr = h->GetBinError(i+1);
4219  cout << " graph error H " << yerrH << "\t histo error " << herr << endl;
4220  cout << " graph error L " << yerrL << "\t histo error " << herr << endl;
4221  yerrH = TMath::Sqrt(yerrH*yerrH+herr*herr);
4222  yerrL = TMath::Sqrt(yerrL*yerrL+herr*herr);
4223  g->SetPointError(i, g->GetErrorX(i), g->GetErrorX(i), yerrL, yerrH);
4224  }
4225  return g;
4226 }
4227 //_____________________________________________________________________________
4229 {
4230 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4231  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4232 #endif
4233  // get rms between a and b for x-axis of histo
4234  if(!h || b < a ) return 0.;
4235 
4236  // a is a lower edge and b an upper edge
4237  Int_t binA(h->GetXaxis()->FindBin(a+.1));
4238  Int_t binB(h->GetXaxis()->FindBin(b+.1)); // on purpose
4239  Double_t RMS(0), mean(0), c(0);
4240 
4241  for(Int_t i(binA); i < binB; i++) {
4242  c++;
4243  mean += h->GetBinContent(i);
4244  }
4245  if( c > 0. ) {
4246  mean /= c;
4247  for(Int_t i(binA); i < binB; i++) RMS += (mean-h->GetBinContent(i))*(mean-h->GetBinContent(i));
4248  return TMath::Sqrt(RMS/c);
4249  }
4250  return 0.;
4251 }
4252 //_____________________________________________________________________________
4254  Float_t pivot, Bool_t subdueError, TString str, Bool_t setContent)
4255 {
4256 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4257  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4258 #endif
4259  // return an error from a fit
4260  TF1* lin = new TF1("lin", Form("(x<%i)*(pol1)+(x>%i)*[2]", (int)pivot, (int)pivot), a, b);
4261  // clone the input
4262  TF1* fit_pol0 = new TF1("fit_pol0", "pol0", pivot, b);
4263  TF1* fit_pol1 = new TF1("fit_pol1", "pol1", a, pivot);
4264 
4265  TH1* h((TH1*)h1->Clone(str.Data()));
4266  // reset the guy and fill it with the maxima
4267  h->Reset();
4268  Double_t _a(0), _b(0);
4269  for(Int_t i(1); i < h->GetNbinsX() + 1; i++) {
4270  _a = TMath::Abs(h1->GetBinContent(i));
4271  _b = TMath::Abs(h2->GetBinContent(i));
4272  if(_a > _b) {
4273  h->SetBinContent(i, _a);
4274  h->SetBinError(i, h1->GetBinError(i));
4275  } else {
4276  h->SetBinContent(i, _b);
4277  h->SetBinError(i, h2->GetBinError(i));
4278  }
4279  }
4280 
4281 
4282 
4283  // fit to full error. root doesn't like fitting a step-function, so fit these
4284  // two components separately ...
4285  h->Fit(fit_pol0, "", "", pivot, b);
4286  lin->SetParameter(2, (subdueError) ? 0. : fit_pol0->GetParameter(0));
4287 
4288  h->Fit(fit_pol1, "", "", a, pivot);
4289 
4290  //check if it makes sense
4291  if(fit_pol1->GetParameter(1) > 0) {
4292  fit_pol1 = new TF1("fit_pol2", "pol0", a, pivot);
4293  h->Fit(fit_pol1, "", "", a, pivot);
4294  lin->SetParameter(0, fit_pol1->GetParameter(0));
4295  lin->SetParameter(1, 0);
4296  } else {
4297  lin->SetParameter(0, fit_pol1->GetParameter(0));
4298  lin->SetParameter(1, fit_pol1->GetParameter(1));
4299  }
4300 
4301 
4302  // .. and then add them together
4303  h->GetListOfFunctions()->Add(lin);
4304 
4305  if(setContent) {
4306  // update the histos with the fit result
4307  for(Int_t i(1); i < h->GetNbinsX() + 1; i++) {
4308  // calculate the integral in a given bin
4309  Double_t dud(lin->Integral(h->GetXaxis()->GetBinLowEdge(i),h->GetXaxis()->GetBinUpEdge(i))/h->GetXaxis()->GetBinWidth(i));
4310  // if it's larger than 0, assign an 'up' error
4311  h1->SetBinContent(i, 0);
4312  h2->SetBinContent(i, 0);
4313  h1->SetBinError(i, 0);
4314  h2->SetBinError(i, 0);
4315  // dont' show errors outside of range of interest
4316  if(a > h->GetXaxis()->GetBinLowEdge(i) || b < h->GetXaxis()->GetBinUpEdge(i)) continue;
4317  // assign them correctly to up or low
4318  if(dud > 0) {
4319  h1->SetBinContent(i, dud);
4320  } else {
4321  h2->SetBinContent(i, dud);
4322  }
4323  }
4324  h->Write();
4325  }
4326  return lin;
4327 }
4328 //_____________________________________________________________________________
4330  TH1* h1, TH1* h2, Double_t r, TString name,
4331  TH1* relativeErrorInUp,
4332  TH1* relativeErrorInLow,
4333  TH1* relativeErrorOutUp,
4334  TH1* relativeErrorOutLow,
4335  Float_t rho) const
4336 {
4337 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4338  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4339 #endif
4340  // get v2 with asymmetric systematic error
4341  // note that this is ONLY the systematic error, no statistical error!
4342  // rho is the pearson correlation coefficient
4343  TGraphErrors* tempV2(GetV2(h1, h2, r, name));
4344  Double_t* ax = new Double_t[fBinsTrue->GetSize()-1];
4345  Double_t* ay = new Double_t[fBinsTrue->GetSize()-1];
4346  Double_t* axh = new Double_t[fBinsTrue->GetSize()-1];
4347  Double_t* axl = new Double_t[fBinsTrue->GetSize()-1];
4348  Double_t* ayh = new Double_t[fBinsTrue->GetSize()-1];
4349  Double_t* ayl = new Double_t[fBinsTrue->GetSize()-1];
4350  Double_t in(0.), out(0.), einUp(0.), einLow(0.), eoutUp(0.), eoutLow(0.), error2Up(0.), error2Low(0.);
4351  // loop through the bins and do error propagation
4352  for(Int_t i(0); i < fBinsTrue->GetSize()-1; i++) {
4353  // extract the absolute errors
4354  in = h1->GetBinContent(i+1);
4355  einUp = TMath::Abs(in*relativeErrorInUp->GetBinContent(i+1));
4356  einLow = TMath::Abs(in*relativeErrorInLow->GetBinContent(1+i));
4357  out = h2->GetBinContent(i+1);
4358  eoutUp = TMath::Abs(out*relativeErrorOutUp->GetBinContent(1+i));
4359  eoutLow = TMath::Abs(out*relativeErrorOutLow->GetBinContent(1+i));
4360  // get the error squared
4361  if(rho <= 0) {
4362  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);
4363  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);
4364  } else {
4365  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);
4366  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);
4367  }
4368  if(error2Up > 0) error2Up = TMath::Sqrt(error2Up);
4369  if(error2Low > 0) error2Low = TMath::Sqrt(error2Low);
4370  // set the errors
4371  ayh[i] = error2Up;
4372  ayl[i] = error2Low;
4373  // get the bin width (which is the 'error' on x)
4374  Double_t binWidth(h1->GetBinWidth(i+1));
4375  axl[i] = binWidth/2.;
4376  axh[i] = binWidth/2.;
4377  // now get the coordinate for the poin
4378  tempV2->GetPoint(i, ax[i], ay[i]);
4379  /* if(rho<=0) {
4380  cout << "[print] pt " << ax[i] << " errorLow" << gReductionFactor*ayl[i] << " errorUp " << gReductionFactor*ayh[i] << endl;
4381  cout << "[print] %" << ax[i] << " low " << (gReductionFactor*ayl[i])/ay[i] << " up " << (gReductionFactor*ayh[i])/ay[i] << endl;
4382  } else {
4383  cout << "[print] pt " << ax[i] << " errorLow" << ayl[i]*gReductionFactorCorr << " errorUp " << ayh[i]*gReductionFactorCorr << endl;
4384  cout << "[print] %" << ax[i] << " low " << (gReductionFactorCorr*ayl[i])/ay[i] << " up " << (gReductionFactorCorr*ayh[i])/ay[i] << endl;
4385  }*/
4386  }
4387  // save the nominal ratio
4388  TGraphAsymmErrors* nominalError(new TGraphAsymmErrors(fBinsTrue->GetSize()-1, ax, ay, axl, axh, ayl, ayh));
4389  nominalError->SetName("v_{2} with shape uncertainty");
4390  // do some memory management
4391  delete tempV2;
4392  delete[] ax;
4393  delete[] ay;
4394  delete[] axh;
4395  delete[] axl;
4396  delete[] ayh;
4397  delete[] ayl;
4398 
4399  return nominalError;
4400 }
4401 //_____________________________________________________________________________
4403 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4404  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4405 #endif
4406  // write object, if a unique identifier is given the object is cloned
4407  // and the clone is saved. setting kill to true will delete the original obect from the heap
4408  if(!object) {
4409  printf(" > WriteObject:: called with NULL arguments \n ");
4410  return;
4411  } else if(!strcmp("", suffix.Data())) object->Write();
4412  else {
4413  TObject* newObject(object->Clone(Form("%s_%s", object->GetName(), suffix.Data())));
4414  newObject->Write();
4415  }
4416  if(kill) delete object;
4417 }
4418 //_____________________________________________________________________________
4420 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4421  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4422 #endif
4423  // construt a delta pt response matrix from supplied dpt distribution
4424  // binning is fine, set fBinsTrue and fBinsRec and call 'RebinTH2D' to
4425  // do a weighted rebinning to a (coarser) dpt distribution
4426  // be careful with the binning of the dpt response: it should be equal to that
4427  // of the response matrix, otherwise dpt and response matrices cannot be multiplied
4428  //
4429  // the response matrix will be square and have the same binning
4430  // (min, max and granularity) of the input histogram
4431  Int_t bins(dpt->GetXaxis()->GetNbins()); // number of bins, will also be no of rows, columns
4432  Double_t _bins[bins+1]; // prepare array with bin borders
4433  for(Int_t i(0); i < bins; i++) _bins[i] = dpt->GetBinLowEdge(i+1);
4434  _bins[bins] = dpt->GetBinLowEdge(bins)+dpt->GetBinWidth(bins+1); // get upper edge
4435  TH2D* res(new TH2D(Form("Response_from_%s", dpt->GetName()), Form("Response_from_%s", dpt->GetName()), bins, _bins, bins, _bins));
4436  for(Int_t j(0); j < bins+1; j++) { // loop on pt true slices j
4437  Bool_t skip = kFALSE;
4438  for(Int_t k(0); k < bins+1; k++) { // loop on pt gen slices k
4439  (skip) ? res->SetBinContent(j, k, 0.) : res->SetBinContent(j, k, dpt->GetBinContent(dpt->GetXaxis()->FindBin(k-j)));
4440  if(AvoidRoundingError && k > j && TMath::AreEqualAbs(dpt->GetBinContent(dpt->GetBinContent(k-j)), 0, 1e-8)) skip = kTRUE;
4441  }
4442  }
4443  return res;
4444 }
4445 //_____________________________________________________________________________
4447 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4448  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4449 #endif
4450  if(!binsTrue || !binsRec) {
4451  printf(" > GetUnityResponse:: function called with NULL arguments < \n");
4452  return 0x0;
4453  }
4454  TString name(Form("unityResponse_%s", suffix.Data()));
4455  TH2D* unity(new TH2D(name.Data(), name.Data(), binsTrue->GetSize()-1, binsTrue->GetArray(), binsRec->GetSize()-1, binsRec->GetArray()));
4456  for(Int_t i(0); i < binsTrue->GetSize(); i++) {
4457  for(Int_t j(0); j < binsRec->GetSize(); j++) {
4458  if(i==j) unity->SetBinContent(1+i, 1+j, 1.);
4459  }
4460  }
4461  return unity;
4462 }
4463 //_____________________________________________________________________________
4464 void AliJetFlowTools::SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const {
4465 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4466  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4467 #endif
4468  // save configuration parameters to histogram
4469  TH1F* summary = new TH1F("UnfoldingConfiguration","UnfoldingConfiguration", 20, -.5, 19.5);
4470  summary->SetBinContent(1, fBetaIn);
4471  summary->GetXaxis()->SetBinLabel(1, "fBetaIn");
4472  summary->SetBinContent(2, fBetaOut);
4473  summary->GetXaxis()->SetBinLabel(2, "fBetaOut");
4474  summary->SetBinContent(3, fCentralityArray->At(0));
4475  summary->GetXaxis()->SetBinLabel(3, "fCentralityArray[0]");
4476  summary->SetBinContent(4, (int)convergedIn);
4477  summary->GetXaxis()->SetBinLabel(4, "convergedIn");
4478  summary->SetBinContent(5, (int)convergedOut);
4479  summary->GetXaxis()->SetBinLabel(5, "convergedOut");
4480  summary->SetBinContent(6, (int)fAvoidRoundingError);
4481  summary->GetXaxis()->SetBinLabel(6, "fAvoidRoundingError");
4482  summary->SetBinContent(7, (int)fUnfoldingAlgorithm);
4483  summary->GetXaxis()->SetBinLabel(7, "fUnfoldingAlgorithm");
4484  summary->SetBinContent(8, (int)fPrior);
4485  summary->GetXaxis()->SetBinLabel(8, "fPrior");
4486  summary->SetBinContent(9, fSVDRegIn);
4487  summary->GetXaxis()->SetBinLabel(9, "fSVDRegIn");
4488  summary->SetBinContent(10, fSVDRegOut);
4489  summary->GetXaxis()->SetBinLabel(10, "fSVDRegOut");
4490  summary->SetBinContent(11, (int)fSVDToy);
4491  summary->GetXaxis()->SetBinLabel(11, "fSVDToy");
4492  summary->SetBinContent(12, fJetRadius);
4493  summary->GetXaxis()->SetBinLabel(12, "fJetRadius");
4494  summary->SetBinContent(13, (int)fNormalizeSpectra);
4495  summary->GetXaxis()->SetBinLabel(13, "fNormalizeSpectra");
4496  summary->SetBinContent(14, (int)fSmoothenPrior);
4497  summary->GetXaxis()->SetBinLabel(14, "fSmoothenPrior");
4498  summary->SetBinContent(15, (int)fTestMode);
4499  summary->GetXaxis()->SetBinLabel(15, "fTestMode");
4500  summary->SetBinContent(16, (int)fUseDetectorResponse);
4501  summary->GetXaxis()->SetBinLabel(16, "fUseDetectorResponse");
4502  summary->SetBinContent(17, fBayesianIterIn);
4503  summary->GetXaxis()->SetBinLabel(17, "fBayesianIterIn");
4504  summary->SetBinContent(18, fBayesianIterOut);
4505  summary->GetXaxis()->SetBinLabel(18, "fBayesianIterOut");
4506  summary->SetBinContent(19, fBayesianSmoothIn);
4507  summary->GetXaxis()->SetBinLabel(19, "fBayesianSmoothIn");
4508  summary->SetBinContent(20, fBayesianSmoothOut);
4509  summary->GetXaxis()->SetBinLabel(20, "fBayesianSmoothOut");
4510 }
4511 //_____________________________________________________________________________
4513 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4514  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4515 #endif
4516  // ugly function: reset all unfolding parameters
4517  TVirtualFitter* fitter(TVirtualFitter::GetFitter());
4518  if(fitter) {
4519  printf(" > Found fitter, will delete it < \n");
4520  delete fitter;
4521  }
4522  if(gMinuit) {
4523  printf(" > Found gMinuit, will re-create it < \n");
4524  delete gMinuit;
4525  gMinuit = new TMinuit;
4526  }
4527  AliUnfolding::fgCorrelationMatrix = 0;
4528  AliUnfolding::fgCorrelationMatrixSquared = 0;
4529  AliUnfolding::fgCorrelationCovarianceMatrix = 0;
4530  AliUnfolding::fgCurrentESDVector = 0;
4531  AliUnfolding::fgEntropyAPriori = 0;
4532  AliUnfolding::fgEfficiency = 0;
4533  AliUnfolding::fgUnfoldedAxis = 0;
4534  AliUnfolding::fgMeasuredAxis = 0;
4535  AliUnfolding::fgFitFunction = 0;
4536  AliUnfolding::fgMaxInput = -1;
4537  AliUnfolding::fgMaxParams = -1;
4538  AliUnfolding::fgOverflowBinLimit = -1;
4539  AliUnfolding::fgRegularizationWeight = 10000;
4540  AliUnfolding::fgSkipBinsBegin = 0;
4541  AliUnfolding::fgMinuitStepSize = 0.1;
4542  AliUnfolding::fgMinuitPrecision = 1e-6;
4543  AliUnfolding::fgMinuitMaxIterations = 1000000;
4544  AliUnfolding::fgMinuitStrategy = 1.;
4545  AliUnfolding::fgMinimumInitialValue = kFALSE;
4546  AliUnfolding::fgMinimumInitialValueFix = -1;
4547  AliUnfolding::fgNormalizeInput = kFALSE;
4548  AliUnfolding::fgNotFoundEvents = 0;
4549  AliUnfolding::fgSkipBin0InChi2 = kFALSE;
4550  AliUnfolding::fgBayesianSmoothing = 1;
4551  AliUnfolding::fgBayesianIterations = 10;
4552  AliUnfolding::fgDebug = kFALSE;
4553  AliUnfolding::fgCallCount = 0;
4554  AliUnfolding::fgPowern = 5;
4555  AliUnfolding::fChi2FromFit = 0.;
4556  AliUnfolding::fPenaltyVal = 0.;
4557  AliUnfolding::fAvgResidual = 0.;
4558  AliUnfolding::fgPrintChi2Details = 0;
4559  AliUnfolding::fgCanvas = 0;
4560  AliUnfolding::fghUnfolded = 0;
4561  AliUnfolding::fghCorrelation = 0;
4562  AliUnfolding::fghEfficiency = 0;
4563  AliUnfolding::fghMeasured = 0;
4564  AliUnfolding::SetMinuitStepSize(1.);
4565  AliUnfolding::SetMinuitPrecision(1e-6);
4566  AliUnfolding::SetMinuitMaxIterations(100000);
4567  AliUnfolding::SetMinuitStrategy(2.);
4568  AliUnfolding::SetDebug(1);
4569 }
4570 //_____________________________________________________________________________
4572 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4573  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4574 #endif
4575  // protect heap by adding unique qualifier to name
4576  if(!protect) return 0x0;
4577  TH1D* p = (TH1D*)protect->Clone();
4578  TString tempString(fActiveString);
4579  tempString+=suffix;
4580  p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
4581  p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4582  if(kill) delete protect;
4583  return p;
4584 }
4585 //_____________________________________________________________________________
4587 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4588  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4589 #endif
4590  // protect heap by adding unique qualifier to name
4591  if(!protect) return 0x0;
4592  TH2D* p = (TH2D*)protect->Clone();
4593  TString tempString(fActiveString);
4594  tempString+=suffix;
4595  p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
4596  p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4597  if(kill) delete protect;
4598  return p;
4599 }
4600 //_____________________________________________________________________________
4602 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4603  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4604 #endif
4605  // protect heap by adding unique qualifier to name
4606  if(!protect) return 0x0;
4607  TGraphErrors* p = (TGraphErrors*)protect->Clone();
4608  TString tempString(fActiveString);
4609  tempString+=suffix;
4610  p->SetName(Form("%s_%s", protect->GetName(), tempString.Data()));
4611  p->SetTitle(Form("%s_%s", protect->GetTitle(), tempString.Data()));
4612  if(kill) delete protect;
4613  return p;
4614 }
4615 //_____________________________________________________________________________
4617 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4618  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4619 #endif
4620  // === azimuthal unfolding ===
4621  //
4622  // unfolds the spectrum in delta phi bins, extracts the yield per bin, and does a fit
4623  // in transverse momentum and azimuthal correlation space to extract v2 params
4624  // settings are equal to the ones used for 'Make()'
4625  //
4626  // basic steps that are followed:
4627  // 1) rebin the raw output of the jet task to the desired binnings
4628  // 2) calls the unfolding routine
4629  // 3) writes output to file
4630  // can be repeated multiple times with different configurations
4631 
4632  Int_t low[] = {1, 6, 11, 16, 21, 26, 31, 36};
4633  Int_t up[] = {5, 10, 15, 20, 25, 30, 35, 40};
4634  TString stringArray[] = {"a", "b", "c", "d", "e", "f", "g", "h"};
4635  const Int_t ptBins(fBinsTrue->GetSize()-1);
4636  const Int_t dPhiBins(8);
4637  TH1D* dPtdPhi[fBinsTrue->GetSize()];
4638  for(Int_t i(0); i < ptBins; i++) dPtdPhi[i] = new TH1D(Form("dPtdPhi_%i", i), Form("dPtdPhi_%i", i), dPhiBins, 0, TMath::Pi());
4639 
4640  for(Int_t i(0); i < dPhiBins; i++) {
4641  // 1) manipulation of input histograms
4642  // check if the input variables are present
4643  if(!PrepareForUnfolding(low[i], up[i])) return;
4644  // 1a) resize the jet spectrum according to the binning scheme in fBinsTrue
4645  // parts of the spectrum can end up in over or underflow bins
4646  TH1D* measuredJetSpectrumIn = RebinTH1D(fSpectrumIn, fBinsRec, Form("resized_%s", stringArray[i].Data()), kFALSE);
4647 
4648  // 1b) resize the jet spectrum to 'true' bins. can serve as a prior and as a template for unfolding
4649  // the template will be used as a prior for the chi2 unfolding
4650  TH1D* measuredJetSpectrumTrueBinsIn = RebinTH1D(fSpectrumIn, fBinsTrue, stringArray[i], kFALSE);
4651 
4652  // get the full response matrix from the dpt and the detector response
4654  // get the full response matrix. if test mode is chosen, the full response is replace by a unity matrix
4655  // so that unfolding should return the initial spectrum
4660  // normalize each slide of the response to one
4662  // resize to desired binning scheme
4663  TH2D* resizedResponseIn = RebinTH2D(fFullResponseIn, fBinsTrue, fBinsRec, stringArray[i]);
4664  // get the kinematic efficiency
4665  TH1D* kinematicEfficiencyIn = resizedResponseIn->ProjectionX();
4666  kinematicEfficiencyIn->SetNameTitle(Form("kin_eff_%s", stringArray[i].Data()), Form("kin_eff_%s", stringArray[i].Data()));
4667  // suppress the errors
4668  for(Int_t j(0); j < kinematicEfficiencyIn->GetXaxis()->GetNbins(); j++) kinematicEfficiencyIn->SetBinError(1+j, 0.);
4669  TH1D* jetFindingEfficiency(0x0);
4670  if(fJetFindingEff) {
4671  jetFindingEfficiency = ProtectHeap(fJetFindingEff);
4672  jetFindingEfficiency->SetNameTitle(Form("%s_coarse", jetFindingEfficiency->GetName()), Form("%s_coarse", jetFindingEfficiency->GetName()));
4673  jetFindingEfficiency = RebinTH1D(jetFindingEfficiency, fBinsTrue);
4674  }
4675  // 2, 3) call the actual unfolding. results and transient objects are stored in a dedicated TDirectoryFile
4676  TH1D* unfoldedJetSpectrumIn(0x0);
4677  fActiveDir->cd(); // select active dir
4678  TDirectoryFile* dirIn = new TDirectoryFile(Form("%s___%s", stringArray[i].Data(), fActiveString.Data()), Form("%s___%s", stringArray[i].Data(), fActiveString.Data()));
4679  dirIn->cd(); // select inplane subdir
4680  // select the unfolding method
4681  unfoldedJetSpectrumIn = UnfoldWrapper(
4682  measuredJetSpectrumIn,
4683  resizedResponseIn,
4684  kinematicEfficiencyIn,
4685  measuredJetSpectrumTrueBinsIn,
4686  TString("dPtdPhiUnfolding"),
4687  jetFindingEfficiency);
4688  // arbitrarily save one of the full outputs (same for all dphi bins, avoid duplicates)
4689  if(i+1 == ptBins) {
4690  resizedResponseIn->SetNameTitle(Form("ResponseMatrix_%s", stringArray[i].Data()), Form("response matrix %s", stringArray[i].Data()));
4691  resizedResponseIn->SetXTitle("p_{T, jet}^{true} [GeV/c]");
4692  resizedResponseIn->SetYTitle("p_{T, jet}^{rec} [GeV/c]");
4693  resizedResponseIn = ProtectHeap(resizedResponseIn);
4694  resizedResponseIn->Write();
4695  kinematicEfficiencyIn->SetNameTitle(Form("KinematicEfficiency_%s", stringArray[i].Data()), Form("Kinematic efficiency, %s", stringArray[i].Data()));
4696  kinematicEfficiencyIn = ProtectHeap(kinematicEfficiencyIn);
4697  kinematicEfficiencyIn->Write();
4698  fDetectorResponse->SetNameTitle("DetectorResponse", "Detector response matrix");
4700  fDetectorResponse->Write();
4701  // optional histograms
4702  if(fSaveFull) {
4703  fSpectrumIn->SetNameTitle("[ORIG]JetSpectrum", Form("[INPUT] Jet spectrum, %s", stringArray[i].Data()));
4704  fSpectrumIn->Write();
4705  fDptInDist->SetNameTitle("[ORIG]DeltaPt", Form("#delta p_{T} distribution, %s", stringArray[i].Data()));
4706  fDptInDist->Write();
4707  fDptIn->SetNameTitle("[ORIG]DeltaPtMatrix", Form("#delta p_{T} matrix, %s", stringArray[i].Data()));
4708  fDptIn->Write();
4709  fFullResponseIn->SetNameTitle("ResponseMatrix", Form("Response matrix, %s", stringArray[i].Data()));
4710  fFullResponseIn->Write();
4711  }
4712  }
4713  fActiveDir->cd();
4714  fDeltaPtDeltaPhi->Write();
4715  fJetPtDeltaPhi->Write();
4716 
4717  TH1D* dud(ProtectHeap(unfoldedJetSpectrumIn, kTRUE, stringArray[i]));;
4718  Double_t integralError(0);
4719  // at this point in the code, the spectrum has been unfolded in a certain region of dPhi space
4720  // next step is splitting it in pt space as well to estimate the yield differentially in pt
4721  for(Int_t j(0); j < ptBins; j++) {
4722  // get the integrated
4723  Double_t integral(dud->IntegralAndError(j+1, j+2, integralError));
4724  dPtdPhi[j]->SetBinContent(i+1, integral);
4725  dPtdPhi[j]->SetBinError(i+1, integralError);
4726  }
4727  dud->Write();
4728  // save the current state of the unfolding object
4729  SaveConfiguration(unfoldedJetSpectrumIn ? kTRUE : kFALSE, kFALSE);
4730  }
4731  TF1* fourier = new TF1("fourier", "[0]*(1.+0.5*[1]*(TMath::Cos(2.*x)))", 0, TMath::Pi());
4732  TH1D* v2(new TH1D("v2FromFit", "v2FromFit", fBinsTrue->GetSize()-1, fBinsTrue->GetArray()));
4733  for(Int_t i(0); i < ptBins; i++) {
4734  dPtdPhi[i]->Fit(fourier, "VI");
4735  Style(gPad, "PEARSON");
4736  Style(dPtdPhi[i], kBlue, kDeltaPhi);
4737  dPtdPhi[i]->DrawCopy();
4739  TString(Form("%.2f #LT p_{T} #LT %.2f", fBinsTrue->At(i), fBinsTrue->At(i+1))),
4740  kBlack,
4741  .38,
4742  .56,
4743  .62,
4744  .65
4745  );
4746  v2->SetBinContent(1+i, fourier->GetParameter(1));
4747  v2->SetBinError(1+i, fourier->GetParError(1));
4748  Style(gPad, "PEARSON");
4749  Style(v2, kBlack, kV2, kTRUE);
4750  v2->DrawCopy();
4751  dPtdPhi[i]->Write();
4752  }
4753 }
4754 //_____________________________________________________________________________
4756  // replace bins
4757 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4758  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4759 #endif
4760  Double_t x(0), y(0);
4761  graph->GetPoint(0, x, y);
4762  graph->SetPoint(array->At(0)-1, fBinsTrue->At(array->At(0)), y);
4763  graph->SetPointError(array->At(0)-1, 10, graph->GetErrorY(0));
4764  graph->SetPoint(array->At(1)-1, -5, -5);
4765 }
4766 //_____________________________________________________________________________
4768 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4769  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4770 #endif
4771  // replace bins
4772  Double_t x(0), y(0);
4773  graph->GetPoint(0, x, y);
4774  graph->SetPoint(array->At(0)-1, fBinsTrue->At(array->At(0)), y);
4775  Double_t yl = graph->GetErrorYlow(0);
4776  Double_t yh = graph->GetErrorYhigh(0);
4777  graph->SetPointError(array->At(0)-1, 10, 10, yl, yh);
4778  graph->SetPoint(array->At(1)-1, -5, -5);
4779 }
4780 //_____________________________________________________________________________
4782  TGraphErrors* n, // points with stat error
4783  TGraphAsymmErrors* shape, // points with shape error
4784  TGraphAsymmErrors* corr, // points with stat error
4785  Int_t low, // lower bin (tgraph starts at 0)
4786  Int_t up // upper bin
4787  )
4788 {
4789 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4790  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4791 #endif
4792  // main use of this function is filling the static buffers
4793  Double_t statE(0), shapeE(0), corrE(0), y(0), x(0);
4794 
4795  // print some stuff
4796  printf(" double v2[] = {\n");
4797  Int_t iterator(0);
4798  for(Int_t i(low); i < up+1; i++) {
4799  n->GetPoint(i, x, y);
4800  if(i==up) printf("%.4f}; \n\n", y);
4801  else printf("%.4f, \n", y);
4802  gV2->SetAt(y, iterator);
4803  iterator++;
4804  }
4805  iterator = 0;
4806  printf(" double stat[] = {\n");
4807  for(Int_t i(low); i < up+1; i++) {
4808  y = n->GetErrorYlow(i);
4809  if(i==up) printf("%.4f}; \n\n", y);
4810  else printf("%.4f, \n", y);
4811  gStat->SetAt(y, iterator);
4812  iterator++;
4813  }
4814  iterator = 0;
4815  printf(" double shape[] = {\n");
4816  for(Int_t i(low); i < up+1; i++) {
4817  y = shape->GetErrorYhigh(i);
4818  if (y < 0.0001) y = 0.0001; // avoid crash in fit routine
4819  y*=gReductionFactor;
4820  shape->SetPointEYhigh(i, y);
4821  y = shape->GetErrorYlow(i);
4822  if ( y < 0.0001) y = 0.0001;
4823  y*=gReductionFactor;
4824  shape->SetPointEYlow(i, y);
4825  if(i==up) printf("%.4f}; \n\n", y);
4826  else printf("%.4f, \n", y);
4827  gShape->SetAt(y, iterator);
4828  iterator++;
4829  }
4830  iterator = 0;
4831  printf(" double corr[] = {\n");
4832  for(Int_t i(low); i < up+1; i++) {
4833  y = gReductionFactorCorr*corr->GetErrorYhigh(i);
4834  corr->SetPointEYhigh(i, y);
4835  y = gReductionFactorCorr*corr->GetErrorYlow(i);
4836  corr->SetPointEYlow(i, y);
4837  if(i==up) printf("%.4f}; \n\n", y);
4838  else printf("%.4f, \n", y);
4839  gCorr->SetAt(y, iterator);
4840  iterator++;
4841  }
4842 
4843  // to plot the average error as function of number of events
4844  Float_t ctr(0);
4845  for(Int_t i(low); i < up+1; i++) {
4846  ctr = ctr + 1.;
4847  // set some flags to 0
4848  x = 0.;
4849  y = 0.;
4850  // get the nominal point
4851  n->GetPoint(i, x, y);
4852  statE += n->GetErrorY(i);
4853  shapeE += shape->GetErrorY(i);
4854  corrE += corr->GetErrorY(i);
4855  // combine the errors
4856  }
4857  printf(" ======================================\n");
4858  printf(" > between %i and %i GeV/c \n", low, up);
4859  cout << " AVERAGE_SHAPE " << shapeE/ctr << endl;
4860  cout << " AVERAGE_CORR " << corrE/ctr << endl;
4861  cout << " AVERAGE_STAT " << statE/ctr << endl;
4862 }
4863 //_____________________________________________________________________________
4865 {
4866 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4867  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4868 #endif
4869  // Choose method upon creation between:
4870  // kMigrad, kSimplex, kCombined,
4871  // kScan, kFumili
4872  ROOT::Minuit2::Minuit2Minimizer min ( ROOT::Minuit2::kMigrad );
4873  min.SetMaxFunctionCalls(1000000);
4874  min.SetMaxIterations(100000);
4875  min.SetTolerance(0.001);
4876 
4877  ROOT::Math::Functor f(&PhenixChi2nd,2);
4878  double step[] = {0.0000001, 0.0000001};
4879  double variable[] = {-1., -1.};
4880 
4881  min.SetFunction(f);
4882  // Set the free variables to be minimized!
4883  min.SetVariable(0,"epsilon_c",variable[0], step[0]);
4884  min.SetVariable(1,"epsilon_b",variable[1], step[1]);
4885 
4886 
4887  min.Minimize();
4888 }
4889 //_____________________________________________________________________________
4891 {
4892 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4893  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4894 #endif
4895  // define arrays with results and errors here, see example at PhenixChi2()
4896 
4897  // return the function value at certain epsilon
4898  const Double_t epsc = xx[0];
4899  const Double_t epsb = xx[1];
4900  Double_t chi2(0);
4901  Int_t counts(gV2->GetSize() + gOffsetStop);
4902 
4903  if(gPwrtTo > -999) {
4904  // altered implemtation of eq 3 of arXiv:0801.1665v2
4905  // see analysis note and QM2014 poster for validation
4906  for(Int_t i(gOffsetStart); i < counts; i++) {
4907 
4908  // quadratic sum of statistical and uncorrelated systematic error
4909  Double_t e = gStat->At(i);
4910 
4911  // sum of v2 plus epsilon times correlated error minus hypothesis (gPwrtTo)
4912  // also the numerator of equation 3 of phenix paper
4913  Double_t numerator = TMath::Power(gV2->At(i)+epsc*gCorr->At(i)+epsb-gPwrtTo, 2);
4914 
4915  // modified denominator of equation 3 of phenix paper
4916  Double_t denominator = e*e;
4917 
4918  // add to the sum
4919  chi2 += numerator/denominator;
4920  }
4921  // add the square of epsilon to the total chi2 as penalty
4922 
4923  Double_t sumEpsb(0);
4924  for(Int_t j(gOffsetStart); j < counts; j++) sumEpsb += (epsb*epsb)/(gShape->At(j)*gShape->At(j));
4925  chi2 += epsc*epsc + sumEpsb/((float)counts);
4926  } else {
4927  // otherwise use the array
4928  for(Int_t i(gOffsetStart); i < counts; i++) {
4929 
4930  // quadratic sum of statistical and uncorrelated systematic error
4931  Double_t e = TMath::Sqrt(gStat->At(i)*gStat->At(i) + gPwrtToStatArray->At(i)*gPwrtToStatArray->At(i));
4932 
4933  // sum of v2 plus epsilon times correlated error minus hypothesis (gPwrtTo)
4934  // also the numerator of equation 3 of phenix paper
4935  Double_t numerator = TMath::Power(gV2->At(i)+epsc*gCorr->At(i)+epsb-gPwrtToArray->At(i), 2);
4936 
4937  // modified denominator of equation 3 of phenix paper
4938  Double_t denominator = e*e;
4939 
4940  // add to the sum
4941  chi2 += numerator/denominator;
4942  }
4943  // add the square of epsilon to the total chi2 as penalty
4944 
4945  Double_t sumEpsb(0);
4946  for(Int_t j(gOffsetStart); j < counts; j++) sumEpsb += (epsb*epsb)/(gShape->At(j)*gShape->At(j));
4947  chi2 += epsc*epsc + sumEpsb/((float)counts);
4948  }
4949  return chi2;
4950 }
4951 //_____________________________________________________________________________
4953 {
4954 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4955  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4956 #endif
4957  // internal use only: evaluate the function at a given point
4958  if(par) SquelchWarning();
4960 }
4961 //_____________________________________________________________________________
4963 {
4964 #ifdef ALIJETFLOWTOOLS_DEBUG_FLAG
4965  printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
4966 #endif
4967  // return the fitting function, pass the p-value w.r.t. 0 by reference
4968  const Int_t DOF(7-2); // dof is n-2
4969  TF2 *f1 = new TF2("ndhist", AliJetFlowTools::ConstructFunctionnd, -100, 100, -100, 100, 0);
4970  printf(" > locating minima < \n");
4971  Double_t x(0), y(0);
4972  f1->GetMinimumXY(x, y);
4973  f1->GetXaxis()->SetTitle("#epsilon{c}");
4974  f1->GetXaxis()->SetTitle("#epsilon_{b}");
4975  f1->GetZaxis()->SetTitle("#chi^{2}");
4976 
4977  printf(" ===============================================================================\n");
4978  printf(" > minimal chi2 f(%.8f, %.8f) = %.8f (i should be ok ... ) \n", x, y, f1->Eval(x, y));
4979  cout << " so the probability of finding data at least as imcompatible with " << gPwrtTo << " as the actually" << endl;
4980  cout << " observed data is " << TMath::Prob(f1->Eval(x, y), DOF) << endl;
4981  cout << " minimization parameters: EPSILON_B " << y << endl;
4982  cout << " EPSILON_C " << x << endl;
4983  printf(" ===============================================================================\n");
4984 
4985  // pass the p-value by reference and return the function
4986  p = TMath::Prob(f1->Eval(x, y), DOF);
4987  return f1;
4988 }
4989 //_____________________________________________________________________________
TH1D * UnfoldSpectrumSVD(const TH1D *measuredJetSpectrum, const TH2D *resizedResponse, const TH1D *kinematicEfficiency, const TH1D *measuredJetSpectrumTrueBins, const TString suffix, const TH1D *jetFindingEfficiency=0x0)
return jsonbuilder str().c_str()
static TH1D * SmoothenPrior(TH1D *spectrum, TF1 *function, Double_t min, Double_t max, Double_t start, Bool_t kill=kTRUE, Bool_t counts=kTRUE)
static void WriteObject(TObject *object, TString suffix="", Bool_t kill=kTRUE)
double Double_t
Definition: External.C:58
static TH2D * ResizeYaxisTH2D(TH2D *histo, TArrayD *x, TArrayD *y, TString suffix="")
TH1D * UnfoldSpectrumBayesian(const TH1D *measuredJetSpectrum, const TH2D *resizedResponse, const TH1D *kinematicEfficiency, const TH1D *measuredJetSpectrumTrueBins, const TString suffix, const TH1D *jetFindingEfficiency=0x0)
Bool_t fUseDetectorResponse
static TH2D * NormalizeTH2D(TH2D *histo, Bool_t noError=kTRUE)
static Float_t gPwrtTo
TArrayD * fBinsRecPrior
static Int_t gOffsetStart
TDirectoryFile * fActiveDir
static Double_t ConstructFunctionnd(Double_t *x, Double_t *par)
TH1D * UnfoldSpectrumBayesianAli(const TH1D *measuredJetSpectrum, const TH2D *resizedResponse, const TH1D *kinematicEfficiency, const TH1D *measuredJetSpectrumTrueBins, const TString suffix, const TH1D *jetFindingEfficiency=0x0)
static TGraphAsymmErrors * AddHistoErrorsToGraphErrors(TGraphAsymmErrors *g, TH1D *h)
static Double_t PhenixChi2nd(const Double_t *xx)
AliAnaChargedJetResponseMaker * fResponseMaker
void GetNominalValues(TH1D *&ratio, TGraphErrors *&v2, TArrayI *in, TArrayI *out, TString inFile="UnfoldedSpectra.root", TString outFile="Nominal.root") const
TCanvas * c
Definition: TestFitELoss.C:172
Float_t fBayesianSmoothOut
static TArrayD * gStat