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